JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
democache.f90
Go to the documentation of this file.
1 !!
2 !! File democache.f90 is part of COLLIER
3 !! - A Complex One-Loop Library In Extended Regularizations
4 !!
5 !! Copyright (C) 2015, 2016 Ansgar Denner, Stefan Dittmaier, Lars Hofer
6 !!
7 !! COLLIER is licenced under the GNU GPL version 3, see COPYING for details.
8 !!
9 
10 program democache
11 
12  use collier
13 
14  implicit none
15 
16  double complex :: mominv6(15), masses2_6(0:5), masses2_6b(0:5)
17  double complex :: mominv5(10), masses2_5(0:4)
18  double complex :: mominv4(6), masses2_4(0:3)
19  double complex :: mominv3(3), masses2_3(0:2)
20  double complex :: mominv2(1), masses2_2(0:1)
21  double complex, allocatable :: t1coeff(:),t1coeffuv(:)
22  double complex, allocatable :: t2coeff(:),t2coeffuv(:)
23  double complex, allocatable :: t3coeff(:),t3coeffuv(:)
24  double complex, allocatable :: t4coeff(:),t4coeffuv(:)
25  double complex, allocatable :: t5coeff(:),t5coeffuv(:)
26  double complex, allocatable :: t6coeff(:),t6coeffuv(:)
27  double precision, allocatable :: t1err(:),t2err(:),t3err(:),t4err(:)
28  double precision, allocatable :: t5err(:),t6err(:)
29  real :: starttime,endtime
30  integer :: rank2,rank3,rank4,rank5,rank6,nmax,rmax
31  integer :: i,k,nevents,nloops,icache,mycase
32 
33 ! Nmax = maximal degree N of N-point function
34  nmax = 6
35  rmax = 4
36 
37 ! Initialize COLLIER
38 ! the default directory for diagnostic messages is "output_cll"
39  call init_cll(nmax,rmax,"output_cll")
40 
41 ! INITIALIZE CACHE nr. "icache" for up to Nmax-point functions
42  icache = 1
43  call initcachesystem_cll(icache,nmax)
44 
45 
46 ! allocate fields for coefficients
47  rank2=3
48  rank3=4
49  rank4=4
50  rank5=4
51  rank6=4
52  allocate(t2coeff(getnc_cll(2,rank2)))
53  allocate(t2coeffuv(getnc_cll(2,rank2)))
54  allocate(t2err(0:rank2))
55  allocate(t3coeff(getnc_cll(3,rank3)))
56  allocate(t3coeffuv(getnc_cll(3,rank3)))
57  allocate(t3err(0:rank3))
58  allocate(t4coeff(getnc_cll(4,rank4)))
59  allocate(t4coeffuv(getnc_cll(4,rank4)))
60  allocate(t4err(0:rank4))
61  allocate(t5coeff(getnc_cll(5,rank5)))
62  allocate(t5coeffuv(getnc_cll(5,rank5)))
63  allocate(t5err(0:rank5))
64  allocate(t6coeff(getnc_cll(6,rank6)))
65  allocate(t6coeffuv(getnc_cll(6,rank6)))
66  allocate(t6err(0:rank6))
67 
68 ! set some values for masses
69  masses2_6(0) = dcmplx(1d2,-1d0)
70  masses2_6(1) = dcmplx(1d2,-1d0)
71  masses2_6(2) = dcmplx(0d0,0d0)
72  masses2_6(3) = dcmplx(1d1,0d0)
73  masses2_6(4) = dcmplx(1d1,0d0)
74  masses2_6(5) = dcmplx(0d0,0d0)
75 
76 
77 ! call InitMonitoring_cll
78 
79 ! issue a bunch of calls of COLLIER for different tensor integrals
80 ! for several "events"
81 
82  nevents = 2000
83 
84  nloops = 2
85 
86  write(*,'(/(a))') ' *** Demonstration of cache system ***'
87  write(*,'(/(a))') ' NOTE: This example merely illustrates the usage of the cache system.'
88  write(*,'((a))') ' The results are not representative for the performance'
89  write(*,'((a))') ' of Collier, neither wrt to speed nor precision.'
90  write(*,'(/a21,i10/)') ' Number of "Events": ',nevents
91  write(*,'(/a39,i4/)') ' Number of "identical calls/integral": ',nloops
92 
93  do mycase=1,4
94 
95  if (mycase.eq.1) then
96 ! use COLI branch
97  call setmode_cll(1)
98 ! initialize cache nr. icache for tensor integrals up to N=6
99  call initcachesystem_cll(icache,6)
100 
101  write(*,'(/(a))') ' COLI branch with cache: '
102 
103  elseif (mycase.eq.2) then
104 ! use COLI branch
105  call setmode_cll(1)
106 ! switch off cache
107  call switchoffcachesystem_cll
108 
109  write(*,'(/(a))') ' COLI branch without cache: '
110 
111  else if (mycase.eq.3) then
112 ! use DD branch
113  call setmode_cll(2)
114 ! reinitialize cache nr. icache
115  call initcachesystem_cll(icache,6)
116 
117  write(*,'(/(a))') ' DD branch with cache: '
118 
119  else if (mycase.eq.4) then
120 ! use DD branch
121  call setmode_cll(2)
122 ! switch off cache
123  call switchoffcachesystem_cll
124 
125  write(*,'(/(a))') ' DD branch without cache '
126 
127  end if
128 
129  call cpu_time(starttime)
130 
131 ! write(*,*) 'start',starttime
132 
133  eventloop: do i=1,nevents
134 
135 ! get momenta and invariants (typically from Monte Carlo generator)
136  call getinvariants(mominv6,i)
137 
138 ! TELL COLLIER START OF NEW EVENT for cache "icache"
139  call initevent_cll(icache)
140 
141 ! call bunch of tensor integrals or coefficients
142 
143 ! loop to generate identical master calls
144  testloop: do k=1,nloops
145 
146  call tn_cll(t6coeff,t6coeffuv,mominv6,masses2_6,6,rank6,t6err)
147 
148 ! pick arguments for tensor integrals of lower rank
149  mominv5 = submominv(6,0,mominv6)
150  masses2_5 = submasses(6,0,masses2_6)
151  call tn_cll(t5coeff,t5coeffuv,mominv5,masses2_5,5,rank5,t5err)
152 
153  mominv4 = submominv(5,0,mominv5)
154  masses2_4 = submasses(5,0,masses2_5)
155  mominv3 = submominv(4,0,mominv4)
156  masses2_3 = submasses(4,0,masses2_4)
157  call tn_cll(t3coeff,t3coeffuv,mominv3,masses2_3,3,rank3,t3err)
158  call tn_cll(t4coeff,t4coeffuv,mominv4,masses2_4,4,rank4,t4err)
159 
160 ! change one of the masses
161  masses2_6b(0) = (0d0,0d0)
162  masses2_6b(1:5) = masses2_6(1:5)
163  call tn_cll(t6coeff,t6coeffuv,mominv6,masses2_6b,6,rank6,t6err)
164 
165  mominv5 = submominv(6,1,mominv6)
166  masses2_5 = submasses(6,1,masses2_6b)
167  mominv4 = submominv(5,2,mominv5)
168  masses2_4 = submasses(5,2,masses2_5)
169  mominv3 = submominv(4,0,mominv4)
170  masses2_3 = submasses(4,0,masses2_4)
171  mominv2 = submominv(3,0,mominv3)
172  masses2_2 = submasses(3,0,masses2_3)
173  call tn_cll(t2coeff,t2coeffuv,mominv2,masses2_2,2,rank2,t2err)
174  call tn_cll(t5coeff,t5coeffuv,mominv5,masses2_5,5,rank5,t5err)
175  call tn_cll(t4coeff,t4coeffuv,mominv4,masses2_4,4,rank4,t4err)
176 
177  end do testloop
178 
179  end do eventloop
180 
181  call cpu_time(endtime)
182 
183  write(*,'(/a29,f10.1/)') ' Needed CPU time in seconds: ',endtime-starttime
184 
185  end do
186 
187 ! call PrintStatistics_cll
188 
189  deallocate(t2coeff,t2coeffuv,t2err)
190  deallocate(t3coeff,t3coeffuv,t3err)
191  deallocate(t4coeff,t4coeffuv,t4err)
192  deallocate(t5coeff,t5coeffuv,t5err)
193  deallocate(t6coeff,t6coeffuv,t6err)
194 
195 contains
196 
197  subroutine getinvariants(MomInv6,ievent)
198 
199  implicit none
200 
201  double precision p(1:6,0:3)
202  double complex MomInv6(1:15)
203  double precision :: s(0:5,0:5),fact=47d0/49d0
204  integer :: N=6,k,i,newmom,ievent
205  save s,fact,n
206 
207  newmom=0
208 
209  if (ievent.eq.1) then
210  newmom=1
211 ! use momentum set from Racoon
212  p(1,0) = -0.9500000000000000d+02
213  p(1,1) = 0.0000000000000000d+00
214  p(1,2) = 0.0000000000000000d+00
215  p(1,3) = 0.9500000000000000d+02
216 
217  p(2,0) = -0.9500000000000000d+02
218  p(2,1) = 0.0000000000000000d+00
219  p(2,2) = 0.0000000000000000d+00
220  p(2,3) = -0.9500000000000000d+02
221 
222  p(3,0) = 0.2046111001757171d+02
223  p(3,1) = 0.1057734233089455d+02
224  p(3,2) = -0.2324961261504543d+01
225  p(3,3) = 0.1736005205921753d+02
226 
227  p(4,0) = 0.3558305163378869d+01
228  p(4,1) = 0.1436222934374051d+01
229  p(4,2) = -0.2174258125294355d+01
230  p(4,3) = -0.2423097382091398d+01
231 
232  p(5,0) = 0.8154540918019539d+02
233  p(5,1) = -0.5230395944682889d+02
234  p(5,2) = 0.3083642435466509d+02
235  p(5,3) = 0.5443403822581044d+02
236 
237  p(6,0) = 0.8443517563885433d+02
238  p(6,1) = 0.4029039418156027d+02
239  p(6,2) = -0.2633720496786619d+02
240  p(6,3) = -0.6937099290293661d+02
241 
242 
243  else if (ievent.eq.50) then
244  newmom=1
245 ! set with small momentum 1.9 versus 3.0 seconds
246  p(1,0) = -250.0000000000000d0
247  p(1,1) = 0.0000000000000000d+00
248  p(1,2) = 0.0000000000000000d+00
249  p(1,3) = -250.0000000000000d0
250 
251  p(2,0) = -250.0000000000000d0
252  p(2,1) = 0.0000000000000000d+00
253  p(2,2) = 0.0000000000000000d+00
254  p(2,3) = 250.0000000000000d0
255 
256  p(3,0) = 239.91697591218110d0
257  p(3,1) = 31.640107125154238d0
258  p(3,2) = -116.12928664566658d0
259  p(3,3) = -207.54047252312847d0
260 
261  p(4,0) = 10.276659578244022d0
262  p(4,1) = 1.2584629023924521d0
263  p(4,2) = 10.131487190720421d0
264  p(4,3) = -1.1742957526504370d0
265 
266  p(5,0) = 181.18800109198216d0
267  p(5,1) = -35.098898474333737d0
268  p(5,2) = 109.26326872331617d0
269  p(5,3) = 140.20947604742324d0
270 
271  p(6,0) = 68.618363417592718d0
272  p(6,1) = 2.2003284467870472d0
273  p(6,2) = -3.2654692683700119d0
274  p(6,3) = 68.505292228355671d0
275 
276  else if (ievent.eq.100) then
277  newmom=1
278  p(1,0) = -100.0000000000000d0
279  p(1,1) = 0.000000000000000d0
280  p(1,2) = 0.000000000000000d0
281  p(1,3) = -100.0000000000000d0
282 
283  p(2,0) = -100.0000000000000d0
284  p(2,1) = 0.000000000000000d0
285  p(2,2) = 0.000000000000000d0
286  p(2,3) = 100.0000000000000d0
287 
288  p(5,0) = 42.13045937773281d0
289  p(5,1) = 38.71045260055064d0
290  p(5,2) = 16.54764683624944d0
291  p(5,3) = 1.628450497276325d0
292 
293  p(6,0) = 73.31699208541762d0
294  p(6,1) = -66.17895647220249d0
295  p(6,2) = 31.53024134467396d0
296  p(6,3) = -1.253367244715983d0
297 
298  p(4,0) = 69.60509150033207d0
299  p(4,1) = 34.42307953372518d0
300  p(4,2) = -60.31848303205667d0
301  p(4,3) = 4.647683605402023d0
302 
303  p(3,0) = 14.94745703651754d0
304  p(3,1) = -6.954575662073326d0
305  p(3,2) = 12.24059485113326d0
306  p(3,3) = -5.022766857962365d0
307  else if (ievent.eq.150) then
308  newmom=1
309  p(1,0) = -100.0000000000000d0
310  p(1,1) = 0.0000000000000000d+00
311  p(1,2) = 0.0000000000000000d+00
312  p(1,3) = -100.0000000000000d0
313 
314  p(2,0) = -100.0000000000000d0
315  p(2,1) = 0.0000000000000000d+00
316  p(2,2) = 0.0000000000000000d+00
317  p(2,3) = 100.0000000000000d0
318 
319  p(3,0) = 72.02905523601244d0
320  p(3,1) = 32.33193950844802d0
321  p(3,2) = 64.10893732125734d0
322  p(3,3) = 5.733641195059672d0
323 
324  p(4,0) = 25.07519219892066d0
325  p(4,1) = 2.906699172283076d0
326  p(4,2) = -24.65281899347465d0
327  p(4,3) = 3.543286526606821d0
328 
329  p(5,0) = 60.19609558319355d0
330  p(5,1) = 2.282219923475870d0
331  p(5,2) = -59.72688022443283d0
332  p(5,3) = -7.145710205299344d0
333 
334  p(6,0) = 42.69965698187335d0
335  p(6,1) = -37.52085860420697d0
336  p(6,2) = 20.27076189665014d0
337  p(6,3) = -2.131217516367149d0
338  else if (ievent.eq.200) then
339  newmom=1
340  p(1,0) = -100.0000000000000d0
341  p(1,1) = 0.0000000000000000d+00
342  p(1,2) = 0.0000000000000000d+00
343  p(1,3) = -100.0000000000000d0
344 
345  p(2,0) = -100.0000000000000d0
346  p(2,1) = 0.0000000000000000d+00
347  p(2,2) = 0.0000000000000000d+00
348  p(2,3) = 100.0000000000000d0
349 
350  p(3,0) = 69.63925440091566d0
351  p(3,1) = -68.33683461638704d0
352  p(3,2) = -3.803463998624776d0
353  p(3,3) = -12.85443307731553d0
354 
355  p(4,0) = 34.49634785978589d0
356  p(4,1) = 21.54290146454439d0
357  p(4,2) = -15.60187514686951d0
358  p(4,3) = 21.96549348532329d0
359 
360  p(5,0) = 66.20628864972379d0
361  p(5,1) = 65.60907499667245d0
362  p(5,2) = -3.500460884404577d0
363  p(5,3) = -8.152834381263276d0
364 
365  p(6,0) = 29.65810908957467d0
366  p(6,1) = -18.81514184482980d0
367  p(6,2) = 22.90580002989887d0
368  p(6,3) = -0.9582260267444891d0
369 
370  else if (ievent.eq.250) then
371  newmom=1
372  p( 1, 0) = -159.2798637127597488d0
373  p( 1, 1) = 0.0000000000000000d0
374  p( 1, 2) = 0.0000000000000000d0
375  p( 1, 3) = -159.2798637127597488d0
376  p( 2, 0) = -237.6261878491357891d0
377  p( 2, 1) = 0.0000000000000000d0
378  p( 2, 2) = 0.0000000000000000d0
379  p( 2, 3) = 237.6261878491357891d0
380  p( 3, 0) = 133.6018149124976446d0
381  p( 3, 1) = 60.5134748218625731d0
382  p( 3, 2) = 112.6312302087531094d0
383  p( 3, 3) = -38.7526811273741174d0
384  p( 4, 0) = 69.9256648527542950d0
385  p( 4, 1) = 66.6463990867762135d0
386  p( 4, 2) = 7.6538437283820429d0
387  p( 4, 3) = -19.7300473909055754d0
388  p( 5, 0) = 85.7314709474216841d0
389  p( 5, 1) = -38.6625297933730749d0
390  p( 5, 2) = -63.4432116127613241d0
391  p( 5, 3) = -42.7791164126088006d0
392  p( 6, 0) = 107.6471008492219283d0
393  p( 6, 1) = -88.4973441152657188d0
394  p( 6, 2) = -56.8418623243738281d0
395  p( 6, 3) = 22.9155207945124531d0
396 
397 
398  end if
399 
400  if (newmom.eq.1) then
401  ! define invariants
402  n=6
403  do k=0,n-1
404  s(modulo(k+1,n),k) = p(k+1,0)**2
405  do i=1,3
406  s(modulo(k+1,n),k) = s(modulo(k+1,n),k) - p(k+1,i)**2
407  end do
408  if(abs(s(modulo(k+1,n),k)).lt.1d-14* abs(p(k+1,0))**2) s(modulo(k+1,n),k) = 0d0
409 
410  s(k,modulo(k+1,n)) = s(modulo(k+1,n),k)
411  end do
412 
413  do k=0,n-1
414  s(modulo(k+2,n),k) = (p(modulo(k+1,n)+1,0) + p(k+1,0))**2
415  do i=1,3
416  s(modulo(k+2,n),k) = s(modulo(k+2,n),k) &
417  - (p(modulo(k+1,n)+1,i) + p(k+1,i))**2
418  end do
419  s(k,modulo(k+2,n)) = s(modulo(k+2,n),k)
420  end do
421  do k=0,(n-1)/2
422  s(modulo(k+3,n),k) = (p(modulo(k+2,n)+1,0) + p(modulo(k+1,n)+1,0) + p(k+1,0))**2
423  do i=1,3
424  s(modulo(k+3,n),k) = s(modulo(k+3,n),k) &
425  - (p(modulo(k+2,n)+1,i) + p(modulo(k+1,n)+1,i) + p(k+1,i))**2
426  end do
427  s(k,modulo(k+3,n)) = s(modulo(k+3,n),k)
428  end do
429 
430  do k=1,n
431  mominv6(k) = s(modulo(k,n),k-1)
432  end do
433  do k=1,n
434  mominv6(k+n) = s(modulo(k+1,n),k-1)
435  end do
436  do k=1,n/2
437  mominv6(k+2*n) = s(modulo(k+2,n),k-1)
438  end do
439  fact = 47d0/97d0
440 
441  else
442 ! rescale invariants
443  fact = 4*fact*(1-fact)
444 
445  do k=1,n
446  mominv6(k) = s(modulo(k,n),k-1)*fact
447  end do
448  do k=1,n
449  mominv6(k+n) = s(modulo(k+1,n),k-1)*fact
450  end do
451  do k=1,n/2
452  mominv6(k+2*n) = s(modulo(k+2,n),k-1)*fact
453  end do
454  end if
455 
456  end subroutine getinvariants
457 
458 
459 
460 
461 
462  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
463  ! function SubMasses
464  !
465  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
466 
467  function submasses(N,k,masses)
468 
469  integer :: n,k,i
470  double complex :: masses(0:n-1), submasses(0:n-2)
471 
472  do i=0,n-2
473  if (i.lt.k) then
474  submasses(i) = masses(i)
475  else
476  submasses(i) = masses(i+1)
477  end if
478  end do
479 
480  end function
481 
482 
483 
484 
485 
486  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
487  ! function SubMomInv
488  ! ******************************************
489  ! * rewritten by Ansgar Denner 25.03.15 *
490  ! ******************************************
491  !
492  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
493 
494  function submominv(N,k,MomInv)
495 
496  integer, intent(in) :: n,k
497  integer :: i,cnt,moms,limit
498  double complex :: mominv(n*(n-1)/2), submominv((n-1)*(n-2)/2)
499 
500  if (n.eq.2) return
501 
502  cnt = 1
503  do moms=1,(n-3)/2
504  do i=1,n
505  if (i.ne.k+1) then
506  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
507  submominv(cnt) = mominv(n*moms+i)
508  else
509  submominv(cnt) = mominv(n*(moms-1)+i)
510  end if
511 
512  cnt = cnt+1
513  end if
514  end do
515  end do
516 
517  if (mod(n,2).eq.1) then
518  moms = (n-1)/2
519  if(k.le.(n-1)/2) then
520  limit = (n+1)/2
521  else
522  limit = (n-1)/2
523  end if
524  do i=1,limit
525  if (i.ne.k+1) then
526  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
527  submominv(cnt) = mominv(n*moms+i-(n-1)/2)
528  else
529  submominv(cnt) = mominv(n*(moms-1)+i)
530  end if
531 
532  cnt = cnt+1
533  end if
534  end do
535 
536  else
537  do i=1,n
538  if (i.ne.k+1) then
539  moms = (n/2-1)
540  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
541  if (i.gt.n/2) then
542  submominv(cnt) = mominv(n*moms-n/2+i)
543  else
544  submominv(cnt) = mominv(n*moms+i)
545  endif
546  else
547  submominv(cnt) = mominv(n*(moms-1)+i)
548  end if
549 
550  cnt = cnt+1
551  end if
552  end do
553 
554  end if
555 
556 
557  end function
558 
559 
560 end program democache
submasses
double complex function, dimension(0:n-2) submasses(N, k, masses)
Definition: democache.f90:468
collier
Definition: COLLIER.F90:30
submominv
double complex function, dimension((n-1) *(n-2)/2) submominv(N, k, MomInv)
Definition: democache.f90:495
democache
program democache
Definition: democache.f90:10
getinvariants
subroutine getinvariants(N, MomInv, MomVec)
Definition: demo.f90:3262