test.f90 12.5 KB
Newer Older
Luca's avatar
Luca committed
1

Luca's avatar
Luca committed
2
! These tests assume that GPLInfinity = 30
Luca's avatar
Luca committed
3

Luca's avatar
Luca committed
4
PROGRAM TEST
5
  use globals
6
  use utils
7
  use shuffle
8
  use maths_functions
Luca's avatar
Luca committed
9 10
  use mpl_module
  use gpl_module
ulrich_y's avatar
ulrich_y committed
11
  use chenreftest, only: do_chen_test
Luca's avatar
Luca committed
12
  implicit none
Luca's avatar
Luca committed
13
  
Luca's avatar
Luca committed
14
  complex(kind=prec) :: res 
ulrich_y's avatar
ulrich_y committed
15
  real :: tol = 8.0e-7
Luca's avatar
Luca committed
16
  logical :: tests_successful = .true. 
17

ulrich_y's avatar
ulrich_y committed
18
#ifdef DEBUG
Luca's avatar
tests  
Luca committed
19
  call parse_cmd_args()
ulrich_y's avatar
ulrich_y committed
20
#endif
ulrich_y's avatar
ulrich_y committed
21 22

  tol = 8e-10
Luca's avatar
tests  
Luca committed
23 24
  call do_MPL_tests() 
  call do_GPL_tests()
25
  ! call do_shuffle_tests() ! put this somewhere else
ulrich_y's avatar
ulrich_y committed
26
  tol = 8.0e-7
ulrich_y's avatar
ulrich_y committed
27
  tests_successful = tests_successful .and. do_chen_test(cmplx(0.3),cmplx(0.1))
Luca's avatar
tests  
Luca committed
28

ulrich_y's avatar
ulrich_y committed
29
#ifdef HAVE_GINAC
ulrich_y's avatar
ulrich_y committed
30 31
#ifdef HAVE_MM
  call do_ginac_tests
ulrich_y's avatar
ulrich_y committed
32
  call do_timing_tests(5)
ulrich_y's avatar
ulrich_y committed
33
#endif
ulrich_y's avatar
ulrich_y committed
34
#endif
Luca's avatar
tests  
Luca committed
35 36 37 38 39 40 41

  if(tests_successful) then
    print*, 'All tests passed. '
  else 
    print*, 'Some tests failed. '
    stop 1
  end if
42

43
CONTAINS
44
   
Luca Naterop's avatar
Luca Naterop committed
45 46 47 48
  subroutine check(res, ref)
    complex(kind=prec) :: res, ref
    real(kind=prec) :: delta

Luca Naterop's avatar
Luca Naterop committed
49
    delta = abs(res-ref)
Luca Naterop's avatar
Luca Naterop committed
50
    if(delta < tol) then
Luca Naterop's avatar
Luca Naterop committed
51
      print*, '       ',' passed with delta = ', delta
Luca Naterop's avatar
Luca Naterop committed
52
    else 
53
      print*, '  ',' FAILED with delta = ', delta
54
      tests_successful = .false.
Luca Naterop's avatar
Luca Naterop committed
55 56
    end if
  end subroutine check
57 58 59 60 61

  subroutine test_one_MPL(m,x,ref, test_id)
    integer :: m(:)
    complex(kind=prec) :: x(:), ref
    character(len=*) :: test_id
Luca's avatar
Luca committed
62
    
63 64 65 66 67
    print*, '  ', 'testing MPL ', test_id, ' ...'
    res = MPL(m,x)
    call check(res,ref)
  end subroutine test_one_MPL

Luca's avatar
Luca committed
68
  subroutine do_MPL_tests()
Luca's avatar
Luca committed
69
    complex(kind=prec) :: ref
Luca Naterop's avatar
Luca Naterop committed
70
    print*, 'doing MPL tests...'
Luca's avatar
Luca committed
71
    
ulrich_y's avatar
ulrich_y committed
72
    ref = cmplx(0.022696600480693277651633)
73
    call test_one_MPL((/ 1,1 /),cmplx((/ 0.3156498673740053, 0.3431255827785649 /)),ref, '1.1')
Luca's avatar
Luca committed
74
    
ulrich_y's avatar
ulrich_y committed
75
    ref = cmplx(0.00023134615630308335448329926098409)
76
    call test_one_MPL((/ 1,1 /),cmplx((/ 0.03, 0.5012562893380046 /)),ref, '1.2')
Luca's avatar
Luca committed
77
    
ulrich_y's avatar
ulrich_y committed
78
    ref = cmplx(0.000023446106415452030937059124671151)
Luca's avatar
Luca committed
79
    call test_one_MPL((/ 2,1,2 /),cmplx((/ 0.03, 0.5012562893380046, 55.3832 /)),ref, '1.3')  
ulrich_y's avatar
ulrich_y committed
80 81 82 83 84

    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')
Luca's avatar
Luca committed
85
  end subroutine do_MPL_tests
Luca's avatar
Luca committed
86

Luca's avatar
Luca committed
87
  subroutine test_one_condensed(m,z,y,k,ref,test_id)
88 89 90
    integer :: m(:), k
    complex(kind=prec) :: z(:), y, res, ref
    character(len=*) :: test_id
Luca Naterop's avatar
Luca Naterop committed
91

92
    print*, '  ', 'testing GPL ', test_id, ' ...'
ulrich_y's avatar
ulrich_y committed
93
    res = G_condensed(m,toinum(z),inum(y,di0),k)
Luca Naterop's avatar
Luca Naterop committed
94
    call check(res,ref)
Luca's avatar
Luca committed
95 96
  end subroutine test_one_condensed

97 98
  subroutine test_one_flat(z,ref,test_id)
    complex(kind=prec) :: z(:), res, ref
Luca's avatar
Luca committed
99 100 101
    character(len=*) :: test_id

    print*, '  ', 'testing GPL ', test_id, ' ...'
ulrich_y's avatar
ulrich_y committed
102
    res = G(z)
Luca's avatar
Luca committed
103 104
    call check(res,ref)
  end subroutine test_one_flat
Luca Naterop's avatar
Luca Naterop committed
105 106

  subroutine do_GPL_tests()
ulrich_y's avatar
ulrich_y committed
107
    complex(kind=prec) :: ref, res
Luca Naterop's avatar
Luca Naterop committed
108
    real(kind=prec) :: z, xchen
Luca Naterop's avatar
Luca Naterop committed
109
    print*, 'doing GPL tests...'
Luca's avatar
Luca committed
110
    
Luca Naterop's avatar
Luca Naterop committed
111
    ref = cmplx(0.0819393734128676)
112
    call test_one_condensed((/ 1,1 /),cmplx((/ 1.3, 1.1 /)),cmplx(0.4),2,ref,'2.1')
Luca's avatar
Luca committed
113
    
Luca Naterop's avatar
Luca Naterop committed
114
    ref = cmplx(0.01592795952537145)
115
    call test_one_condensed((/ 3,2 /),cmplx((/ 1.3, 1.1 /)),cmplx(0.4),2,ref,'2.2')
116
    
Luca Naterop's avatar
Luca Naterop committed
117
    ref = cmplx(0.0020332632172573974)
Luca's avatar
tests  
Luca committed
118
    call test_one_condensed((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3')
Luca's avatar
Luca committed
119

Luca Naterop's avatar
Luca Naterop committed
120
    ! requires making convergent
121
    ref = (0.09593041677639341, -0.8829351795197851)
122
    call test_one_flat(cmplx([0,1,3,2]),ref,'2.4')
Luca Naterop's avatar
Luca Naterop committed
123

124
    ref = (0.009947947789928621,0.0)
Luca Naterop's avatar
Luca Naterop committed
125 126 127
    call test_one_flat(cmplx([0.0, 0.0, -3.3333333333333335, -3.3333333333333335, 1.0]),ref,'2.5')

    ! requires hoelder convolution
128
    ref = (-0.012709942828250949,0.0)
Luca Naterop's avatar
Luca Naterop committed
129 130 131
    call test_one_flat(cmplx([0.0, 3.3333333333333335, 1.0, 3.3333333333333335, 1.0]),ref,'2.6')

    ! here the tests from the mathematica nb start
132 133
    ! --------------------------------------------------------------------------

Luca Naterop's avatar
Luca Naterop committed
134 135
    z = 1./200; xchen = 0.3;

136
    ref = (-0.0050125418235441935,0.0)
Luca Naterop's avatar
Luca Naterop committed
137 138
    call test_one_flat(cmplx([1./z,1.0]),ref,'3.1')

139
    ref = (-0.0015011261262671913,0.0)
Luca Naterop's avatar
Luca Naterop committed
140 141
    call test_one_flat(cmplx([1./(xchen*z),1.0]),ref,'3.2')

142
    ref = (-0.0007502860817810596,0.0)
Luca Naterop's avatar
Luca Naterop committed
143 144
    call test_one_flat(cmplx([(1+sqrt(1-z**2))/(xchen*z),1.0]),ref,'3.3')

145
    ref = (0.0074335969894765335,0.0)
Luca Naterop's avatar
Luca Naterop committed
146 147
    call test_one_flat(cmplx([-1./xchen,-1./xchen,1.,1.,1.0]),ref,'3.4')
    
148
    ref = (-8.403785974849544e-6,0.0)
Luca Naterop's avatar
Luca Naterop committed
149 150
    call test_one_flat(cmplx([-1./xchen,0.,-1./xchen,1./(xchen*z),1.0]),ref,'3.5')
    
151 152
    ref = cmplx((0.4925755847450199,2.6389214054743295))
    call test_one_flat(cmplx([-1.,-1.,z,z,1.]),ref,'3.6')
Luca Naterop's avatar
Luca Naterop committed
153
    
Luca Naterop's avatar
Luca Naterop committed
154 155
    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')
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191


    ! here the chen integral tests start
    ! ----------------------------------------------------------------------

    ref = (-1.2039728043259361,0.0)
    call test_one_flat(cmplx([0., 0.3]),ref,'4.1')

    ref = (10.603796220956799,0.0)
    call test_one_flat(cmplx([0., 0., 0.01]),ref,'4.2')

    ref = (0.0005042990065180765,0.0)
    call test_one_flat(cmplx([100., 1., 0.3]),ref,'4.3')

    ref = (0.05630861877185141,0.0)
    call test_one_flat(cmplx([1., 0., 0.01]),ref,'4.4')

    ref = (0.04032895150872735,0.2922709647568842)
    call test_one_flat([(0.01,0.9999499987499375), (0.3,0)],ref,'4.5')

    ref = cmplx(0.000025210645098340354)
    call test_one_flat(cmplx([100., 199.99499987499374, 1.]),ref,'4.6')

    ref = (0.0556470547632135,0.)
    call test_one_flat(cmplx([-1., 0.01, 0., 0.01]),ref,'4.7')

    ref = (0.00003794895846844715,0.)
    call test_one_flat(cmplx([100., 100., 1., 0., 1.]),ref,'4.8')

    ref = (0.00007574284433216497,0.)
    call test_one_flat(cmplx([100., 1., 100., 0., 1.]),ref,'4.9')

    ref = (1.8801911586012443,2.509434138598062)
    call test_one_flat(cmplx([0.01, -1.0, 0.01, 1.]),ref,'4.10')

    ref = (-0.012539108315054982, -0.015414250168437678)
Luca Naterop's avatar
Luca Naterop committed
192
    call test_one_flat(cmplx([0.01, 199.99499987499374, 0.01, 1.]),ref,'4.11')
193

ulrich_y's avatar
ulrich_y committed
194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238
    ! 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 = (0.00003794895846844715,0.)
    print*, '  ', 'testing GPL ', '5.3', ' ...'
    res = G((/100., 100., 1., 0., 1./))
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.4', ' ...'
    res = G((/100,100,1,0,1/))
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.5', ' ...'
    res = G_superflatn(toinum((/100,100,1,0,1/)), 5)
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.6', ' ...'
    res = G((/100.,100.,1.,0./), 1.)
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.7', ' ...'
    res = G(cmplx((/100.,100.,1.,0./)), cmplx(1.))
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.8', ' ...'
    res = G( (/1,1,1,1/) , (/ 100., 100., 1., 0. /), 1.)
    call check(res,ref)

    !ref = cmplx(0.01592795952537145)
    ref = cmplx(0.048503553116818435, - 2.8232688511893116)
    print*, '  ', 'testing GPL ', '5.9', ' ...'
    res = G( (/ 3, 2 /), toinum((/ 1.3, 1.1 /)), toinum(4.) )
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.10', ' ...'
    res = G( (/ 3, 2 /), (/ 1.3, 1.1 /), 4. )
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.11', ' ...'
    res = G( (/ 3, 2 /), cmplx((/ 1.3, 1.1 /)), cmplx(4.) )
    call check(res,ref)

Luca Naterop's avatar
Luca Naterop committed
239

Luca Naterop's avatar
Luca Naterop committed
240 241 242

  end subroutine do_GPL_tests
  
ulrich_y's avatar
ulrich_y committed
243 244 245



ulrich_y's avatar
ulrich_y committed
246
#ifdef HAVE_GINAC
ulrich_y's avatar
ulrich_y committed
247 248 249 250 251 252 253 254 255

  function evalt(arr, what)
    implicit none
    complex(kind=prec) :: arr(:), evalt, geval
    integer what, i, l

    evalt =0.

    do i=1,size(arr)
256
      if (abs(arr(i)) .gt. 1.e10) then ! isnan?
ulrich_y's avatar
ulrich_y committed
257 258 259 260 261 262 263 264 265 266 267 268 269
        l = i-1
        goto 123
      endif
    enddo
    l = size(arr)

123 continue
    if (l==0) return

    if (what .eq. 0) then
      evalt = G(arr(1:l))
    elseif (what.eq.1) then
      evalt = geval(arr(1:l),l)
ulrich_y's avatar
ulrich_y committed
270
    endif
ulrich_y's avatar
ulrich_y committed
271 272 273
  end function

  subroutine perform_ginacv(n, args)
ulrich_y's avatar
ulrich_y committed
274
    use maths_functions, only:clearcache
ulrich_y's avatar
ulrich_y committed
275 276
    complex(kind=prec) :: args(:,:)
    integer i,n
ulrich_y's avatar
ulrich_y committed
277
    call clearcache
ulrich_y's avatar
ulrich_y committed
278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298

    do i=1,size(args,1)
      write(*,900) n,i
      call check(evalt(args(i,:),0),evalt(args(i,:),1))
    enddo

900 format("   testing GPL ",I1,".",I4.4," ...")

  end subroutine

  subroutine do_ginac_tests
    use maths_functions, only:clearcache
    use gtestchen   , only: inichen   =>args
    use gtestmuone  , only: inimuone  =>args
    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)) )

ulrich_y's avatar
ulrich_y committed
299
  end subroutine
ulrich_y's avatar
ulrich_y committed
300

ulrich_y's avatar
ulrich_y committed
301 302 303 304
  subroutine do_one_speed_test(args, u, msg)
    use maths_functions, only:clearcache
    implicit none
    complex(kind=prec) :: args(:,:,:),res
ulrich_y's avatar
ulrich_y committed
305 306
    integer(kind=8) cstart, cend, count_rate
    real(kind=prec) :: time(2), ttime(2)
ulrich_y's avatar
ulrich_y committed
307 308 309 310 311
    integer i,j, u
    character,parameter :: cr = achar(13)
    character(len=*) msg
    do j=1,size(args,1)
      ! try function a bunch of times
ulrich_y's avatar
ulrich_y committed
312
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
313 314 315
      do i=1,size(args,3)
        res=evalt(args(j,:,i),0)
      enddo
ulrich_y's avatar
ulrich_y committed
316 317 318
      call system_clock(cend, count_rate=count_rate)
      time(1) = real(cend-cstart)/count_rate/size(args,3)
      if (time(1).lt.zero) print*,j, cend-cstart,count_rate
ulrich_y's avatar
ulrich_y committed
319 320
      ttime(1) = ttime(1) + time(1)

ulrich_y's avatar
ulrich_y committed
321
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
322 323 324
      do i=1,size(args,3)
        res=evalt(args(j,:,i),1)
      enddo
ulrich_y's avatar
ulrich_y committed
325 326
      call system_clock(cend, count_rate=count_rate)
      time(2) = real(cend-cstart)/count_rate/size(args,3)
ulrich_y's avatar
ulrich_y committed
327 328 329 330 331 332 333 334 335 336
      if (time(2).lt.zero) print*,j
      ttime(2) = ttime(2) + time(2)

      write(u,*) time

      write( * , 900, advance='no' ) cr, j, size(args,1), msg
    enddo
    print*,
    write(*,901) msg, size(args,1)/ttime(2)/1000., size(args,1)/ttime(1)/1000., int(ttime(2)/ttime(1))

ulrich_y's avatar
ulrich_y committed
337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353
    ! Lets to another, fair comparison
    call system_clock(cstart, count_rate=count_rate)
    do j=1,size(args,1)
      res=evalt(args(j,:,1),0)
    enddo
    call system_clock(cend, count_rate=count_rate)
    ttime(1) = real(cend-cstart)/count_rate

    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

    write(*,901) msg, size(args,1)/ttime(2)/1000., size(args,1)/ttime(1)/1000., int(ttime(2)/ttime(1))

ulrich_y's avatar
ulrich_y committed
354 355 356 357 358 359
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)')
  end subroutine

  subroutine do_timing_tests(n)
    use gtestchen   , only: inichen   =>args
ulrich_y's avatar
ulrich_y committed
360
    use gtestchenff , only: inichenff =>args
ulrich_y's avatar
ulrich_y committed
361 362 363 364 365 366
    use gtestmuone  , only: inimuone  =>args
    use gtestmuonenp, only: inimuonenp=>args
    implicit none
    integer, intent(in) :: n
    integer i
    complex(kind=prec) :: cargs( 1399,5,n)
ulrich_y's avatar
ulrich_y committed
367
    complex(kind=prec) :: fargs(  540,5,n)
ulrich_y's avatar
ulrich_y committed
368 369 370 371 372 373 374 375 376
    complex(kind=prec) :: pargs(  198,5,n)
    complex(kind=prec) :: nargs( 1733,5,n)
    real(kind=prec) :: z, x, y, w
    integer ranseed
    ranseed = 233123

    do i=1,n
      z = ran2(ranseed) / 2.
      x = ran2(ranseed)*(1-z) + z
ulrich_y's avatar
ulrich_y committed
377 378
      cargs(:,:,i) = inichen  (cmplx(x), cmplx(z))
      fargs(:,:,i) = inichenff(cmplx(x), cmplx(z))
ulrich_y's avatar
ulrich_y committed
379 380 381 382 383 384 385 386 387 388

      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))

      x = ran2(ranseed)
      y = ran2(ranseed)
      pargs(:,:,i) = inimuone(cmplx(x), cmplx(y))
    enddo

ulrich_y's avatar
ulrich_y committed
389 390
    cargs(1181,:,:)=1.e15
    fargs(367,:,:)=1.e15
ulrich_y's avatar
ulrich_y committed
391 392 393


    open(unit=9, file="stats.txt")
ulrich_y's avatar
ulrich_y committed
394 395
    write(9,*) "Chen form factor"
    call do_one_speed_test(fargs,9,"Chen FF")
ulrich_y's avatar
ulrich_y committed
396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423
    write(9,*) "Chen"
    call do_one_speed_test(cargs,9,"Chen")
    write(9,*) "MUonE-planar"
    call do_one_speed_test(pargs,9,"Muone planar")
    write(9,*) "MUonE-non-planar"
    call do_one_speed_test(nargs,9,"Muone non planar")
    close(unit=9)

  end subroutine

      FUNCTION RAN2(randy)

            ! This is the usual "random"

      implicit none
      real(kind=prec) :: MINV,RAN2
      integer m,a,Qran,r,hi,lo,randy
      PARAMETER(M=2147483647,A=16807,Qran=127773,R=2836)
      PARAMETER(MINV=0.46566128752458e-09)
      HI = RANDY/Qran
      LO = MOD(RANDY,Qran)
      RANDY = A*LO - R*HI
      IF(RANDY.LE.0) RANDY = RANDY + M
      RAN2 = RANDY*MINV
      END FUNCTION RAN2



ulrich_y's avatar
ulrich_y committed
424

ulrich_y's avatar
ulrich_y committed
425
#endif
426 427 428
  ! subroutine do_shuffle_tests() 
  !   complex(kind=prec) :: v(2) = cmplx((/1,2/))
  !   complex(kind=prec) :: w(2) = cmplx((/3,4/))
429

430 431
  !   call print_matrix(shuffle_product(v,w))
  ! end subroutine do_shuffle_tests
432

Luca's avatar
Luca committed
433
END PROGRAM TEST
Luca's avatar
Luca committed
434 435

! In terminal kann man den exit code bekommen via echo $?