Commit 5c8d3f84 authored by ulrich_y's avatar ulrich_y
Browse files

Tests are now precise enough for quad precision

This follows up on 9dbf60e8
parent 9dbf60e8
This diff is collapsed.
......@@ -11,8 +11,7 @@ PROGRAM TEST
use chenreftest, only: do_chen_test
implicit none
complex(kind=prec) :: res
real :: tol = 8.0e-7
real(kind=prec) :: tol = 8.0e-7
logical :: tests_successful = .true.
#ifdef DEBUG
......@@ -24,7 +23,7 @@ PROGRAM TEST
call do_GPL_tests()
! call do_shuffle_tests() ! put this somewhere else
tol = 8.0e-7
tests_successful = tests_successful .and. do_chen_test(cmplx(0.3),cmplx(0.1))
tests_successful = tests_successful .and. do_chen_test()
#ifdef HAVE_GINAC
#ifdef HAVE_MM
......@@ -57,7 +56,7 @@ CONTAINS
subroutine test_one_MPL(m,x,ref, test_id)
integer :: m(:)
complex(kind=prec) :: x(:), ref
complex(kind=prec) :: x(:), ref, res
character(len=*) :: test_id
print*, ' ', 'testing MPL ', test_id, ' ...'
......@@ -69,19 +68,19 @@ CONTAINS
complex(kind=prec) :: ref
print*, 'doing MPL tests...'
ref = cmplx(0.022696600480693277651633)
call test_one_MPL((/ 1,1 /),cmplx((/ 0.3156498673740053, 0.3431255827785649 /)),ref, '1.1')
ref = (0.022696600480693277651633_prec,0.)
call test_one_MPL((/ 1,1 /),(/ (0.3156498673740053_prec,0.), (0.3431255827785649_prec,0.) /),ref, '1.1')
ref = cmplx(0.00023134615630308335448329926098409)
call test_one_MPL((/ 1,1 /),cmplx((/ 0.03, 0.5012562893380046 /)),ref, '1.2')
ref = (0.00023134615630308335448329926098409_prec,0.)
call test_one_MPL((/ 1,1 /),(/ (0.03_prec,0.), (0.5012562893380046_prec,0.) /),ref, '1.2')
ref = cmplx(0.000023446106415452030937059124671151)
call test_one_MPL((/ 2,1,2 /),cmplx((/ 0.03, 0.5012562893380046, 55.3832 /)),ref, '1.3')
ref = (0.000023446106415452030937059124671151_prec,0.)
call test_one_MPL((/ 2,1,2 /),(/ (0.03_prec,0.), (0.5012562893380046_prec,0.), (55.3832_prec,0.) /),ref, '1.3')
ref = cmplx(-0.06565799418838372)
call test_one_MPL((/1, 1/), cmplx((/-0.25,-2. /)), ref, '1.4')
ref = cmplx(-0.03199896396564833)
call test_one_MPL((/2, 1/), cmplx((/-0.25,-2. /)), ref, '1.4')
ref = (-0.06565799418838388900030511939327151809905663800733251335678_prec,0.)
call test_one_MPL((/1, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.4')
ref = (-0.0319989639656484125354488353677231989818697987217232462808_prec,0.)
call test_one_MPL((/2, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.4')
end subroutine do_MPL_tests
subroutine test_one_condensed(m,z,y,k,ref,test_id)
......@@ -108,98 +107,100 @@ CONTAINS
real(kind=prec) :: z, xchen
print*, 'doing GPL tests...'
ref = cmplx(0.0819393734128676)
call test_one_condensed((/ 1,1 /),cmplx((/ 1.3, 1.1 /)),cmplx(0.4),2,ref,'2.1')
ref = (0.0819393734128676670196050250748145117414_prec,0.)
call test_one_condensed((/ 1,1 /),(/ (1.3_prec,0.), (1.1_prec,0.) /),(0.4_prec,0.),2,ref,'2.1')
ref = cmplx(0.01592795952537145)
call test_one_condensed((/ 3,2 /),cmplx((/ 1.3, 1.1 /)),cmplx(0.4),2,ref,'2.2')
ref = cmplx(0.0020332632172573974)
call test_one_condensed((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3')
ref = (0.0159279595253714815409530431586662402173_prec,0.)
call test_one_condensed((/ 3,2 /),(/ (1.3_prec,0.), (1.1_prec,0.) /),(0.4_prec,0.),2,ref,'2.2')
ref = (0.0020332632172573960663410568998748341439_prec,0.)
call test_one_condensed((/ 4 /),(/ (0._prec,0.) /),(1.6_prec,0.),1,ref,'2.3')
! requires making convergent
ref = (0.09593041677639341, -0.8829351795197851)
call test_one_flat(cmplx([0,1,3,2]),ref,'2.4')
ref = (0.0959304167763934268889451723647597269243_prec,-0.882935179519785044294098901083860819793_prec)
call test_one_flat(cmplx([0,1,3,2],kind=prec),ref,'2.4')
ref = (0.009947947789928621,0.0)
call test_one_flat(cmplx([0.0, 0.0, -3.3333333333333335, -3.3333333333333335, 1.0]),ref,'2.5')
ref = (0.00994794778992862285848386211876302459484_prec,0.0)
call test_one_flat([ (0._prec,0.), (0._prec,0.), cmplx(-10._prec/3._prec,kind=prec), cmplx(-10._prec/3._prec,kind=prec), (1._prec,0.) ],ref,'2.5')
! requires hoelder convolution
ref = (-0.012709942828250949,0.0)
call test_one_flat(cmplx([0.0, 3.3333333333333335, 1.0, 3.3333333333333335, 1.0]),ref,'2.6')
ref = (-0.01270994282825128082040127705282437854604_prec,0.0)
call test_one_flat([ (0._prec,0.), cmplx(10._prec/3._prec,kind=prec), (1._prec,0.), cmplx(10._prec/3._prec,kind=prec), (1._prec,0.) ],ref,'2.6')
! here the tests from the mathematica nb start
! --------------------------------------------------------------------------
z = 1./200; xchen = 0.3;
z = 1._prec/200; xchen = 0.3_prec;
ref = (-0.0050125418235441935,0.0)
call test_one_flat(cmplx([1./z,1.0]),ref,'3.1')
ref = (-0.00501254182354428204309373895836778138658_prec,0.0)
call test_one_flat(cmplx([1./z,1.0_prec],kind=prec),ref,'3.1')
ref = (-0.0015011261262671913,0.0)
call test_one_flat(cmplx([1./(xchen*z),1.0]),ref,'3.2')
ref = (-0.001501126126267145650881556104825717433567_prec,0.0)
call test_one_flat(cmplx([1./(xchen*z),1.0_prec],kind=prec),ref,'3.2')
ref = (-0.0007502860817810596,0.0)
call test_one_flat(cmplx([(1+sqrt(1-z**2))/(xchen*z),1.0]),ref,'3.3')
ref = (-7.50286081781062398521044689345496255877E-4_prec,0.0)
call test_one_flat(cmplx([(1+sqrt(1-z**2))/(xchen*z),1.0_prec],kind=prec),ref,'3.3')
ref = (0.0074335969894765335,0.0)
call test_one_flat(cmplx([-1./xchen,-1./xchen,1.,1.,1.0]),ref,'3.4')
ref = (0.0074335969894769006767497160582533005912_prec,0.0)
call test_one_flat(cmplx([-1._prec/xchen,-1._prec/xchen,1._prec,1._prec,1._prec],kind=prec),ref,'3.4')
ref = (-8.403785974849544e-6,0.0)
call test_one_flat(cmplx([-1./xchen,0.,-1./xchen,1./(xchen*z),1.0]),ref,'3.5')
ref = (-8.40378597494480117619893114520705717404E-6_prec,0.0)
call test_one_flat(cmplx([-1._prec/xchen,0._prec,-1._prec/xchen,1._prec/(xchen*z),1.0_prec],kind=prec),ref,'3.5')
ref = cmplx((0.4925755847450199,2.6389214054743295))
call test_one_flat(cmplx([-1.,-1.,z,z,1.]),ref,'3.6')
ref = (0.492575584745015188956099120118429131062_prec,2.63892140547433046711653044177665570864_prec)
call test_one_flat(cmplx([-1._prec,-1._prec,z,z,1._prec],kind=prec),ref,'3.6')
ref = cmplx((-0.0015317713178859967,-0.00045003911367000565))
call test_one_flat(cmplx([0.,0.,(1-sqrt(1-z**2))/(xchen*z), 1./(xchen*z),1.]),ref,'3.7')
ref = (-0.001531771317885848316868878120481266022195_prec,-4.50039113668843745227527662921010599E-4_prec)
call test_one_flat(cmplx([0._prec,0._prec,(1._prec-sqrt(1-z**2))/(xchen*z), 1._prec/(xchen*z),1._prec],kind=prec),ref,'3.7')
! here the chen integral tests start
! ----------------------------------------------------------------------
ref = (-1.2039728043259361,0.0)
call test_one_flat(cmplx([0., 0.3]),ref,'4.1')
z = 1/100._prec
ref = (-1.203972804325935992622746217761838502954_prec,0.0)
call test_one_flat((/ (0._prec,0.), (0.3_prec,0.) /),ref,'4.1')
ref = (10.603796220956799,0.0)
call test_one_flat(cmplx([0., 0., 0.01]),ref,'4.2')
ref = (10.6037962209567960211233327771880353832_prec,0.0)
call test_one_flat( (/ (0._prec,0.), (0._prec,0.), (0.01_prec,0.) /),ref,'4.2')
ref = (0.0005042990065180765,0.0)
call test_one_flat(cmplx([100., 1., 0.3]),ref,'4.3')
ref = (5.0429900651807654463744449762857296353E-4_prec,0.0)
call test_one_flat((/(100._prec,0.), (1._prec,0.), (0.3_prec,0.) /),ref,'4.3')
ref = (0.05630861877185141,0.0)
call test_one_flat(cmplx([1., 0., 0.01]),ref,'4.4')
ref = (0.05630861877185110944118569266547272128_prec,0.0)
call test_one_flat( (/ (1._prec,0.), (0._prec,0.), (0.01_prec,0.) /) ,ref,'4.4')
ref = (0.04032895150872735,0.2922709647568842)
call test_one_flat([(0.01,0.9999499987499375), (0.3,0)],ref,'4.5')
ref = (0.0403289515087272334386799760345612761274_prec,0.29227096475688417627233177065883250904_prec)
call test_one_flat([cmplx(z,sqrt(1._prec-z**2),kind=prec), (0.3_prec,0._prec)],ref,'4.5')
ref = cmplx(0.000025210645098340354)
call test_one_flat(cmplx([100., 199.99499987499374, 1.]),ref,'4.6')
ref = (2.52106450983403521007132254962389653479E-5_prec,0.)
call test_one_flat(cmplx([1._prec/z, (1+sqrt(1-z**2))/z, 1._prec], kind=prec),ref,'4.6')
ref = (0.0556470547632135,0.)
call test_one_flat(cmplx([-1., 0.01, 0., 0.01]),ref,'4.7')
ref = (0.0556470547634805907802456701296765199105_prec,0.)
call test_one_flat((/ (-1._prec,0.), (0.01_prec,0.), (0._prec,0.), (0.01_prec,0.) /),ref,'4.7')
ref = (0.00003794895846844715,0.)
call test_one_flat(cmplx([100., 100., 1., 0., 1.]),ref,'4.8')
ref = (3.79489584700899518570534417045601144619E-5_prec,0.)
call test_one_flat((/ (100._prec,0.), (100._prec,0.), (1._prec,0.), (0._prec,0.), (1._prec,0.) /),ref,'4.8')
ref = (0.00007574284433216497,0.)
call test_one_flat(cmplx([100., 1., 100., 0., 1.]),ref,'4.9')
ref = (7.57428443370789293292663421167820669347E-5_prec,0.)
call test_one_flat((/ (100._prec,0.), (1._prec,0.), (100._prec,0.), (0._prec,0.), (1._prec,0.) /),ref,'4.9')
ref = (1.8801911586012443,2.509434138598062)
call test_one_flat(cmplx([0.01, -1.0, 0.01, 1.]),ref,'4.10')
ref = (1.88019115860102242751826218248903485831_prec,2.50943413859806280406885100241705929527_prec)
call test_one_flat((/ (0.01_prec,0.), (-1.0_prec,0.), (0.01_prec,0.), (1._prec,0.) /),ref,'4.10')
ref = (-0.012539108315054982, -0.015414250168437678)
call test_one_flat(cmplx([0.01, 199.99499987499374, 0.01, 1.]),ref,'4.11')
ref = (-0.0125391083150538283259376060156453327489_prec,-0.0154142501684376766660782277248386124052_prec)
call test_one_flat(cmplx([z, (1+sqrt(1-z**2))/z, z, 1._prec],kind=prec),ref,'4.11')
! Increase test coverage
ref = (-1.2039728043259361, pi)
call test_one_flat(cmplx([0., -0.3]),ref,'5.1')
ref = (-0.5394048306859651, 1.0303768265243125)
call test_one_flat([(0.,0.), (0.3, 0.5)],ref,'5.2')
ref = (-1.203972804325935992622746217761838502953611_prec, pi)
call test_one_flat((/ (0._prec,0.), (-0.3_prec,0.) /),ref,'5.1')
ref = (0.00003794895846844715,0.)
ref = (-0.539404830685964989184608085018712655769250_prec, 1.0303768265243124637877433270311515319634364_prec)
call test_one_flat([(0._prec,0._prec), (0.3_prec, 0.5_prec)],ref,'5.2')
ref = (3.79489584700899518570534417045601144619E-5_prec,0.)
print*, ' ', 'testing GPL ', '5.3', ' ...'
res = G((/100., 100., 1., 0., 1./))
res = G((/100._prec, 100._prec, 1._prec, 0._prec, 1._prec/))
call check(res,ref)
print*, ' ', 'testing GPL ', '5.4', ' ...'
......@@ -211,29 +212,34 @@ CONTAINS
call check(res,ref)
print*, ' ', 'testing GPL ', '5.6', ' ...'
res = G((/100.,100.,1.,0./), 1.)
res = G((/100._prec,100._prec,1._prec,0._prec/), 1._prec)
call check(res,ref)
print*, ' ', 'testing GPL ', '5.7', ' ...'
res = G(cmplx((/100.,100.,1.,0./)), cmplx(1.))
res = G((/(100._prec,0.),(100._prec,0.),(1._prec,0.),(0._prec,0.)/), (1._prec,0.))
call check(res,ref)
print*, ' ', 'testing GPL ', '5.8', ' ...'
res = G( (/1,1,1,1/) , (/ 100., 100., 1., 0. /), 1.)
res = G( (/1,1,1,1/) , (/ 100._prec, 100._prec, 1._prec, 0._prec /), 1._prec)
call check(res,ref)
!ref = cmplx(0.01592795952537145)
ref = cmplx(0.048503553116818435, - 2.8232688511893116)
ref = (0.0485035531168477720320998551314479928648_prec,-2.82326885118931135859594797186107474877_prec)
print*, ' ', 'testing GPL ', '5.9', ' ...'
res = G( (/ 3, 2 /), toinum((/ 1.3, 1.1 /)), toinum(4.) )
res = G( (/ 3, 2 /), toinum((/ 1.3_prec, 1.1_prec /)), toinum(4._prec) )
call check(res,ref)
print*, ' ', 'testing GPL ', '5.10', ' ...'
res = G( (/ 3, 2 /), (/ 1.3, 1.1 /), 4. )
res = G( (/ 3, 2 /), (/ 1.3_prec, 1.1_prec /), 4._prec )
call check(res,ref)
print*, ' ', 'testing GPL ', '5.11', ' ...'
res = G( (/ 3, 2 /), cmplx((/ 1.3, 1.1 /)), cmplx(4.) )
res = G( (/ 3, 2 /), (/ (1.3_prec,0.), (1.1_prec,0.) /), (4._prec,0.) )
call check(res,ref)
print*, ' ', 'testing GPL ', '5.12', ' ...'
ref = (0.107961231635965271810236308594279564021882373573100109007684973130605145011141_prec, 0.)
res = G((/ 0._prec, 1.3_prec, 0._prec, 0._prec, 4._prec, 1.1_prec /))
call check(res,ref)
......@@ -244,6 +250,7 @@ CONTAINS
#ifdef HAVE_GINAC
#ifdef HAVE_MM
function evalt(arr, what)
implicit none
......@@ -292,9 +299,9 @@ CONTAINS
use gtestmuonenp, only: inimuonenp=>args
implicit none
tol = 6.0e-6
call perform_ginacv( 6, inichen (cmplx(0.3), cmplx(0.1)) )
call perform_ginacv( 7, inimuone (cmplx(0.5), cmplx(0.6)) )
call perform_ginacv( 8, inimuonenp(cmplx(0.3), cmplx(0.6)) )
call perform_ginacv( 6, inichen ((0.3_prec,0.), (0.1_prec,0.)) )
call perform_ginacv( 7, inimuone ((0.5_prec,0.), (0.6_prec,0.)) )
call perform_ginacv( 8, inimuonenp((0.3_prec,0.), (0.6_prec,0.)) )
end subroutine
......@@ -314,7 +321,7 @@ CONTAINS
res=evalt(args(j,:,i),0)
enddo
call system_clock(cend, count_rate=count_rate)
time(1) = real(cend-cstart)/count_rate/size(args,3)
time(1) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
if (time(1).lt.zero) print*,j, cend-cstart,count_rate
ttime(1) = ttime(1) + time(1)
......@@ -323,7 +330,7 @@ CONTAINS
res=evalt(args(j,:,i),1)
enddo
call system_clock(cend, count_rate=count_rate)
time(2) = real(cend-cstart)/count_rate/size(args,3)
time(2) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
if (time(2).lt.zero) print*,j
ttime(2) = ttime(2) + time(2)
......@@ -332,6 +339,7 @@ CONTAINS
write( * , 900, advance='no' ) cr, j, size(args,1), msg
enddo
print*,
print*,ttime
write(*,901) msg, size(args,1)/ttime(2)/1000., size(args,1)/ttime(1)/1000., int(ttime(2)/ttime(1))
! Lets to another, fair comparison
......@@ -340,19 +348,19 @@ CONTAINS
res=evalt(args(j,:,1),0)
enddo
call system_clock(cend, count_rate=count_rate)
ttime(1) = real(cend-cstart)/count_rate
ttime(1) = real(cend-cstart)/real(count_rate,kind=prec)
call system_clock(cstart, count_rate=count_rate)
do j=1,size(args,1)
res=evalt(args(j,:,1),1)
enddo
call system_clock(cend, count_rate=count_rate)
ttime(2) = real(cend-cstart)/count_rate
ttime(2) = real(cend-cstart)/real(count_rate,kind=prec)
write(*,901) msg, size(args,1)/ttime(2)/1000., size(args,1)/ttime(1)/1000., int(ttime(2)/ttime(1))
900 FORMAT(a , 'Function ',i4,'/',i4,' for ',a)
901 format('Evaluating ',A,' using GiNaC at ',F9.2,'kG/s and GPL at ',F9.2,'kG/s (',I3,'x)')
901 format('Evaluating ',A,' using GiNaC at ',F9.2,'kG/s and handyG at ',F9.2,'kG/s (',I3,'x)')
end subroutine
subroutine do_timing_tests(n)
......@@ -374,16 +382,16 @@ CONTAINS
do i=1,n
z = ran2(ranseed) / 2.
x = ran2(ranseed)*(1-z) + z
cargs(:,:,i) = inichen (cmplx(x), cmplx(z))
fargs(:,:,i) = inichenff(cmplx(x), cmplx(z))
cargs(:,:,i) = inichen (cmplx(x,kind=prec), cmplx(z,kind=prec))
fargs(:,:,i) = inichenff(cmplx(x,kind=prec), cmplx(z,kind=prec))
w = ran2(ranseed) ! 0<w<1
z = ran2(ranseed) * (sqrt(1-w+w**2)-sqrt(w)) + sqrt(w)
nargs(:,:,i) = inimuonenp(cmplx(w), cmplx(z))
nargs(:,:,i) = inimuonenp(cmplx(w,kind=prec), cmplx(z,kind=prec))
x = ran2(ranseed)
y = ran2(ranseed)
pargs(:,:,i) = inimuone(cmplx(x), cmplx(y))
pargs(:,:,i) = inimuone(cmplx(x,kind=prec), cmplx(y,kind=prec))
enddo
cargs(1181,:,:)=1.e15
......@@ -422,6 +430,7 @@ CONTAINS
#endif
#endif
! subroutine do_shuffle_tests()
! complex(kind=prec) :: v(2) = cmplx((/1,2/))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment