JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
Functions/Subroutines
reductiontn Module Reference

Functions/Subroutines

subroutine calctn (TN, TNuv, MomInv, masses2, N, rmax, id, TNerr, TNerr2)
 
recursive subroutine calctnint (TN, TNuv, MomInv, masses2, N, rmax, id, TNerr, TNerr2)
 
recursive subroutine calctnred (TN, TNuv, MomInv, masses2, N, rmax, id, TNerr, TNerr2)
 
double complex function, dimension(0:n-2) submasses (N, k, masses)
 
double complex function, dimension(0:3, n-2) submomvec (N, k, MomVec)
 
double complex function, dimension(binomtable(2, n-1)) submominv (N, k, MomInv)
 

Function/Subroutine Documentation

◆ calctn()

subroutine reductiontn::calctn ( double complex, dimension(ncoefs(rmax,n)), intent(inout)  TN,
double complex, dimension(ncoefs(rmax,n)), intent(inout)  TNuv,
double complex, dimension(binomtable(2,n)), intent(in)  MomInv,
double complex, dimension(0:n-1), intent(in)  masses2,
integer, intent(in)  N,
integer, intent(in)  rmax,
integer, intent(in)  id,
double precision, dimension(0:rmax), intent(out)  TNerr,
double precision, dimension(0:rmax), intent(out)  TNerr2 
)

Definition at line 45 of file reductionTN.F90.

45 
46  integer, intent(in) :: N,rmax,id
47  double complex, intent(in) :: MomInv(BinomTable(2,N)), masses2(0:N-1)
48  double complex, intent(inout) :: TN(NCoefs(rmax,N))
49  double complex, intent(inout) :: TNuv(NCoefs(rmax,N))
50  double precision, intent(out) :: TNerr(0:rmax),TNerr2(0:rmax)
51  double complex, allocatable :: TNaux(:,:,:),TNuvaux(:,:,:)
52  double complex :: B(0:rmax,0:rmax), Buv(0:rmax,0:rmax)
53  double complex :: C(0:rmax,0:rmax,0:rmax), Cuv(0:rmax,0:rmax,0:rmax)
54  double complex :: D(0:rmax,0:rmax,0:rmax,0:rmax), Duv(0:rmax,0:rmax,0:rmax,0:rmax)
55  double complex :: E(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
56  double complex :: Euv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
57  double complex :: x(BinomTable(2,N)+N)
58  double complex, allocatable :: fct(:)
59  integer :: r,n0,n1,n2,n3,n4,i,rank,bino,bino_m2,cnt
60  logical :: nocalc,wrica
61 
62 
63  select case (n)
64 
65  case(1)
66  call calca(tn,tnuv,masses2(0),rmax,tnerr)
67  tnerr2=tnerr
68  return
69 
70  case(2)
71  call calcb(b,buv,mominv(1),masses2(0),masses2(1),rmax,id,tnerr)
72 
73  cnt = 1
74  do r=0,rmax
75  do n0=r/2,0,-1
76  n1 = r-2*n0
77  tn(cnt) = b(n0,n1)
78  cnt = cnt+1
79  end do
80  end do
81  cnt = 1
82  do r=0,rmax
83  do n0=r/2,0,-1
84  n1 = r-2*n0
85  tnuv(cnt) = buv(n0,n1)
86  cnt = cnt+1
87  end do
88  end do
89  tnerr2=tnerr
90  return
91 
92  case(3)
93  call calcc(c,cuv,mominv(1),mominv(2),mominv(3), &
94  masses2(0),masses2(1),masses2(2),rmax,id,tnerr,tnerr2)
95 
96  cnt = 1
97  do r=0,rmax
98  do n0=r/2,0,-1
99  do n1=r-2*n0,0,-1
100  n2 = r-2*n0-n1
101  tn(cnt) = c(n0,n1,n2)
102  cnt = cnt+1
103  end do
104  end do
105  end do
106  tnuv = 0d0
107  do r=2,rmax
108  cnt = ncoefs(r-1,3)+1
109  do n0=r/2,1,-1
110  do n1=r-2*n0,0,-1
111  n2 = r-2*n0-n1
112  tnuv(cnt) = cuv(n0,n1,n2)
113  cnt = cnt+1
114  end do
115  end do
116  end do
117 
118 ! write(*,*) 'CalcTN 3 TNerr ',TNerr
119 ! write(*,*) 'CalcTN 3 TNerr2',TNerr2
120 
121  return
122 
123 
124 
125  case(4)
126  call calcd(d,duv,mominv(1),mominv(2),mominv(3),mominv(4),mominv(5),mominv(6), &
127  masses2(0),masses2(1),masses2(2),masses2(3),rmax,id,tnerr,tnerr2)
128 
129  cnt = 1
130  do r=0,rmax
131  do n0=r/2,0,-1
132  do n1=r-2*n0,0,-1
133  do n2=r-2*n0-n1,0,-1
134  n3 = r-2*n0-n1-n2
135  tn(cnt) = d(n0,n1,n2,n3)
136  cnt = cnt+1
137  end do
138  end do
139  end do
140  end do
141  tnuv = 0d0
142  do r=4,rmax
143  cnt = ncoefs(r-1,4)+1
144  do n0=r/2,2,-1
145  do n1=r-2*n0,0,-1
146  do n2=r-2*n0-n1,0,-1
147  n3 = r-2*n0-n1-n2
148  tnuv(cnt) = duv(n0,n1,n2,n3)
149  cnt = cnt+1
150  end do
151  end do
152  end do
153  end do
154 
155 ! write(*,*) 'CalcTN 4 TNerr ',TNerr
156 ! write(*,*) 'CalcTN 4 TNerr2',TNerr2
157 
158  return
159 
160  case(5)
161  call calce(e,euv,mominv(1),mominv(2),mominv(3),mominv(4),mominv(5), &
162  mominv(6),mominv(7),mominv(8),mominv(9),mominv(10), &
163  masses2(0),masses2(1),masses2(2),masses2(3),masses2(4) &
164  ,rmax,id,tnerr,tnerr2)
165 
166  cnt = 1
167  do r=0,rmax
168  do n0=r/2,0,-1
169  do n1=r-2*n0,0,-1
170  do n2=r-2*n0-n1,0,-1
171  do n3=r-2*n0-n1-n2,0,-1
172  n4 = r-2*n0-n1-n2-n3
173  tn(cnt) = e(n0,n1,n2,n3,n4)
174  cnt = cnt+1
175  end do
176  end do
177  end do
178  end do
179  end do
180  tnuv = 0d0
181  do r=6,rmax
182  cnt = ncoefs(r-1,5)+1
183  do n0=r/2,3,-1
184  do n1=r-2*n0,0,-1
185  do n2=r-2*n0-n1,0,-1
186  do n3=r-2*n0-n1-n2,0,-1
187  n4 = r-2*n0-n1-n2-n3
188  tnuv(cnt) = euv(n0,n1,n2,n3,n4)
189  cnt = cnt+1
190  end do
191  end do
192  end do
193  end do
194  end do
195  return
196 
197  case default
198  allocate(tnaux(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax))
199  allocate(tnuvaux(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax))
200  call calctnint(tnaux,tnuvaux,mominv,masses2,n,rmax,id,tnerr,tnerr2)
201 
202  cnt = 1
203  do r=0,rmax
204  do n0=r/2,0,-1
205  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
206  tn(cnt) = tnaux(i,n0,r)
207  cnt = cnt+1
208  end do
209  end do
210  end do
211  tnuv = 0d0
212  do r=2*n-4,rmax
213  cnt = ncoefs(r-1,n)+1
214  do n0=r/2,n-2,-1
215  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
216  tnuv(cnt) = tnuvaux(i,n0,r)
217  cnt = cnt+1
218  end do
219  end do
220  end do
221  return
222 
223  end select
224 
225 

◆ calctnint()

recursive subroutine reductiontn::calctnint ( double complex, dimension(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax), intent(inout)  TN,
double complex, dimension(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax), intent(inout)  TNuv,
double complex, dimension(binomtable(2,n)), intent(in)  MomInv,
double complex, dimension(0:n-1), intent(in)  masses2,
integer, intent(in)  N,
integer, intent(in)  rmax,
integer, intent(in)  id,
double precision, dimension(0:rmax), intent(out)  TNerr,
double precision, dimension(0:rmax), intent(out)  TNerr2 
)

Definition at line 238 of file reductionTN.F90.

238 
239  integer, intent(in) :: N,rmax,id
240  double complex, intent(in) :: MomInv(BinomTable(2,N)), masses2(0:N-1)
241  double complex, intent(inout) :: TN(BinomTable(rmax,N+rmax-2),0:rmax/2,0:rmax)
242  double complex, intent(inout) :: TNuv(BinomTable(rmax,N+rmax-2),0:rmax/2,0:rmax)
243  double precision, intent(out) :: TNerr(0:rmax),TNerr2(0:rmax)
244  double complex, allocatable :: TNaux(:,:,:), TNuvaux(:,:,:)
245  double complex :: E(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
246  double complex :: Euv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
247  double precision, allocatable :: TNerraux(:),TNerr2aux(:)
248  double complex :: x(BinomTable(2,N)+N)
249  double complex, allocatable :: fct(:)
250  integer :: r,n0,n1,n2,n3,n4,i,rank,bino,bino_m2,cnt
251  logical :: nocalc,wrica,noten
252 
253 
254 
255  if (n.eq.5) then
256  call calce(e,euv,mominv(1),mominv(2),mominv(3),mominv(4),mominv(5), &
257  mominv(6),mominv(7),mominv(8),mominv(9),mominv(10), &
258  masses2(0),masses2(1),masses2(2),masses2(3),masses2(4),rmax,id,tnerr,tnerr2)
259 
260  do r=0,rmax
261  do n0=r/2,0,-1
262  cnt = 1
263  do n1=r-2*n0,0,-1
264  do n2=r-2*n0-n1,0,-1
265  do n3=r-2*n0-n1-n2,0,-1
266  n4 = r-2*n0-n1-n2-n3
267  tn(cnt,n0,r) = e(n0,n1,n2,n3,n4)
268  cnt = cnt+1
269  end do
270  end do
271  end do
272  end do
273  end do
274  do r=6,rmax
275  do n0=r/2,3,-1
276  cnt = 1
277  do n1=r-2*n0,0,-1
278  do n2=r-2*n0-n1,0,-1
279  do n3=r-2*n0-n1-n2,0,-1
280  n4 = r-2*n0-n1-n2-n3
281  tnuv(cnt,n0,r) = euv(n0,n1,n2,n3,n4)
282  cnt = cnt+1
283  end do
284  end do
285  end do
286  end do
287  end do
288  return
289  end if
290 
291  if ((use_cache_system).and.(tencache.gt.n)) then
292  if ((ncache.gt.0).and.(ncache.le.ncache_max)) then
293  bino = binomtable(2,n)
294  x(1:bino) = mominv
295  x(bino+1:bino+n) = masses2
296  rank = rmax
297 
298  if(rmax.ge.2*n-4) then
299  allocate(fct(ncoefs(rmax,n)+ncoefs(rmax-2*n+4,n)+2*(rmax+1)))
300  call readcache(fct,ncoefs(rmax,n)+ncoefs(rmax-2*n+4,n)+2*(rmax+1),x,bino+n,1,id,n,rank,nocalc,wrica)
301  else
302  allocate(fct(ncoefs(rmax,n)+2*(rmax+1)))
303  call readcache(fct,ncoefs(rmax,n)+2*(rmax+1),x,bino+n,1,id,n,rank,nocalc,wrica)
304  end if
305 
306  if(nocalc)then
307  cnt = 0
308  do r=0,rmax
309  do n0=r/2,0,-1
310  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
311  cnt = cnt+1
312  tn(i,n0,r) = fct(cnt)
313  end do
314  end do
315  do n0=r/2,n-2,-1
316  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
317  cnt = cnt+1
318  tnuv(i,n0,r) = fct(cnt)
319  end do
320  end do
321  cnt = cnt+1
322  tnerr(r) = real(fct(cnt))
323  cnt = cnt+1
324  tnerr2(r) = real(fct(cnt))
325  end do
326  return
327  endif
328 
329 
330  if(rank.eq.rmax) then
331 
332  call calctnred(tn,tnuv,mominv,masses2,n,rank,id,tnerr,tnerr2)
333 
334  if (wrica) then
335  cnt = 0
336  do r=0,rank
337  do n0=r/2,0,-1
338  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
339 
340  cnt = cnt+1
341  fct(cnt) = tn(i,n0,r)
342 
343  end do
344  end do
345  do n0=r/2,n-2,-1
346  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
347 
348  cnt = cnt+1
349  fct(cnt) = tnuv(i,n0,r)
350 
351  end do
352  end do
353  cnt = cnt+1
354  fct(cnt) = tnerr(r)
355  cnt = cnt+1
356  fct(cnt) = tnerr2(r)
357  end do
358 
359  if(rank.ge.2*n-4) then
360  call writecache(fct,ncoefs(rank,n)+ncoefs(rank-2*n+4,n)+2*(rank+1),id,n,rank)
361  else
362  call writecache(fct,ncoefs(rank,n)+2*(rank+1),id,n,rank)
363  end if
364 
365  end if
366 
367  return
368 
369 
370  else
371  allocate(tnaux(binomtable(rank,n+rank-2),0:rank/2,0:rank))
372  allocate(tnuvaux(binomtable(rank,n+rank-2),0:rank/2,0:rank))
373  allocate(tnerraux(0:rank))
374  allocate(tnerr2aux(0:rank))
375 
376  call calctnred(tnaux,tnuvaux,mominv,masses2,n,rank,id,tnerraux,tnerr2aux)
377 
378  if (wrica) then
379  cnt = 0
380  deallocate(fct)
381  if(rank.ge.2*n-4) then
382  allocate(fct(ncoefs(rank,n)+ncoefs(rank-2*n+4,n)+2*(rank+1)))
383  else
384  allocate(fct(ncoefs(rank,n)+2*(rank+1)))
385  end if
386  do r=0,rank
387  do n0=r/2,0,-1
388  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
389 
390  cnt = cnt+1
391  fct(cnt) = tnaux(i,n0,r)
392 
393  end do
394  end do
395  do n0=r/2,n-2,-1
396  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
397 
398  cnt = cnt+1
399  fct(cnt) = tnuvaux(i,n0,r)
400 
401  end do
402  end do
403  cnt = cnt+1
404  fct(cnt) = tnerraux(r)
405  cnt = cnt+1
406  fct(cnt) = tnerr2aux(r)
407  end do
408 
409  if(rank.ge.2*n-4) then
410  call writecache(fct,ncoefs(rank,n)+ncoefs(rank-2*n+4,n)+2*(rank+1),id,n,rank)
411  else
412  call writecache(fct,ncoefs(rank,n)+2*(rank+1),id,n,rank)
413  end if
414 
415  end if
416 
417  tn = tnaux(1:binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax)
418  tnuv = tnuvaux(1:binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax)
419  tnerr = tnerraux(0:rmax)
420  tnerr2 = tnerr2aux(0:rmax)
421  return
422 
423  end if
424  end if
425  end if
426 
427  call calctnred(tn,tnuv,mominv,masses2,n,rmax,id,tnerr,tnerr2)
428 
429 

◆ calctnred()

recursive subroutine reductiontn::calctnred ( double complex, dimension(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax), intent(inout)  TN,
double complex, dimension(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax), intent(inout)  TNuv,
double complex, dimension(binomtable(2,n)), intent(in)  MomInv,
double complex, dimension(0:n-1), intent(in)  masses2,
integer, intent(in)  N,
integer, intent(in)  rmax,
integer, intent(in)  id,
double precision, dimension(0:rmax), intent(out)  TNerr,
double precision, dimension(0:rmax), intent(out)  TNerr2 
)

Definition at line 442 of file reductionTN.F90.

442 
443  integer, intent(in) :: N,rmax,id
444  double complex, intent(in) :: MomInv(BinomTable(2,N)), masses2(0:N-1)
445  double precision, intent(out) :: TNerr(0:rmax),TNerr2(0:rmax)
446  double complex :: q10,q21,q32,q43,q54,q20,q31,q42,q53,q50,q30,q41,q52,q40,q51
447  double complex :: mm02,mm12,mm22,mm32,mm42,mm52
448  double complex :: mx(0:5,0:5), mxinv(0:5,0:5),mx0k(5,5),mx0kinv(5,5),ff(N-1)
449  double complex :: det,newdet,TNaux,mxinvs,mx0kinvs(5)
450  double complex :: zmx0kinv(5,5),Z(5,5),zmx0kinvs(5)
451  double precision :: maxZ,maxzmx0kinv(5),maxzmx0kinvs
452  double complex, intent(inout) :: TN(BinomTable(rmax,N+rmax-2),0:rmax/2,0:rmax)
453  double complex, intent(inout) :: TNuv(BinomTable(rmax,N+rmax-2),0:rmax/2,0:rmax)
454  double complex, allocatable :: TNm1_0(:,:,:), TNm1UV_0(:,:,:), TNm1_i(:,:,:,:)
455  double complex, allocatable :: TNm1_0aux(:,:,:), TNm1UV_0aux(:,:,:), TNm1UV_i(:,:,:,:)
456  double complex :: S(5), elimminf2_coli,chdet,gram(1:5,1:5),gramdet
457  double precision :: TNm1err(0:5,0:rmax),TNm1err2(0:5,0:rmax)
458  integer :: r,n0,n1,n2,n3,n4,n5,n6,k,i,j,m,nid(0:5),r0,bin,kbest,Nm1,inds(N-1)
459  integer :: bino,bino_0,bino_i,bino_m1,bino_m2,cnt,iaux,rmax_m1,shift,cind,rBCD
460  integer :: mia1,mia2,mia3
461  logical :: errorwriteflag
462  character(len=*),parameter :: fmt10 = "(A18,'(',d25.18,' , ',d25.18,' )')"
463 
464 #ifdef TNtest
465  character(len=*),parameter :: fmt999 = "(' gm(',i1,',',i1,') = ',d23.16,' , ',d23.16)"
466 #endif
467 
468 ! write(*,*) 'CalcTNred in N',N,rmax
469 ! write(*,*) 'CalcTNred in s',MomInv
470 ! write(*,*) 'CalcTNred in m',masses2
471 
472  ! allocation of TN functions
473  nm1 = n-1
474  rmax_m1 = max(rmax-1,0)
475  bino_0 = binomtable(rmax_m1,n+rmax_m1-2)
476  bino_i = binomtable(rmax_m1,nm1+rmax_m1-2)
477  allocate(tnm1_0(bino_0,0:rmax_m1/2,0:rmax_m1))
478  allocate(tnm1_0aux(bino_i,0:rmax_m1/2,0:rmax_m1))
479  allocate(tnm1_i(bino_i,0:rmax_m1/2,0:rmax_m1,5))
480  allocate(tnm1uv_0(bino_0,0:rmax_m1/2,0:rmax_m1))
481  allocate(tnm1uv_0aux(bino_i,0:rmax_m1/2,0:rmax_m1))
482  allocate(tnm1uv_i(bino_i,0:rmax_m1/2,0:rmax_m1,5))
483 
484  ! determine binaries for TN-coefficients
485  k=0
486  bin = 1
487  do while (k.le.5)
488  if (mod(id/bin,2).eq.0) then
489  nid(k) = id+bin
490  k = k+1
491  end if
492  bin = 2*bin
493  end do
494 
495  call calctnint(tnm1_0aux,tnm1uv_0aux,submominv(n,0,mominv),submasses(n,0,masses2), &
496  nm1,rmax_m1,nid(0),tnerr=tnm1err(0,0:rmax_m1), &
497  tnerr2=tnm1err2(0,0:rmax_m1))
498  do k=1,5
499  call calctnint(tnm1_i(:,:,:,k),tnm1uv_i(:,:,:,k),submominv(n,k,mominv), &
500  submasses(n,k,masses2),nm1,rmax_m1,nid(k), &
501  tnerr=tnm1err(k,0:rmax_m1),tnerr2=tnm1err2(k,0:rmax_m1))
502  end do
503 
504  ! shift of integration momentum in TN\{0}
505  do n0=0,rmax_m1/2
506  tnm1_0(1,n0,2*n0) = tnm1_0aux(1,n0,2*n0)
507  do r=1,rmax_m1-2*n0
508  mia1 = binomtable(r-1,n+r-3)+1
509  mia2 = binomtable(r,n+r-3)
510  mia3 = binomtable(r,n+r-2)
511  tnm1_0(mia1:mia3,n0,r+2*n0) = tnm1_0aux(1:mia2,n0,r+2*n0)
512  end do
513  end do
514  do r=1,rmax_m1
515  do n0=0,(r-1)/2
516  do i=binomtable(r-2*n0-1,n+r-2*n0-3),1,-1
517  tnm1_0(i,n0,r) = -tnm1_0(i,n0,r-1)
518  do j=2,n-1
519  tnm1_0(i,n0,r) = tnm1_0(i,n0,r) &
520  - tnm1_0(addtocind(j,i,r-2*n0-1,n-1),n0,r)
521  end do
522  end do
523  end do
524  end do
525  do n0=n-3,rmax_m1/2
526  tnm1uv_0(1,n0,2*n0) = tnm1uv_0aux(1,n0,2*n0)
527  do r=1,rmax_m1-2*n0
528  mia1 = binomtable(r-1,n+r-3)+1
529  mia2 = binomtable(r,n+r-3)
530  mia3 = binomtable(r,n+r-2)
531  tnm1uv_0(mia1:mia3,n0,r+2*n0) = tnm1uv_0aux(1:mia2,n0,r+2*n0)
532  end do
533  end do
534  do r=2*n-6,rmax_m1
535  do n0=n-3,(r-1)/2
536  do i=binomtable(r-2*n0-1,n+r-2*n0-3),1,-1
537  tnm1uv_0(i,n0,r) = -tnm1uv_0(i,n0,r-1)
538  do j=2,n-1
539  tnm1uv_0(i,n0,r) = tnm1uv_0(i,n0,r) &
540  - tnm1uv_0(addtocind(j,i,r-2*n0-1,n-1),n0,r)
541  end do
542  end do
543  end do
544  end do
545 
546 
547  ! determine inverse modified Caley matrix
548  mm02 = elimminf2_coli(masses2(0))
549  mm12 = elimminf2_coli(masses2(1))
550  mm22 = elimminf2_coli(masses2(2))
551  mm32 = elimminf2_coli(masses2(3))
552  mm42 = elimminf2_coli(masses2(4))
553  mm52 = elimminf2_coli(masses2(5))
554  q10 = elimminf2_coli(mominv(1))
555  q21 = elimminf2_coli(mominv(2))
556  q32 = elimminf2_coli(mominv(3))
557  q43 = elimminf2_coli(mominv(4))
558  q54 = elimminf2_coli(mominv(5))
559  if (n.le.9) then
560  q50 = elimminf2_coli(mominv((n-6)*n+6))
561  else
562  q50 = elimminf2_coli(mominv(4*n+1))
563  end if
564  q20 = elimminf2_coli(mominv(n+1))
565  q31 = elimminf2_coli(mominv(n+2))
566  q42 = elimminf2_coli(mominv(n+3))
567  q53 = elimminf2_coli(mominv(n+4))
568  if (n.le.7) then
569  q40 = elimminf2_coli(mominv((n-5)*n+5))
570  q51 = elimminf2_coli(mominv((n-5)*n+6))
571  else
572  q40 = elimminf2_coli(mominv(3*n+1))
573  q51 = elimminf2_coli(mominv(3*n+2))
574  end if
575  q30 = elimminf2_coli(mominv(2*n+1))
576  q41 = elimminf2_coli(mominv(2*n+2))
577  q52 = elimminf2_coli(mominv(2*n+3))
578  ! write(*,*) 'q50', q50
579  ! write(*,*) 'q20', q20
580  ! write(*,*) 'q31', q31
581  ! write(*,*) 'q42', q42
582  ! write(*,*) 'q53', q53
583  ! write(*,*) 'q40', q40
584  ! write(*,*) 'q51', q51
585  ! write(*,*) 'q30', q30
586  ! write(*,*) 'q41', q41
587  ! write(*,*) 'q52', q52
588 
589  cnt = 1
590  do i=1,n-1
591  ff(i) = elimminf2_coli(mominv(cnt)) + mm02 &
592  - elimminf2_coli(masses2(i))
593  cnt = cnt+n-i
594  end do
595 
596  mx(0,0) = 2d0*mm02
597  mx(1,0) = q10 - mm12 + mm02
598  mx(2,0) = q20 - mm22 + mm02
599  mx(3,0) = q30 - mm32 + mm02
600  mx(4,0) = q40 - mm42 + mm02
601  mx(5,0) = q50 - mm52 + mm02
602  mx(0,1) = mx(1,0)
603  mx(1,1) = 2d0*q10
604  mx(2,1) = q10+q20-q21
605  mx(3,1) = q10+q30-q31
606  mx(4,1) = q10+q40-q41
607  mx(5,1) = q10+q50-q51
608  mx(0,2) = mx(2,0)
609  mx(1,2) = mx(2,1)
610  mx(2,2) = 2d0*q20
611  mx(3,2) = q20+q30-q32
612  mx(4,2) = q20+q40-q42
613  mx(5,2) = q20+q50-q52
614  mx(0,3) = mx(3,0)
615  mx(1,3) = mx(3,1)
616  mx(2,3) = mx(3,2)
617  mx(3,3) = 2d0*q30
618  mx(4,3) = q30+q40-q43
619  mx(5,3) = q30+q50-q53
620  mx(0,4) = mx(4,0)
621  mx(1,4) = mx(4,1)
622  mx(2,4) = mx(4,2)
623  mx(3,4) = mx(4,3)
624  mx(4,4) = 2d0*q40
625  mx(5,4) = q40+q50-q54
626  mx(0,5) = mx(5,0)
627  mx(1,5) = mx(5,1)
628  mx(2,5) = mx(5,2)
629  mx(3,5) = mx(5,3)
630  mx(4,5) = mx(5,4)
631  mx(5,5) = 2d0*q50
632 
633  call chinv(6,mx,mxinv)
634 
635  ! determine X_(0,5)
636  do j=1,5
637  do i=1,5
638  mx0k(i,j) = mx(i,j-1)
639  end do
640  end do
641 
642  det = chdet(5,mx0k)
643  kbest = 5
644 
645 #ifdef TNtest
646  write(*,*) 'det',5,det
647 #endif
648 
649  do j=5,2,-1
650  do i=1,5
651  mx0k(i,j) = mx(i,j)
652  end do
653 
654  newdet = chdet(5,mx0k)
655  if (abs(newdet).gt.abs(det)) then
656  kbest = j-1
657 #ifdef TNtest
658  write(*,*) 'det',j-1,newdet
659 #endif
660  det = newdet
661  end if
662 
663  end do
664 
665 #ifdef TNtest
666 ! kbest=5
667  write(*,*) 'kbest',kbest
668 #endif
669 
670 
671 
672 
673  do i=1,5
674  mx0k(i,1) = mx(i,1)
675  mx0k(i,kbest) = mx(i,0)
676  end do
677 
678  call chinv(5,mx0k,mx0kinv)
679  do i=1,5
680  mx0kinv(kbest,i) = 0d0
681  end do
682 
683  mxinvs = sum(mxinv(0:5,0))
684  do i=1,5
685  mx0kinvs(i) = sum(mx0kinv(i,1:5))
686  end do
687 
688  ! for alternative error estimate
689  z(1:5,1:5) = mx(1:5,1:5)
690 
691  zmx0kinv = matmul(z,mx0kinv)
692 
693  do i=1,5
694  maxzmx0kinv(i) = maxval(abs(zmx0kinv(1:5,i)))
695  zmx0kinvs(i) = sum(zmx0kinv(i,1:5))
696  end do
697 
698  maxzmx0kinvs = maxval(abs(zmx0kinvs(1:5)))
699 
700  maxz = maxval(abs(z))
701 
702  ! calculation of UV-divergent parts
703  ! TN UV-finite for n0<N-2
704  ! TNuv(:,0:min(rmax/2,N-3),:) = 0d0
705 
706  ! PV reduction (5.10)
707  do r=2*n-4,rmax
708  do n0=n-2,r/2
709  do cind=1,binomtable(r-2*n0,n+r-2*n0-2)
710 
711  tnuv(cind,n0,r) = tnm1uv_0(cind,n0-1,r-2) + 2*mm02*tnuv(cind,n0-1,r-2)
712  do i=1,nm1
713  tnuv(cind,n0,r) = tnuv(cind,n0,r) &
714  + ff(i)*tnuv(addtocind(i,cind,r-2*n0,n-1),n0-1,r-1)
715  end do
716  tnuv(cind,n0,r) = tnuv(cind,n0,r)/(2*(r+3-n))
717 
718  end do
719  end do
720  end do
721 
722  ! scalar coefficient
723  tn = 0d0
724  tn(1,0,0) = -mxinv(0,0)*tnm1_0(1,0,0)
725  do k=1,5
726  tn(1,0,0) = tn(1,0,0) &
727  + mxinv(k,0)*(tnm1_i(1,0,0,k)-tnm1_0(1,0,0))
728  end do
729 
730 ! TNerr(0) = max( maxval(abs(mxinv(0:4,0)))*TNm1err(0,0), &
731  tnerr(0) = max(abs(mxinvs)*tnm1err(0,0), &
732  abs(mxinv(1,0))*tnm1err(1,0) , &
733  abs(mxinv(2,0))*tnm1err(2,0) , &
734  abs(mxinv(3,0))*tnm1err(3,0) , &
735  abs(mxinv(4,0))*tnm1err(4,0) , &
736  abs(mxinv(5,0))*tnm1err(5,0) )
737 
738  tnerr2(0) = max(abs(mxinvs)*tnm1err2(0,0), &
739  abs(mxinv(1,0))*tnm1err2(1,0) , &
740  abs(mxinv(2,0))*tnm1err2(2,0) , &
741  abs(mxinv(3,0))*tnm1err2(3,0) , &
742  abs(mxinv(4,0))*tnm1err2(4,0) , &
743  abs(mxinv(5,0))*tnm1err2(5,0) )
744 
745  if (rmax.eq.0) return
746 
747 
748  ! formula (7.13) extended to arbitrary N
749  do r=0,rmax_m1
750  do n0=0,max((r-n+6)/2,0)
751  bino_m1 = binomtable(r-2*n0,n+r-2*n0-2)
752  do i=1,bino_m1
753  inds = calccindarr(nm1,r-2*n0,i)
754 
755  do m=1,5
756  s(m) = -tnm1_0(i,n0,r)
757  end do
758 
759 #ifdef TNtest
760  write(*,*) 'Tnred s(m)',s
761 #endif
762 
763  if (inds(1).eq.0) then
764  s(1) = s(1) + tnm1_i(dropcind(1,i,r-2*n0,nm1),n0,r,1)
765  end if
766  if (inds(2).eq.0) then
767  s(2) = s(2) + tnm1_i(dropcind(2,i,r-2*n0,nm1),n0,r,2)
768  end if
769  if (inds(3).eq.0) then
770  s(3) = s(3) + tnm1_i(dropcind(3,i,r-2*n0,nm1),n0,r,3)
771  end if
772  if (inds(4).eq.0) then
773  s(4) = s(4) + tnm1_i(dropcind(4,i,r-2*n0,nm1),n0,r,4)
774  end if
775  if (inds(5).eq.0) then
776  s(5) = s(5) + tnm1_i(dropcind(5,i,r-2*n0,nm1),n0,r,5)
777  end if
778 
779 #ifdef TNtest
780  write(*,*) 'Tnred s(m)',s
781 #endif
782 
783  do k=1,5
784  tnaux = mx0kinv(k,1)*s(1)+mx0kinv(k,2)*s(2) &
785  + mx0kinv(k,3)*s(3)+mx0kinv(k,4)*s(4)+mx0kinv(k,5)*s(5)
786  iaux = addtocind(k,i,r,n-1)
787  tn(iaux,n0,r+1) = tn(iaux,n0,r+1) + (inds(k)+1)*tnaux/(r+1)
788 #ifdef TNtest
789  write(*,*) 'Tnred comb',k,mx0kinv(k,1)*s(1),mx0kinv(k,2)*s(2) &
790  , mx0kinv(k,3)*s(3),mx0kinv(k,4)*s(4),mx0kinv(k,5)*s(5), tnaux
791 
792  if(abs(tnaux).ne.0d0) then
793  write(*,*) 'Tnred cancel', k, max(abs(mx0kinv(k,1)*s(1)),abs(mx0kinv(k,2)*s(2)) &
794  , abs(mx0kinv(k,3)*s(3)),abs(mx0kinv(k,4)*s(4)),abs(mx0kinv(k,5)*s(5))) &
795  /abs(tnaux)
796  endif
797 #endif
798 
799  end do
800 
801 
802 
803  end do
804  end do
805 
806  if (r.le.rmax-1) then
807 
808  tnerr(r+1) = max( maxval(abs(mx0kinvs(1:5)))*tnm1err(0,r), &
809  maxval(abs(mx0kinv(1:5,1)))*tnm1err(1,r) , &
810  maxval(abs(mx0kinv(1:5,2)))*tnm1err(2,r) , &
811  maxval(abs(mx0kinv(1:5,3)))*tnm1err(3,r) , &
812  maxval(abs(mx0kinv(1:5,4)))*tnm1err(4,r) , &
813  maxval(abs(mx0kinv(1:5,5)))*tnm1err(5,r) )
814 
815  tnerr2(r+1) = max( abs(maxzmx0kinvs)*tnm1err2(0,r), &
816  abs(maxzmx0kinv(1))*tnm1err2(1,r) , &
817  abs(maxzmx0kinv(2))*tnm1err2(2,r) , &
818  abs(maxzmx0kinv(3))*tnm1err2(3,r) , &
819  abs(maxzmx0kinv(4))*tnm1err2(4,r) , &
820  abs(maxzmx0kinv(5))*tnm1err2(5,r) )/maxz
821 
822 
823  if (mode_coli.lt.1) then
824  gram= mx(1:5,1:5)
825 
826  gramdet = chdet(5,gram)
827 
828 #ifdef TNtest
829  write(*,*)
830  do i = 1,5
831  do j = 1,5
832  write(*,fmt999) i,j,gram(i,j)
833  enddo
834  enddo
835  write(*,*) 'TNred relGramdet=',gramdet/det
836 #endif
837 
838  if (max(abs(tn(1,0,0)),maxval(abs(tn(1:6,0,1))))*abs(gramdet/det).gt. &
839  tnerr(r+1)) then
840 
841  tnerr(r+1)= max(tnerr(r+1), &
842  max(abs(tn(1,0,0)),maxval(abs(tn(1:6,0,1))))*abs(gramdet/det) )
843  tnerr2(r+1)= max(tnerr2(r+1), &
844  max(abs(tn(1,0,0)),maxval(abs(tn(1:6,0,1))))*abs(gramdet/det) )
845 
846  if (abs(gramdet/det).gt.reqacc_coli) then
847  call seterrflag_coli(-6)
848  call errout_coli('CalcTNred', &
849  'input momenta inconsistent! (not 4-dimensional)', &
850  errorwriteflag)
851  if (errorwriteflag) then
852  write(nerrout_coli,fmt10) ' CalcTNred: q10 = ',q10
853  write(nerrout_coli,fmt10) ' CalcTNred: q21 = ',q21
854  write(nerrout_coli,fmt10) ' CalcTNred: q32 = ',q32
855  write(nerrout_coli,fmt10) ' CalcTNred: q43 = ',q43
856  write(nerrout_coli,fmt10) ' CalcTNred: q54 = ',q54
857  write(nerrout_coli,fmt10) ' CalcTNred: q50 = ',q50
858  write(nerrout_coli,fmt10) ' CalcTNred: q20 = ',q10
859  write(nerrout_coli,fmt10) ' CalcTNred: q31 = ',q31
860  write(nerrout_coli,fmt10) ' CalcTNred: q42 = ',q42
861  write(nerrout_coli,fmt10) ' CalcTNred: q53 = ',q53
862  write(nerrout_coli,fmt10) ' CalcTNred: q40 = ',q40
863  write(nerrout_coli,fmt10) ' CalcTNred: q51 = ',q51
864  write(nerrout_coli,fmt10) ' CalcTNred: q30 = ',q30
865  write(nerrout_coli,fmt10) ' CalcTNred: q41 = ',q41
866  write(nerrout_coli,fmt10) ' CalcTNred: q52 = ',q52
867  write(nerrout_coli,fmt10) ' CalcTNred: mm02 = ',mm02
868  write(nerrout_coli,fmt10) ' CalcTNred: mm12 = ',mm12
869  write(nerrout_coli,fmt10) ' CalcTNred: mm22 = ',mm22
870  write(nerrout_coli,fmt10) ' CalcTNred: mm32 = ',mm32
871  write(nerrout_coli,fmt10) ' CalcTNred: mm42 = ',mm42
872  write(nerrout_coli,fmt10) ' CalcTNred: mm52 = ',mm52
873  write(nerrout_coli,fmt10) ' CalcTNred: gram = ',gramdet/det
874  end if
875  end if
876 
877  end if
878 
879  end if
880  end if
881 
882 
883 #ifdef TNtest
884  write(*,*) 'Tnred mn err',tnm1err(1:5,r)
885  write(*,*) 'Tnred coef err', (maxval(abs(mx0kinv(1:5,k))),k=1,5)
886  write(*,*) 'Tnred err',tnerr(r+1)
887 #endif
888 
889  end do
890 
891 

◆ submasses()

double complex function, dimension(0:n-2) reductiontn::submasses ( integer  N,
integer  k,
double complex, dimension(0:n-1)  masses 
)

Definition at line 904 of file reductionTN.F90.

904 
905  integer :: N,k,i
906  double complex :: masses(0:N-1), SubMasses(0:N-2)
907 
908  if ((k.gt.n-1).or.(k.lt.0)) then
909  write(nerrout_coli,*) 'SubMasses:'
910  write(nerrout_coli,*) 'inkonsistent argument k', k
911  write(nerrout_coli,*) 0, '<= k <=', n-1, 'required!'
912  stop
913  end if
914 
915  do i=0,n-2
916  if (i.lt.k) then
917  submasses(i) = masses(i)
918  else
919  submasses(i) = masses(i+1)
920  end if
921  end do
922 

◆ submominv()

double complex function, dimension(binomtable(2,n-1)) reductiontn::submominv ( integer, intent(in)  N,
integer, intent(in)  k,
double complex, dimension(binomtable(2,n))  MomInv 
)

Definition at line 977 of file reductionTN.F90.

977 
978  integer, intent(in) :: N,k
979  integer :: i,cnt,moms,limit
980  double complex :: MomInv(BinomTable(2,N)), SubMomInv(BinomTable(2,N-1))
981 
982  ! write(*,*) 'SubMomIv', N,k
983  if ((k.gt.n-1).or.(k.lt.0)) then
984  write(nerrout_coli,*) 'SubMomInv:'
985  write(nerrout_coli,*) 'inkonsistent argument k', k
986  write(nerrout_coli,*) 0, '<= k <=', n-1, 'required!'
987  stop
988  end if
989 
990  if (n.eq.2) return
991 
992  cnt = 1
993  do moms=1,(n-3)/2
994  do i=1,n
995  if (i.ne.k+1) then
996  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
997  submominv(cnt) = mominv(n*moms+i)
998  else
999  submominv(cnt) = mominv(n*(moms-1)+i)
1000  end if
1001 
1002  cnt = cnt+1
1003  end if
1004  end do
1005  end do
1006 
1007  if (mod(n,2).eq.1) then
1008  moms = (n-1)/2
1009  if(k.le.(n-1)/2) then
1010  limit = (n+1)/2
1011  else
1012  limit = (n-1)/2
1013  end if
1014  do i=1,limit
1015  if (i.ne.k+1) then
1016  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
1017  submominv(cnt) = mominv(n*moms+i-(n-1)/2)
1018  else
1019  submominv(cnt) = mominv(n*(moms-1)+i)
1020  end if
1021 
1022  cnt = cnt+1
1023  end if
1024  end do
1025 
1026  else
1027  do i=1,n
1028  if (i.ne.k+1) then
1029  moms = (n/2-1)
1030  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
1031  if (i.gt.n/2) then
1032  submominv(cnt) = mominv(n*moms-n/2+i)
1033  else
1034  submominv(cnt) = mominv(n*moms+i)
1035  endif
1036  else
1037  submominv(cnt) = mominv(n*(moms-1)+i)
1038  end if
1039 
1040  cnt = cnt+1
1041  end if
1042  end do
1043 
1044  end if
1045 
1046 

◆ submomvec()

double complex function, dimension(0:3,n-2) reductiontn::submomvec ( integer  N,
integer  k,
double complex, dimension(0:3,n-1)  MomVec 
)

Definition at line 935 of file reductionTN.F90.

935 
936  integer :: N,k,i
937  double complex :: MomVec(0:3,N-1), SubMomVec(0:3,N-2)
938 
939  if ((k.gt.n-1).or.(k.lt.0)) then
940  write(nerrout_coli,*) 'SubMomVec:'
941  write(nerrout_coli,*) 'inkonsistent argument k', k
942  write(nerrout_coli,*) 0, '<= k <=', n-1, 'required!'
943  stop
944  end if
945 
946  if (k.eq.0) then
947  do i=1,n-2
948  submomvec(0:3,i) = momvec(0:3,i+1)-momvec(0:3,1)
949  end do
950 
951  else
952  do i=1,n-2
953  if (i.lt.k) then
954  submomvec(0:3,i) = momvec(0:3,i)
955  else
956  submomvec(0:3,i) = momvec(0:3,i+1)
957  end if
958  end do
959 
960  end if
961 
submasses
double complex function, dimension(0:n-2) submasses(N, k, masses)
Definition: democache.f90:468
submominv
double complex function, dimension((n-1) *(n-2)/2) submominv(N, k, MomInv)
Definition: democache.f90:495