test.f90 18.3 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

ulrich_y's avatar
ulrich_y committed
21
  tol = zero * 1.e5_prec
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
    
ulrich_y's avatar
ulrich_y committed
71 72 73 74 75
    ref = (0.0226974036143118403367048874190211851348_prec,0.)
    call test_one_MPL((/ 1,1 /),cmplx([ 119._prec/377._prec,  35._prec / 102._prec ],kind=prec),ref, '1.1')

    ref = (2.31346156303083313941538625891818239608127E-4_prec,0.)
    call test_one_MPL((/ 1,1 /),cmplx([0.03_prec, 100._prec-30._prec*sqrt(11._prec) ], kind=prec),ref, '1.2')
Luca's avatar
Luca committed
76
    
ulrich_y's avatar
ulrich_y committed
77 78
    ref = (2.34617486908178608355021488772249510275409004880846912378356e-5_prec,0.)
    call test_one_MPL((/ 2,1,2 /),cmplx([0.03_prec, 100._prec-30._prec*sqrt(11._prec), 250._prec/9._prec + sqrt(6875._prec)/3._prec ], kind=prec),ref, '1.3')
Luca's avatar
Luca committed
79
    
80 81
    ref = (-0.06565799418838388900030511939327151809905663800733251335678_prec,0.)
    call test_one_MPL((/1, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.4')
ulrich_y's avatar
ulrich_y committed
82

83
    ref = (-0.0319989639656484125354488353677231989818697987217232462808_prec,0.)
ulrich_y's avatar
ulrich_y committed
84
    call test_one_MPL((/2, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.5')
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
    
111 112
    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
113
    
114 115 116 117 118
    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
119

Luca Naterop's avatar
Luca Naterop committed
120
    ! requires making convergent
121 122
    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
123

124 125
    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
126 127

    ! requires hoelder convolution
128 129
    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
130 131

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

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

136 137
    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
138

139 140
    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
141

142 143
    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
144

145 146
    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
147
    
148 149
    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
150
    
151 152
    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
153
    
154 155
    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')
156 157 158 159 160


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

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

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

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

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

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

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

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

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

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

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

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

ulrich_y's avatar
ulrich_y committed
195
    ! Increase test coverage
196 197
    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
198

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

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

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

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

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

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

ulrich_y's avatar
ulrich_y committed
246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296
    print*, '  ', 'testing GPL ', '5.13', ' ...'
    ref = (-1.229846637025798984989235266565630713100_prec,-1.05705870670612913651142424798975000765_prec)
    res = G([inum(2._prec, +1), inum(7._prec, +1)], inum(5._prec, +1))
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.14', ' ...'
    ref = (-1.229846637025798984989235266565630713100_prec,+1.05705870670612913651142424798975000765_prec)
    res = G([inum(2._prec, -1), inum(7._prec, +1)], inum(5._prec, +1))
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.15', ' ...'
    ref = (0.2982254208688088675254638762780704094718_prec,0.)
    res = G([inum(2._prec, -1), inum(7._prec, +1)], inum(-5._prec, -1))
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.17', ' ...'
    ref = (0.190800137777535619036913153766083992418_prec, 0.)
    res = G((/ 6, 1, 1 /))
    call check(res,ref)

    print*, '  ', 'testing GPL ', '5.18', ' ...'
    ref = (0.058192342415778512650117048874978455691_prec, 0.)
    res = G((/ 6, 1, -1 /))
    call check(res,ref)

    ref = cmplx(log(0.5_prec), pi,kind=prec)
    print*, '  ', 'testing GPL ', '5.19a', ' ...'
    res = G((/ inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
    print*, '  ', 'testing GPL ', '5.19b', ' ...'
    res = G((/ inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
    print*, '  ', 'testing GPL ', '5.19c', ' ...'
    res = G((/ inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
    print*, '  ', 'testing GPL ', '5.19d', ' ...'
    res = G((/ inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

    print*, '  ', 'testing GPL ', '5.19e', ' ...'
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
    print*, '  ', 'testing GPL ', '5.19f', ' ...'
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
    print*, '  ', 'testing GPL ', '5.19g', ' ...'
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
    print*, '  ', 'testing GPL ', '5.19h', ' ...'
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
Luca Naterop's avatar
Luca Naterop committed
297

Luca Naterop's avatar
Luca Naterop committed
298

ulrich_y's avatar
ulrich_y committed
299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324
    ref = -(2.374395270272480200677499763071638424_prec, - 1.273806204919600530933131685580471698_prec)
    print*, '  ', 'testing GPL ', '5.20a', ' ...'
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
    print*, '  ', 'testing GPL ', '5.20b', ' ...'
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
    print*, '  ', 'testing GPL ', '5.20c', ' ...'
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
    print*, '  ', 'testing GPL ', '5.20d', ' ...'
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

    print*, '  ', 'testing GPL ', '5.20e', ' ...'
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
    print*, '  ', 'testing GPL ', '5.20f', ' ...'
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
    print*, '  ', 'testing GPL ', '5.20g', ' ...'
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
    print*, '  ', 'testing GPL ', '5.20h', ' ...'
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
Luca Naterop's avatar
Luca Naterop committed
325 326
  end subroutine do_GPL_tests
  
ulrich_y's avatar
ulrich_y committed
327 328 329



ulrich_y's avatar
ulrich_y committed
330
#ifdef HAVE_GINAC
331
#ifdef HAVE_MM
ulrich_y's avatar
ulrich_y committed
332 333 334 335 336 337 338 339 340

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

    evalt =0.

    do i=1,size(arr)
341
      if (abs(arr(i)) .gt. 1.e10) then ! isnan?
ulrich_y's avatar
ulrich_y committed
342 343 344 345 346 347 348 349 350 351 352 353 354
        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
355
    endif
ulrich_y's avatar
ulrich_y committed
356 357 358
  end function

  subroutine perform_ginacv(n, args)
ulrich_y's avatar
ulrich_y committed
359
    use maths_functions, only:clearcache
ulrich_y's avatar
ulrich_y committed
360 361
    complex(kind=prec) :: args(:,:)
    integer i,n
ulrich_y's avatar
ulrich_y committed
362
    call clearcache
ulrich_y's avatar
ulrich_y committed
363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379

    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
380 381 382
    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
383

ulrich_y's avatar
ulrich_y committed
384
  end subroutine
ulrich_y's avatar
ulrich_y committed
385

ulrich_y's avatar
ulrich_y committed
386 387 388 389
  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
390 391
    integer(kind=8) cstart, cend, count_rate
    real(kind=prec) :: time(2), ttime(2)
ulrich_y's avatar
ulrich_y committed
392 393 394 395 396
    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
397
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
398 399 400
      do i=1,size(args,3)
        res=evalt(args(j,:,i),0)
      enddo
ulrich_y's avatar
ulrich_y committed
401
      call system_clock(cend, count_rate=count_rate)
402
      time(1) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
403
      if (time(1).lt.zero) print*,j, cend-cstart,count_rate
ulrich_y's avatar
ulrich_y committed
404 405
      ttime(1) = ttime(1) + time(1)

ulrich_y's avatar
ulrich_y committed
406
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
407 408 409
      do i=1,size(args,3)
        res=evalt(args(j,:,i),1)
      enddo
ulrich_y's avatar
ulrich_y committed
410
      call system_clock(cend, count_rate=count_rate)
411
      time(2) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
412 413 414 415 416 417 418 419
      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*,
420
    print*,ttime
ulrich_y's avatar
ulrich_y committed
421 422
    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
423 424 425 426 427 428
    ! 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)
429
    ttime(1) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
430 431 432 433 434 435

    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)
436
    ttime(2) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
437 438 439

    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
440
900 FORMAT(a , 'Function ',i4,'/',i4,' for ',a)
441
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
442 443 444 445
  end subroutine

  subroutine do_timing_tests(n)
    use gtestchen   , only: inichen   =>args
ulrich_y's avatar
ulrich_y committed
446
    use gtestchenff , only: inichenff =>args
ulrich_y's avatar
ulrich_y committed
447 448 449 450 451 452
    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
453
    complex(kind=prec) :: fargs(  540,5,n)
ulrich_y's avatar
ulrich_y committed
454 455 456 457 458 459 460 461 462
    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
463 464
      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
465 466 467

      w = ran2(ranseed) ! 0<w<1
      z = ran2(ranseed) * (sqrt(1-w+w**2)-sqrt(w)) + sqrt(w)
468
      nargs(:,:,i) = inimuonenp(cmplx(w,kind=prec), cmplx(z,kind=prec))
ulrich_y's avatar
ulrich_y committed
469 470 471

      x = ran2(ranseed)
      y = ran2(ranseed)
472
      pargs(:,:,i) = inimuone(cmplx(x,kind=prec), cmplx(y,kind=prec))
ulrich_y's avatar
ulrich_y committed
473 474
    enddo

ulrich_y's avatar
ulrich_y committed
475 476
    cargs(1181,:,:)=1.e15
    fargs(367,:,:)=1.e15
ulrich_y's avatar
ulrich_y committed
477 478 479


    open(unit=9, file="stats.txt")
ulrich_y's avatar
ulrich_y committed
480 481
    write(9,*) "Chen form factor"
    call do_one_speed_test(fargs,9,"Chen FF")
ulrich_y's avatar
ulrich_y committed
482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509
    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
510

511
#endif
ulrich_y's avatar
ulrich_y committed
512
#endif
513 514 515
  ! subroutine do_shuffle_tests() 
  !   complex(kind=prec) :: v(2) = cmplx((/1,2/))
  !   complex(kind=prec) :: w(2) = cmplx((/3,4/))
516

517 518
  !   call print_matrix(shuffle_product(v,w))
  ! end subroutine do_shuffle_tests
519

Luca's avatar
Luca committed
520
END PROGRAM TEST
Luca's avatar
Luca committed
521 522

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