JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
reductionAB.F90
Go to the documentation of this file.
1 !!
2 !! File reductionAB.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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
11 !
12 ! ************************
13 ! * module reductionAB *
14 ! * by Lars Hofer *
15 ! ************************
16 !
17 ! global variables:
18 ! dprec_coli
19 !
20 ! functions and subroutines:
21 ! CalcA, CalcB, CalcDB
22 !
23 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
24 
25 
26 
27 
28 module reductionab
29 
30  use combinatorics
31  use cache
32  use coli_aux2
33  use collier_global
34 
35 
36  implicit none
37 
38 ! double complex, allocatable :: BCDuv_cp(:,:,:,:,:), BCD_cp(:,:,:,:,:)
39 ! double complex, allocatable :: Euv_cp(:,:,:,:,:,:), E_cp(:,:,:,:,:,:)
40 ! double complex, allocatable :: Fuv_cp(:,:,:,:,:,:,:), F_cp(:,:,:,:,:,:,:)
41 ! double complex, allocatable :: TN_cp(:,:,:,:), TNuv_cp(:,:,:,:), TeN_cp(:,:,:,:,:)
42 ! integer, allocatable :: rmax_cp(:), scheme_cp(:)
43 
44 
45 contains
46 
47  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
48  ! subroutine CalcA(A,Auv,m02,rmax,Aerr)
49  !
50  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
51 
52  subroutine calca(A,Auv,m02,rmax,Aerr)
53 
54  integer, intent(in) :: rmax
55  double complex, intent(in) :: m02
56  double complex, intent(out) :: A(0:rmax/2), Auv(0:rmax/2)
57  double precision, optional, intent(out) :: Aerr(0:rmax)
58  double complex :: A0_coli, mm02, elimminf2_coli
59  double complex :: fac, ccc
60  integer :: n0
61 
62  if(present(aerr)) then
63  aerr = 0d0
64  end if
65  auv(0) = m02
66  a(0) = a0_coli(m02)
67 
68  ccc = 2d0
69  mm02 = elimminf2_coli(m02)
70  do n0=1,rmax/2
71  fac = m02/(2*(n0+1))
72  ccc = ccc*fac
73  auv(n0) = fac*auv(n0-1)
74  a(n0) = fac*(a(n0-1)+ccc)
75  end do
76  if(present(aerr)) then
77  aerr = abs(a(0))*acc_def_b
78  end if
79 
80  end subroutine calca
81 
82 
83 
84 
85 
86  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
87  ! subroutine CalcB(B,Buv,p10,m02,m12,rmax,id,Berr)
88  !
89  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
90 
91  subroutine calcb(B,Buv,p10,m02,m12,rmax,id,Berr)
92 
93  integer, intent(in) :: rmax, id
94  double complex, intent(in) :: p10,m02,m12
95  double complex, intent(inout) :: B(0:rmax,0:rmax), Buv(0:rmax,0:rmax)
96  double precision, optional, intent(out) :: Berr(0:rmax)
97  double complex, allocatable :: Baux(:,:), Buvaux(:,:),fct(:)
98  double precision, allocatable :: Berraux(:)
99  double complex :: x(3)
100  integer :: rank,switch,cnt,n0,n1,r
101  logical :: nocalc,wrica
102 
103  if (use_cache_system) then
104  if ((ncache.gt.0).and.(ncache.le.ncache_max)) then
105 ! if (use_cache(ncache).ge.2) then
106  x(1)=p10
107  x(2)=m02
108  x(3)=m12
109  rank = rmax
110  switch = 0
111 
112  allocate(fct(2*ncoefsg(rmax,2)+rmax+1))
113  call readcache(fct,2*ncoefsg(rmax,2)+rmax+1,x,3,1,id,2,rank,nocalc,wrica)
114 
115  if(nocalc)then
116  cnt = 0
117  do r=0,rmax
118  do n0=0,r
119  n1 = r-n0
120 
121  cnt = cnt+1
122  b(n0,n1) = fct(cnt)
123  cnt = cnt+1
124  buv(n0,n1) = fct(cnt)
125 
126  end do
127  cnt = cnt+1
128  if(present(berr))then
129  berr(r) = real(fct(cnt))
130  end if
131  end do
132  return
133  endif
134 
135 
136  if(rank.eq.rmax) then
137 
138  allocate(berraux(0:rank))
139  call calcbred(b,buv,p10,m02,m12,rank,id,berraux)
140  if(present(berr)) berr=berraux
141 
142  if (wrica) then
143  cnt = 0
144  do r=0,rank
145  do n0=0,r
146  n1 = r-n0
147 
148  cnt = cnt+1
149  fct(cnt) = b(n0,n1)
150  cnt = cnt+1
151  fct(cnt) = buv(n0,n1)
152 
153  end do
154  cnt = cnt+1
155  fct(cnt) = berraux(r)
156  end do
157 
158  call writecache(fct,2*ncoefsg(rank,2)+rank+1,id,2,rank)
159 
160  end if
161 
162  return
163 
164 
165  else
166  allocate(baux(0:rank,0:rank))
167  allocate(buvaux(0:rank,0:rank))
168  allocate(berraux(0:rank))
169 
170  call calcbred(baux,buvaux,p10,m02,m12,rank,id,berraux)
171 
172  if (wrica) then
173  cnt = 0
174  deallocate(fct)
175  allocate(fct(2*ncoefsg(rank,2)+rank+1))
176  do r=0,rank
177  do n0=0,r
178  n1 = r-n0
179 
180  cnt = cnt+1
181  fct(cnt) = baux(n0,n1)
182  cnt = cnt+1
183  fct(cnt) = buvaux(n0,n1)
184 
185  end do
186  cnt = cnt+1
187  fct(cnt) = berraux(r)
188  end do
189 
190  call writecache(fct,2*ncoefsg(rank,2)+rank+1,id,2,rank)
191 
192  end if
193 
194  b(0:rmax,0:rmax) = baux(0:rmax,0:rmax)
195  buv(0:rmax,0:rmax) = buvaux(0:rmax,0:rmax)
196  if(present(berr))then
197  berr = berraux(0:rmax)
198  end if
199  return
200 
201  end if
202 ! end if
203  end if
204  end if
205 
206 
207  if(present(berr)) then
208  call calcbred(b,buv,p10,m02,m12,rmax,id,berr)
209  else
210  allocate(berraux(0:rmax))
211  call calcbred(b,buv,p10,m02,m12,rmax,id,berraux)
212  deallocate(berraux)
213  end if
214 
215  end subroutine calcb
216 
217 
218 
219 
220 
221  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
222  ! subroutine CalcBred(B,Buv,p10,m02,m12,rmax,id,Berr)
223  !
224  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
225 
226  subroutine calcbred(B,Buv,p10,m02,m12,rmax,id,Berr)
227 
228  integer, intent(in) :: rmax, id
229  double complex, intent(in) :: p10,m02,m12
230  double complex, intent(inout) :: B(0:rmax,0:rmax), Buv(0:rmax,0:rmax)
231  double precision, intent(out) :: Berr(0:rmax)
232 ! double complex, allocatable :: A(:), Auv(:)
233  double complex :: A(0:rmax-1), Auv(0:rmax-1)
234  double complex :: q2, mm02, mm12, f1
235  double complex :: Bn_coli,elimminf2_coli
236  integer :: rmaxA,n0,n1,r,sgn,k,rid,r0
237  logical :: finarg
238  double complex, parameter :: cd0 = dcmplx(0d0,0d0)
239 
240  r0 = 0
241 
242  ! r=0
243  if (r0.eq.0) then
244  buv(0,0) = 1d0
245  b(0,0) = bn_coli(0,p10,m02,m12)
246  end if
247 
248  ! calculate B(0,n1)
249  sgn = (-1)**r0
250  do n1=max(r0,1),rmax
251  sgn = -sgn
252  buv(0,n1) = sgn/(n1+1d0)
253  b(0,n1) = bn_coli(n1,p10,m02,m12)
254  end do
255 
256 
257  ! calculate B(n0,n1) for n0=/=0
258  if (rmax.ge.1) then
259 
260  call elminf2iv_coli(p10,m02,m12,q2,mm02,mm12,finarg)
261 
262 ! q2 = elimminf2_coli(p10)
263 ! mm12 = elimminf2_coli(m12)
264 ! mm02 = elimminf2_coli(m02)
265 ! finarg=q2.ne.cd0.or.mm12.ne.cd0.or.mm02.ne.cd0
266 
267  ! at least one of the invariants p10,m02,m12 different from zero
268  if (finarg) then
269 
270  ! calculation of A-functions
271  rmaxa = max(rmax-1,0)
272 ! allocate(A(0:rmaxA))
273 ! allocate(Auv(0:rmaxA))
274  call calca(a,auv,mm12,2*rmaxa)
275 
276  ! p10 non-negligible
277  if (abs(q2)/(abs(q2)+abs(mm02+mm12)).gt.dprec_coli) then
278 
279  f1 = q2-mm12+mm02
280 
281 ! do n0=1,rmax/2
282 ! k=rmax-2*n0
283  do r=2,2*rmax
284  do n0=max(1,r-rmax),r/2
285  k=r-2*n0
286  do n1=max(r0-2*n0,0),k
287  buv(n0,n1) = ((-1)**n1*auv(n0-1) + 2*mm02*buv(n0-1,n1) &
288  + f1*buv(n0-1,n1+1))/(2*n0+n1+1)/2
289  b(n0,n1) = ((-1)**n1*a(n0-1) + 2*mm02*b(n0-1,n1) &
290  + f1*b(n0-1,n1+1) + 4*buv(n0,n1))/(2*n0+n1+1)/2
291 ! write(*,*) 'CalcBred B(n0,n1)',n0,n1,B(n0,n1)
292  end do
293  end do
294  end do
295 
296  ! p10 negligible
297  else
298 
299 ! do n0=1,rmax/2
300 ! k=rmax-2*n0
301  do r=2,2*rmax
302  do n0=max(1,r-rmax),r/2
303  k=r-2*n0
304  do n1=max(r0-2*n0,0),k
305 
306  buv(n0,n1) = ((-1)**n1*auv(n0-1) + mm02*buv(n0-1,n1))/(n0+n1+1)/2
307  b(n0,n1) = ((-1)**n1*a(n0-1) + mm02*b(n0-1,n1) &
308  + 2*buv(n0,n1))/(n0+n1+1)/2
309  end do
310  end do
311  end do
312 
313  end if
314 
315  ! p10=m02=m12=0
316  else
317 
318  buv(1:,:) = 0d0
319  b(1:,:) = 0d0
320 
321  end if
322  end if
323 
324  berr(0) = abs(b(0,0))
325  do r=1,rmax
326  berr(r) = max(berr(r-1),abs(b(0,r)))
327  end do
328  berr = berr*acc_def_b
329 
330  end subroutine calcbred
331 
332 
333 
334 
335  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
336  ! subroutine CalcDB(DB,DBuv,p10,m02,m12,rmax,id,DBerr)
337  !
338  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
339 
340  subroutine calcdb(DB,DBuv,p10,m02,m12,rmax,id,DBerr)
341 
342  integer, intent(in) :: rmax, id
343  double complex, intent(in) :: p10,m02,m12
344  double complex, intent(inout) :: DB(0:rmax/2,0:rmax), DBuv(0:rmax/2,0:rmax)
345  double precision, optional, intent(out) :: DBerr(0:rmax)
346  double precision :: Berr(0:rmax)
347  double complex :: B(0:rmax,0:rmax), Buv(0:rmax,0:rmax)
348 ! double complex, allocatable :: A(:), Auv(:)
349  double complex :: A(0:rmax-1), Auv(0:rmax-1)
350  double complex :: q2, mm02, mm12, f1
351  double complex :: DBn_coli,elimminf2_coli
352  integer :: rmaxA,n0,n1,r,sgn,k,rid
353  logical :: finarg
354  double complex, parameter :: cd0 = dcmplx(0d0,0d0)
355 
356 ! write(*,*) 'CalcDB in',p10,m02,m12,rmax,id
357 
358  ! r=0
359  dbuv(0,0) = 0d0
360  db(0,0) = dbn_coli(0,p10,m02,m12)
361 
362  call calcbred(b,buv,p10,m02,m12,rmax,id,berr)
363 
364  ! calculate DB(0,n1)
365  do n1=1,rmax
366  sgn = -sgn
367  dbuv(0,n1) = 0d0
368  db(0,n1) = dbn_coli(n1,p10,m02,m12)
369  end do
370 
371 
372  ! calculate DB(n0,n1) for n0=/=0
373  if (rmax.ge.2) then
374 
375  call elminf2iv_coli(p10,m02,m12,q2,mm02,mm12,finarg)
376 
377 ! q2 = elimminf2_coli(p10)
378 ! mm12 = elimminf2_coli(m12)
379 ! mm02 = elimminf2_coli(m02)
380 ! finarg=q2.ne.cd0.or.mm12.ne.cd0.or.mm02.ne.0d0
381 
382  ! at least one of the invariants p10,m02,m12 different from zero
383  if (q2.ne.0d0.or.mm02.ne.0d0.or.mm12.ne.0d0) then
384 
385  ! calculation of A-functions
386  rmaxa = max(rmax-1,0)
387 ! allocate(A(0:rmaxA))
388 ! allocate(Auv(0:rmaxA))
389  call calca(a,auv,mm12,2*rmaxa)
390 
391  f1 = q2-mm12+mm02
392 
393  do r=2,rmax
394  do n0=max(1,r-rmax),r/2
395  k=r-2*n0
396  do n1=0,k
397  dbuv(n0,n1) = (buv(n0-1,n1+1) + 2*mm02*dbuv(n0-1,n1) &
398  + f1*dbuv(n0-1,n1+1))/(2*n0+n1+1)/2
399 
400  db(n0,n1) = b(n0-1,n1+1) + 4*dbuv(n0,n1)
401  if(mm02.ne.0d0) then
402  db(n0,n1) = db(n0,n1) + 2*mm02*db(n0-1,n1)
403  end if
404  if(f1.ne.0d0) then
405  db(n0,n1) = db(n0,n1) + f1*db(n0-1,n1+1)
406  end if
407 
408  db(n0,n1) = db(n0,n1) /(2*(2*n0+n1+1))
409 
410  end do
411  end do
412  end do
413 
414  ! p10=m02=m12=0
415  else
416 
417  dbuv(1:,:) = 0d0
418 
419  do r=2,rmax
420  do n0=max(1,r-rmax),r/2
421  k=r-2*n0
422  do n1=0,k
423  if (.not.ir_rational_terms_cll) then
424  dbuv(n0,n1) = buv(n0-1,n1+1)/(2*n0+n1+1)/2
425  end if
426  db(n0,n1) = (b(n0-1,n1+1) + 4*dbuv(n0,n1))/(2*n0+n1+1)/2
427  end do
428  end do
429  end do
430 
431  end if
432  end if
433 
434  dberr(0) = abs(db(0,0))
435  do r=1,rmax
436  dberr(r) = max(dberr(r-1),abs(db(0,r)))
437  end do
438  dberr = dberr*acc_def_b
439 
440  end subroutine calcdb
441 
442 
443 end module reductionab
reductionab::calcb
subroutine calcb(B, Buv, p10, m02, m12, rmax, id, Berr)
Definition: reductionAB.F90:92
reductionab::calcdb
subroutine calcdb(DB, DBuv, p10, m02, m12, rmax, id, DBerr)
Definition: reductionAB.F90:341
coli_aux2::dprec_coli
double precision dprec_coli
Definition: coli_aux2.F90:33
coli_aux2
Definition: coli_aux2.F90:23
collier_global
Definition: collier_global.F90:23
collier_global::ir_rational_terms_cll
logical ir_rational_terms_cll
Definition: collier_global.F90:39
reductionab::calcbred
subroutine calcbred(B, Buv, p10, m02, m12, rmax, id, Berr)
Definition: reductionAB.F90:227
reductionab::calca
subroutine calca(A, Auv, m02, rmax, Aerr)
Definition: reductionAB.F90:53
reductionab
Definition: reductionAB.F90:28
coli_aux2::acc_def_b
double precision acc_def_b
Definition: coli_aux2.F90:35