...
 
Commits (6)
......@@ -2,6 +2,11 @@
!/**/
!*.*
makefile
src/hplog.f
src/tdhpl.f
adpg_v1_0.gz
adpu_v1_0.tar.gz
# gcc coverage testing tool files
*.gcno
......
......@@ -34,6 +34,14 @@ PROGRAM TEST
case('-gpl-test')
tol = zero * 1.e5_prec
call do_GPL_tests
#ifdef HAVE_TDHPL
case('-hpl-test')
tol = 10*zero
call do_HPL_tests
#else
case('-hpl-test')
call iprint("Argument "//trim(arg)//" is not available, compile with --2dhpl", 2)
#endif
case('-chen-test')
tol = 8.0e-7
tests_successful = tests_successful .and. do_chen_test()
......@@ -468,6 +476,88 @@ CONTAINS
end subroutine do_GPL_tests
#ifdef HAVE_TDHPL
subroutine do_HPL_tests
use HPL_interface, only: hpl
implicit none
complex(kind=prec) :: ref
ref = (0.26236426446749105203549598688095440_prec,0.)
call test_one_HPL((/-1/), 0.3_prec, ref, "H.1")
ref = (-1.2039728043259359926227462177618385_prec,0.)
call test_one_HPL((/0/), 0.3_prec, ref, "H.2")
ref = (0.35667494393873237891263871124118448_prec,0.)
call test_one_HPL((/1/), 0.3_prec, ref, "H.3")
ref = (0.28007433375958290423021697230542536_prec,0.)
call test_one_HPL((/0,-1/), 0.3_prec, ref, "H.4")
ref = (-0.59595377300541962996225138433054502_prec,0.)
call test_one_HPL((/-1,0/), 0.3_prec, ref, "H.5")
ref = (0.72477525677822927901564821362260275_prec,0.)
call test_one_HPL((/0,0/), 0.3_prec, ref, "H.6")
ref = (-0.75555644256218769787140889183866340_prec,0.)
call test_one_HPL((/1,0/), 0.3_prec, ref, "H.7")
ref = (-0.91648194931990515581968426452348372_prec,0.)
call test_one_HPL((/0,-1,0/), 0.3_prec, ref, "H.8")
ref = (-1.0174514166047997621715595600723809_prec,0.)
call test_one_HPL((/0,1,0/), 0.3_prec, ref, "H.9")
ref = (-0.096156176800657414704535610841285475_prec,0.)
call test_one_HPL((/-1,-1,0/), 0.3_prec, ref, "H.10")
ref = (0.81699704232693125973188535221842666_prec,0.)
call test_one_HPL((/-1,0,0/), 0.3_prec, ref, "H.11")
ref = (-0.29086989946977835131926209119273790_prec,0.)
call test_one_HPL((/0,0,0/), 0.3_prec, ref, "H.12")
ref = (0.96356041279146243412056107045850765_prec,0.)
call test_one_HPL((/1,0,0/), 0.3_prec, ref, "H.13")
ref = (-1.2327590115110687817725163762469362_prec,0.)
call test_one_HPL((/0,0,-1,0/), 0.3_prec, ref, "H.14")
ref = (-1.2941049241728956868728872939899452_prec,0.)
call test_one_HPL((/0,0,1,0/), 0.3_prec, ref, "H.15")
ref = (-0.061811027943749597866741119227363730_prec,0.)
call test_one_HPL((/0,-1,-1,0/), 0.3_prec, ref, "H.16")
ref = (1.7844686828294620604418638073704301_prec,0.)
call test_one_HPL((/0,-1,0,0/), 0.3_prec, ref, "H.17")
ref = (-0.13615497727913000601405538432901756_prec,0.)
call test_one_HPL((/-1,0,-1,0/), 0.3_prec, ref, "H.18")
ref = (-0.20093319914956897041204850097815064_prec,0.)
call test_one_HPL((/1,0,-1,0/), 0.3_prec, ref, "H.19")
ref = (-0.21577152810861711268331714964024090_prec,0.)
call test_one_HPL((/1,0,1,0/), 0.3_prec, ref, "H.20")
ref = (1.9065968418304341703412524922325282_prec,0.)
call test_one_HPL((/0,1,0,0/), 0.3_prec, ref, "H.21")
ref = (0.15686771352941380971916881529487140_prec,0.)
call test_one_HPL((/-1,-1,0,0/), 0.3_prec, ref, "H.22")
ref = (-0.92270363433527097196573008319186789_prec,0.)
call test_one_HPL((/-1,0,0,0/), 0.3_prec, ref, "H.23")
ref = (0.087549862139658031076027422464290803_prec,0.)
call test_one_HPL((/0,0,0,0/), 0.3_prec, ref, "H.24")
ref = (-1.0222324580521425613350726928441460_prec,0.)
call test_one_HPL((/1,0,0,0/), 0.3_prec, ref, "H.25")
end subroutine do_HPL_tests
#endif
#ifdef HAVE_GINAC
......
......@@ -119,4 +119,24 @@ contains
endif
end subroutine test_one_flat
#ifdef HAVE_TDHPL
subroutine test_one_HPL(w,x,ref,test_id)
use hpl_interface, only: HPL
integer :: w(:)
real(kind=prec) :: x
complex(kind=prec) :: res, ref
character(len=*) :: test_id
call iprint(' testing HPL '//test_id//' against ref ...',-1)
res = HPL(w, x)
call check(res,ref)
call iprint(' testing HPL '//test_id//' against handyG ...',-1)
ref = (-1)**count(w.eq.1) * G([real(w,kind=prec), x])
call check(res,ref)
end subroutine test_one_HPL
#endif
END MODULE TTOOLS
......@@ -111,6 +111,8 @@ for arg in "$@" ; do
CONF_PREFIX="${arg#--prefix=}" ;;
--moduledir=*)
CONF_MODDIR="${arg#--moduledir=}" ;;
--with-2dhpl)
HAVE_2DHPL=true ;;
--with-ginac)
HAVE_GINAC=true ;;
--with-mcc)
......@@ -701,6 +703,16 @@ _EOF_
fi
fi
if $HAVE_2DHPL ; then
eval addflag FFLAGS -DHAVE_TDHPL
echo -n "do we have 2DHPL... " 1>&3
test -f "adpu_v1_0.tar.gz" && test -f "adpg_v1_0.gz" && \
echo " yes" 1>&3 || {
echo " no. Please donwload adpu and adpg from CPC." 1>&3 ;
exit 1
}
fi
if $HAVE_GINAC && $HAVE_MCC; then
eval addflag LFLAGS -lrt
......@@ -770,8 +782,11 @@ PREFIX=$CONF_PREFIX
MMPREFIX=$CONF_MMPREFIX
MODDIR=$CONF_MODDIR
files=globals.o ieps.o utils.o shuffle.o maths_functions.o mpl_module.o gpl_module.o handyG.o
EOF
echo -n "files=globals.o ieps.o" >> makefile
$HAVE_2DHPL && echo -n " tdhpl.o hplog.o hpl.o" >> makefile
echo " utils.o shuffle.o maths_functions.o mpl_module.o gpl_module.o handyG.o" >> makefile
cat >> makefile <<EOF
objects = \$(addprefix build/,\$(files))
......@@ -804,6 +819,21 @@ geval: libhandyg.a build/geval.o
EOF
$HAVE_2DHPL && cat >> makefile <<EOF
src/hplog.f: adpg_v1_0.gz
@echo "GZ \$<"
@cat \$< | gzip -d | tail -n +4 > \$@
src/tdhpl.f: adpu_v1_0.tar.gz
@echo "TAR x \$<"
@tar xzfO \$< tdhpl.f > \$@
build/%.o: src/%.f
@echo "F \$@"
@\$(FC) -c \$< -o \$@
EOF
$HAVE_MCC && cat >> makefile <<EOF
build/%.tm.c: src/%.tm
@echo "MPREP \$@"
......@@ -902,6 +932,10 @@ check-chen: test
./$< -chen-test
EOF
$HAVE_2DHPL && cat >> makefile <<EOF
check-hpl: test
./$< -hpl-test
EOF
$HAVE_GINAC && $HAVE_MCC && cat >> makefile <<EOF
check-ginac: test
./$< -ginac-tests
......
MODULE HPL_INTERFACE
use globals
type el
real(kind=prec) :: x
integer nw
complex(kind=8), dimension(4,-1:1,-1:1,-1:1,-1:1) :: mat
end type el
type(el) :: cache(100)
integer :: cachesize = 0
CONTAINS
FUNCTION SELECTANS(w, mat)
integer w(:)
complex(kind=8) :: mat(4,-1:1,-1:1,-1:1,-1:1), selectans
select case(size(w))
case(1)
selectans = mat(1,0,0,0,w(1))
case(2)
selectans = mat(2,0,0,w(1),w(2))
case(3)
selectans = mat(3,0,w(1),w(2),w(3))
case(4)
selectans = mat(4,w(1),w(2),w(3),w(4))
end select
END FUNCTION
FUNCTION HPL(w, x, inw)
integer :: w(:)
real(kind=prec) :: x
integer nw
integer, optional :: inw
complex(kind=prec) :: hpl
complex(kind=8) :: mat(4,-1:1,-1:1,-1:1,-1:1)
complex(kind=8) hc1(-1:1)
complex(kind=8) hc2(-1:1,-1:1)
complex(kind=8) hc3(-1:1,-1:1,-1:1)
complex(kind=8) hc4(-1:1,-1:1,-1:1,-1:1)
real(kind=8) hr1(-1:1)
real(kind=8) hr2(-1:1,-1:1)
real(kind=8) hr3(-1:1,-1:1,-1:1)
real(kind=8) hr4(-1:1,-1:1,-1:1,-1:1)
real(kind=8) hi1(-1:1)
real(kind=8) hi2(-1:1,-1:1)
real(kind=8) hi3(-1:1,-1:1,-1:1)
real(kind=8) hi4(-1:1,-1:1,-1:1,-1:1)
nw = size(w)
if (present(inw)) then
nw = max(inw,size(w))
endif
if (nw.gt.4) then
print*,"Can't calculate HPL with w>4"
stop
endif
#ifndef NOCACHE
do i=1,cachesize
if ((nw.le.cache(i)%nw).and.(abs(x-cache(i)%x) .lt. zero)) then
hpl = cmplx(selectans(w, cache(i)%mat),kind=prec)
return
endif
enddo
#endif
call hplog(x, nw, hc1,hc2,hc3,hc4, &
hr1,hr2,hr3,hr4, hi1, hi2, hi3, hi4, -1, 1)
mat = 0.
if (nw>=1) mat(1,0,0,0,:) = hc1
if (nw>=2) mat(2,0,0,:,:) = hc2
if (nw>=3) mat(3,0,:,:,:) = hc3
if (nw>=4) mat(4,:,:,:,:) = hc4
#ifndef NOCACHE
if (cachesize.lt.size(cache)) then
cachesize = cachesize + 1
cache(cachesize) = el(x,nw,mat)
endif
#endif
hpl = cmplx(selectans(w, mat),kind=prec)
END FUNCTION HPL
END MODULE HPL_INTERFACE