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
#ifndef NOSPEED
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
ulrich_y's avatar
ulrich_y committed
35
#endif
Luca's avatar
tests  
Luca committed
36 37 38 39 40 41 42

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

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

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

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

Luca's avatar
Luca committed
69
  subroutine do_MPL_tests()
Luca's avatar
Luca committed
70
    complex(kind=prec) :: ref
Luca Naterop's avatar
Luca Naterop committed
71
    print*, 'doing MPL tests...'
Luca's avatar
Luca committed
72
    
ulrich_y's avatar
ulrich_y committed
73 74 75 76 77
    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
78
    
ulrich_y's avatar
ulrich_y committed
79 80
    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
81
    
82 83
    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
84

85
    ref = (-0.0319989639656484125354488353677231989818697987217232462808_prec,0.)
ulrich_y's avatar
ulrich_y committed
86
    call test_one_MPL((/2, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.5')
Luca's avatar
Luca committed
87
  end subroutine do_MPL_tests
Luca's avatar
Luca committed
88

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

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

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

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

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

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

126 127
    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
128 129

    ! requires hoelder convolution
130 131
    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
132 133

    ! here the tests from the mathematica nb start
134 135
    ! --------------------------------------------------------------------------

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

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

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

144 145
    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
146

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


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

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

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

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

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

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

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

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

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

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

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

194 195
    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')
196

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

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

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

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

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

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

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

ulrich_y's avatar
ulrich_y committed
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 297 298
    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
299

Luca Naterop's avatar
Luca Naterop committed
300

ulrich_y's avatar
ulrich_y committed
301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326
    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
327 328
  end subroutine do_GPL_tests
  
ulrich_y's avatar
ulrich_y committed
329 330 331



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

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

    evalt =0.

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

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

    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
382 383 384
    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
385

ulrich_y's avatar
ulrich_y committed
386
  end subroutine
ulrich_y's avatar
ulrich_y committed
387

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

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

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

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

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

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

      x = ran2(ranseed)
      y = ran2(ranseed)
475
      pargs(:,:,i) = inimuone(cmplx(x,kind=prec), cmplx(y,kind=prec))
ulrich_y's avatar
ulrich_y committed
476 477
    enddo

ulrich_y's avatar
ulrich_y committed
478 479
    cargs(1181,:,:)=1.e15
    fargs(367,:,:)=1.e15
ulrich_y's avatar
ulrich_y committed
480 481 482


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

      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
514

515
#endif
ulrich_y's avatar
ulrich_y committed
516
#endif
517 518 519
  ! subroutine do_shuffle_tests() 
  !   complex(kind=prec) :: v(2) = cmplx((/1,2/))
  !   complex(kind=prec) :: w(2) = cmplx((/3,4/))
520

521 522
  !   call print_matrix(shuffle_product(v,w))
  ! end subroutine do_shuffle_tests
523

Luca's avatar
Luca committed
524
END PROGRAM TEST
Luca's avatar
Luca committed
525 526

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