test.f90 14.8 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
  
14
  real(kind=prec) :: tol = 8.0e-7
Luca's avatar
Luca committed
15
  logical :: tests_successful = .true. 
16

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

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

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

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

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

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

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

Luca's avatar
Luca committed
67
  subroutine do_MPL_tests()
Luca's avatar
Luca committed
68
    complex(kind=prec) :: ref
Luca Naterop's avatar
Luca Naterop committed
69
    print*, 'doing MPL tests...'
Luca's avatar
Luca committed
70
    
71 72
    ref = (0.022696600480693277651633_prec,0.)
    call test_one_MPL((/ 1,1 /),(/ (0.3156498673740053_prec,0.), (0.3431255827785649_prec,0.) /),ref, '1.1')
Luca's avatar
Luca committed
73
    
74 75
    ref = (0.00023134615630308335448329926098409_prec,0.)
    call test_one_MPL((/ 1,1 /),(/ (0.03_prec,0.), (0.5012562893380046_prec,0.) /),ref, '1.2')
Luca's avatar
Luca committed
76
    
77 78
    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')
ulrich_y's avatar
ulrich_y committed
79

80 81 82 83
    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')
Luca's avatar
Luca committed
84
  end subroutine do_MPL_tests
Luca's avatar
Luca committed
85

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

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

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

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

  subroutine do_GPL_tests()
ulrich_y's avatar
ulrich_y committed
106
    complex(kind=prec) :: ref, res
Luca Naterop's avatar
Luca Naterop committed
107
    real(kind=prec) :: z, xchen
Luca Naterop's avatar
Luca Naterop committed
108
    print*, 'doing GPL tests...'
Luca's avatar
Luca committed
109
    
110 111
    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')
Luca's avatar
Luca committed
112
    
113 114 115 116 117
    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')
Luca's avatar
Luca committed
118

Luca Naterop's avatar
Luca Naterop committed
119
    ! requires making convergent
120 121
    ref = (0.0959304167763934268889451723647597269243_prec,-0.882935179519785044294098901083860819793_prec)
    call test_one_flat(cmplx([0,1,3,2],kind=prec),ref,'2.4')
Luca Naterop's avatar
Luca Naterop committed
122

123 124
    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')
Luca Naterop's avatar
Luca Naterop committed
125 126

    ! requires hoelder convolution
127 128
    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')
Luca Naterop's avatar
Luca Naterop committed
129 130

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

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

135 136
    ref = (-0.00501254182354428204309373895836778138658_prec,0.0)
    call test_one_flat(cmplx([1./z,1.0_prec],kind=prec),ref,'3.1')
Luca Naterop's avatar
Luca Naterop committed
137

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

141 142
    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')
Luca Naterop's avatar
Luca Naterop committed
143

144 145
    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')
Luca Naterop's avatar
Luca Naterop committed
146
    
147 148
    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')
Luca Naterop's avatar
Luca Naterop committed
149
    
150 151
    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')
Luca Naterop's avatar
Luca Naterop committed
152
    
153 154
    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')
155 156 157 158 159


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

160 161 162
    z = 1/100._prec
    ref = (-1.203972804325935992622746217761838502954_prec,0.0)
    call test_one_flat((/ (0._prec,0.), (0.3_prec,0.) /),ref,'4.1')
163

164 165
    ref = (10.6037962209567960211233327771880353832_prec,0.0)
    call test_one_flat( (/ (0._prec,0.), (0._prec,0.), (0.01_prec,0.) /),ref,'4.2')
166

167 168
    ref = (5.0429900651807654463744449762857296353E-4_prec,0.0)
    call test_one_flat((/(100._prec,0.), (1._prec,0.), (0.3_prec,0.) /),ref,'4.3')
169

170 171
    ref = (0.05630861877185110944118569266547272128_prec,0.0)
    call test_one_flat( (/ (1._prec,0.), (0._prec,0.), (0.01_prec,0.) /) ,ref,'4.4')
172

173 174
    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')
175

176 177
    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')
178

179 180
    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')
181

182 183
    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')
184

185 186
    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')
187

188 189
    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')
190

191 192
    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')
193

ulrich_y's avatar
ulrich_y committed
194
    ! Increase test coverage
195 196
    ref = (-1.203972804325935992622746217761838502953611_prec, pi)
    call test_one_flat((/ (0._prec,0.), (-0.3_prec,0.) /),ref,'5.1')
ulrich_y's avatar
ulrich_y committed
197

198 199 200 201
    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.)
ulrich_y's avatar
ulrich_y committed
202
    print*, '  ', 'testing GPL ', '5.3', ' ...'
203
    res = G((/100._prec, 100._prec, 1._prec, 0._prec, 1._prec/))
ulrich_y's avatar
ulrich_y committed
204 205 206 207 208 209 210 211 212 213 214
    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', ' ...'
215
    res = G((/100._prec,100._prec,1._prec,0._prec/), 1._prec)
ulrich_y's avatar
ulrich_y committed
216 217 218
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.7', ' ...'
219
    res = G((/(100._prec,0.),(100._prec,0.),(1._prec,0.),(0._prec,0.)/), (1._prec,0.))
ulrich_y's avatar
ulrich_y committed
220 221 222
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.8', ' ...'
223
    res = G( (/1,1,1,1/) , (/ 100._prec, 100._prec, 1._prec, 0._prec /), 1._prec)
ulrich_y's avatar
ulrich_y committed
224 225 226
    call check(res,ref)

    !ref = cmplx(0.01592795952537145)
227
    ref = (0.0485035531168477720320998551314479928648_prec,-2.82326885118931135859594797186107474877_prec)
ulrich_y's avatar
ulrich_y committed
228
    print*, '  ', 'testing GPL ', '5.9', ' ...'
229
    res = G( (/ 3, 2 /), toinum((/ 1.3_prec, 1.1_prec /)), toinum(4._prec) )
ulrich_y's avatar
ulrich_y committed
230 231 232
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.10', ' ...'
233
    res = G( (/ 3, 2 /), (/ 1.3_prec, 1.1_prec /), 4._prec )
ulrich_y's avatar
ulrich_y committed
234 235 236
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.11', ' ...'
237 238 239 240 241 242
    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 /))
ulrich_y's avatar
ulrich_y committed
243 244
    call check(res,ref)

Luca Naterop's avatar
Luca Naterop committed
245

Luca Naterop's avatar
Luca Naterop committed
246 247 248

  end subroutine do_GPL_tests
  
ulrich_y's avatar
ulrich_y committed
249 250 251



ulrich_y's avatar
ulrich_y committed
252
#ifdef HAVE_GINAC
253
#ifdef HAVE_MM
ulrich_y's avatar
ulrich_y committed
254 255 256 257 258 259 260 261 262

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

    evalt =0.

    do i=1,size(arr)
263
      if (abs(arr(i)) .gt. 1.e10) then ! isnan?
ulrich_y's avatar
ulrich_y committed
264 265 266 267 268 269 270 271 272 273 274 275 276
        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
277
    endif
ulrich_y's avatar
ulrich_y committed
278 279 280
  end function

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

    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
302 303 304
    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.)) )
ulrich_y's avatar
ulrich_y committed
305

ulrich_y's avatar
ulrich_y committed
306
  end subroutine
ulrich_y's avatar
ulrich_y committed
307

ulrich_y's avatar
ulrich_y committed
308 309 310 311
  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
312 313
    integer(kind=8) cstart, cend, count_rate
    real(kind=prec) :: time(2), ttime(2)
ulrich_y's avatar
ulrich_y committed
314 315 316 317 318
    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
319
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
320 321 322
      do i=1,size(args,3)
        res=evalt(args(j,:,i),0)
      enddo
ulrich_y's avatar
ulrich_y committed
323
      call system_clock(cend, count_rate=count_rate)
324
      time(1) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
325
      if (time(1).lt.zero) print*,j, cend-cstart,count_rate
ulrich_y's avatar
ulrich_y committed
326 327
      ttime(1) = ttime(1) + time(1)

ulrich_y's avatar
ulrich_y committed
328
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
329 330 331
      do i=1,size(args,3)
        res=evalt(args(j,:,i),1)
      enddo
ulrich_y's avatar
ulrich_y committed
332
      call system_clock(cend, count_rate=count_rate)
333
      time(2) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
334 335 336 337 338 339 340 341
      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*,
342
    print*,ttime
ulrich_y's avatar
ulrich_y committed
343 344
    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
345 346 347 348 349 350
    ! 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)
351
    ttime(1) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
352 353 354 355 356 357

    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)
358
    ttime(2) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
359 360 361

    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
362
900 FORMAT(a , 'Function ',i4,'/',i4,' for ',a)
363
901 format('Evaluating ',A,' using GiNaC at ',F9.2,'kG/s and handyG at ',F9.2,'kG/s (',I3,'x)')
ulrich_y's avatar
ulrich_y committed
364 365 366 367
  end subroutine

  subroutine do_timing_tests(n)
    use gtestchen   , only: inichen   =>args
ulrich_y's avatar
ulrich_y committed
368
    use gtestchenff , only: inichenff =>args
ulrich_y's avatar
ulrich_y committed
369 370 371 372 373 374
    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
375
    complex(kind=prec) :: fargs(  540,5,n)
ulrich_y's avatar
ulrich_y committed
376 377 378 379 380 381 382 383 384
    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
385 386
      cargs(:,:,i) = inichen  (cmplx(x,kind=prec), cmplx(z,kind=prec))
      fargs(:,:,i) = inichenff(cmplx(x,kind=prec), cmplx(z,kind=prec))
ulrich_y's avatar
ulrich_y committed
387 388 389

      w = ran2(ranseed) ! 0<w<1
      z = ran2(ranseed) * (sqrt(1-w+w**2)-sqrt(w)) + sqrt(w)
390
      nargs(:,:,i) = inimuonenp(cmplx(w,kind=prec), cmplx(z,kind=prec))
ulrich_y's avatar
ulrich_y committed
391 392 393

      x = ran2(ranseed)
      y = ran2(ranseed)
394
      pargs(:,:,i) = inimuone(cmplx(x,kind=prec), cmplx(y,kind=prec))
ulrich_y's avatar
ulrich_y committed
395 396
    enddo

ulrich_y's avatar
ulrich_y committed
397 398
    cargs(1181,:,:)=1.e15
    fargs(367,:,:)=1.e15
ulrich_y's avatar
ulrich_y committed
399 400 401


    open(unit=9, file="stats.txt")
ulrich_y's avatar
ulrich_y committed
402 403
    write(9,*) "Chen form factor"
    call do_one_speed_test(fargs,9,"Chen FF")
ulrich_y's avatar
ulrich_y committed
404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431
    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
432

433
#endif
ulrich_y's avatar
ulrich_y committed
434
#endif
435 436 437
  ! subroutine do_shuffle_tests() 
  !   complex(kind=prec) :: v(2) = cmplx((/1,2/))
  !   complex(kind=prec) :: w(2) = cmplx((/3,4/))
438

439 440
  !   call print_matrix(shuffle_product(v,w))
  ! end subroutine do_shuffle_tests
441

Luca's avatar
Luca committed
442
END PROGRAM TEST
Luca's avatar
Luca committed
443 444

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