JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
TensorReduction.F90
Go to the documentation of this file.
1 !!
2 !! File TensorReduction.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 TensorReduction *
14 ! * by Lars Hofer *
15 ! ****************************
16 !
17 ! functions and subroutines:
18 ! init_tables, CalcTensorNr, CalcTNrPVco, CalcTNrPVco_nd, CalcTtilde, CalcTgtilde,
19 ! CalcTensorEr, EreduD0, EreduD1, EreduDr, CalcTensorFr, FreduE0, FreduEr
20 !
21 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
22 
23 
24 
25 module tensorreduction
26 
27  use collier_coefs
28  use buildtensors
29 
30  implicit none
31 
32 
33 contains
34 
35 
36 
37 
38 
39  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
40  ! subroutine CalcTensorFr(TF,TFuv,TFerr,MomVec,MomInv,masses2,rmax,CFuv)
41  !
42  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
43 
44  subroutine calctensorfr(TF,TFuv,TFerr,MomVec,MomInv,masses2,rmax)
45 
46  integer, intent(in) :: rmax
47  double complex, intent(in) :: MomVec(0:3,5), masses2(0:5), MomInv(15)
48  double complex, intent(out) :: TF(0:rmax,0:rmax,0:rmax,0:rmax)
49  double complex, intent(out) :: TFuv(0:rmax,0:rmax,0:rmax,0:rmax)
50  double precision, intent(out) :: TFerr(0:rmax)
51  double complex :: CFuv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
52  double complex, allocatable :: TFaux(:,:,:,:), CFuvaux(:,:,:,:,:,:)
53  double precision, allocatable :: TFerraux(:)
54  double complex :: x(41)
55  double complex, allocatable :: fct(:)
56  integer :: r,n0,n1,n2,n3,n4,n5,i,rank,bino,cnt
57  logical :: nocalc,wrica
58 
59 
60  if (tenred_cll.gt.6) then
61  if (erroutlev_cll.ge.1) then
62  write(nerrout_cll,*) 'inconsistent call of CalcTensorFr'
63  end if
64  stop
65  endif
66 
67  if (use_cache_system) then
68  if ((ncache.gt.0).and.(ncache.le.ncache_max)) then
69 ! if (use_cache(ncache).ge.6) then
70  x(1:15) = mominv
71  x(16:21) = masses2
72  do i=1,5
73  x(22+(i-1)*4:21+i*4) = momvec(0:3,i)
74  end do
75  rank = rmax
76 
77  if (rmax.ge.8) then
78  allocate(fct(rts(rmax)+ncoefs(rmax-8,6)+rmax+1))
79  call readcache(fct,rts(rmax)+ncoefs(rmax-8,6)+rmax+1,x,41,10+mode_cll,0,6,rank,nocalc,wrica)
80  else
81  allocate(fct(rts(rmax)+rmax+1))
82  call readcache(fct,rts(rmax)+rmax+1,x,41,10+mode_cll,0,6,rank,nocalc,wrica)
83  end if
84 
85  if(nocalc)then
86  tfuv = 0d0
87  cnt = 1
88  do r=0,rmax
89  do n0=0,r
90  do n1=0,r-n0
91  do n2=0,r-n0-n1
92  n3 = r-n0-n1-n2
93  tf(n0,n1,n2,n3) = fct(cnt)
94  cnt = cnt+1
95  end do
96  end do
97  end do
98  do n0=4,r/2
99  do n1=0,r-2*n0
100  do n2=0,r-2*n0-n1
101  do n3=0,r-2*n0-n1-n2
102  do n4=0,r-2*n0-n1-n2-n3
103  n5 = r-2*n0-n1-n2-n3-n4
104 
105  cfuv(n0,n1,n2,n3,n4,n5) = fct(cnt)
106  cnt = cnt+1
107 
108  end do
109  end do
110  end do
111  end do
112  end do
113  tferr(r) = real(fct(cnt))
114  cnt = cnt+1
115  end do
116  tfuv = 0d0
117  if (rmax.ge.8) call calctensorfuv(tfuv,cfuv,momvec,rmax)
118  return
119  endif
120 
121 
122 
123  if(rank.eq.rmax) then
124 
125  call calctensorfrred(tf,cfuv,tferr,momvec,mominv,masses2,rank)
126 
127  if (wrica) then
128  cnt = 0
129  do r=0,rank
130  do n0=0,r
131  do n1=0,r-n0
132  do n2=0,r-n0-n1
133  n3 = r-n0-n1-n2
134 
135  cnt = cnt+1
136  fct(cnt) = tf(n0,n1,n2,n3)
137 
138  end do
139  end do
140  end do
141  do n0=4,r/2
142  do n1=0,r-2*n0
143  do n2=0,r-2*n0-n1
144  do n3=0,r-2*n0-n1-n2
145  do n4=0,r-2*n0-n1-n2-n3
146  n5 = r-2*n0-n1-n2-n3-n4
147 
148  cnt = cnt+1
149  fct(cnt) = cfuv(n0,n1,n2,n3,n4,n5)
150 
151  end do
152  end do
153  end do
154  end do
155  end do
156  cnt = cnt+1
157  fct(cnt) = tferr(r)
158  end do
159 
160  if (rank.ge.8) then
161  call writecache(fct,rts(rank)+ncoefs(rank-8,6)+rank+1,0,6,rank)
162  else
163  call writecache(fct,rts(rank)+rank+1,0,6,rank)
164  end if
165 
166  end if
167 
168  tfuv = 0d0
169  if (rmax.ge.8) call calctensorfuv(tfuv,cfuv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax),momvec,rmax)
170  return
171 
172 
173  else
174  allocate(tfaux(0:rank,0:rank,0:rank,0:rank))
175  allocate(cfuvaux(0:rank/2,0:rank,0:rank,0:rank,0:rank,0:rank))
176  allocate(tferraux(0:rank))
177 
178  call calctensorfrred(tfaux,cfuvaux,tferraux,momvec,mominv,masses2,rank)
179 
180  if (wrica) then
181  cnt = 0
182  deallocate(fct)
183  if (rank.ge.8) then
184  allocate(fct(rts(rank)+ncoefs(rank-8,6)+rank+1))
185  else
186  allocate(fct(rts(rank)+rank+1))
187  end if
188  do r=0,rank
189  do n0=0,r
190  do n1=0,r-n0
191  do n2=0,r-n0-n1
192  n3 = r-n0-n1-n2
193 
194  cnt = cnt+1
195  fct(cnt) = tfaux(n0,n1,n2,n3)
196 
197  end do
198  end do
199  end do
200  do n0=4,r/2
201  do n1=0,r-2*n0
202  do n2=0,r-2*n0-n1
203  do n3=0,r-2*n0-n1-n2
204  do n4=0,r-2*n0-n1-n2-n3
205  n5 = r-2*n0-n1-n2-n3-n4
206 
207  cnt = cnt+1
208  fct(cnt) = cfuvaux(n0,n1,n2,n3,n4,n5)
209 
210  end do
211  end do
212  end do
213  end do
214  end do
215  cnt = cnt+1
216  fct(cnt) = tferraux(r)
217  end do
218 
219  if (rank.ge.8) then
220  call writecache(fct,rts(rank)+ncoefs(rank-8,6)+rank+1,0,6,rank)
221  else
222  call writecache(fct,rts(rank)+rank+1,0,6,rank)
223  end if
224 
225  end if
226 
227  tf = tfaux(0:rmax,0:rmax,0:rmax,0:rmax)
228  tfuv = 0d0
229  if (rmax.ge.8) call calctensorfuv(tfuv,cfuvaux(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax),momvec,rmax)
230  tferr = tferraux(0:rmax)
231  return
232 
233  end if
234 
235 ! end if
236  end if
237  end if
238 
239  allocate(cfuvaux(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax))
240  call calctensorfrred(tf,cfuvaux,tferr,momvec,mominv,masses2,rmax)
241  tfuv = 0d0
242  if (rmax.ge.8) call calctensorfuv(tfuv,cfuvaux,momvec,rmax)
243 
244 
245  end subroutine calctensorfr
246 
247 
248 
249 
250 
251  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
252  ! subroutine CalcTensorFrRed(TF,TFuv,TFerr,MomVec,MomInv,masses2,rmax,CFuv)
253  !
254  ! calculate tensor integral TFr for 6-point functions by direct reduction
255  ! from the tensor integrals TEr for 5-point functions.
256  !
257  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
258 
259  subroutine calctensorfrred(TF,CFuv,TFerr,MomVec,MomInv,masses2,rmax)
260 
261  integer, intent(in) :: rmax
262  double complex, intent(in) :: MomVec(0:3,5), masses2(0:5), MomInv(15)
263  double complex, intent(out) :: TF(0:rmax,0:rmax,0:rmax,0:rmax)
264  double precision, intent(out) :: TFerr(0:rmax)
265  double complex,intent(out) :: CFuv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
266  double complex, allocatable :: TE(:,:,:,:,:), TE_0aux(:,:,:,:)
267  double complex, allocatable :: TEuv(:,:,:,:,:), TEuv_0aux(:,:,:,:)
268  double precision, allocatable :: TEerr(:,:)
269  double precision :: p1max
270  double complex, allocatable :: CE(:,:,:,:,:), CEuv(:,:,:,:,:), CE0uv(:,:,:,:,:,:)
271  double precision, allocatable :: CEerr(:)
272  double complex :: MomVecE(0:3,4), masses2E(0:4), MomInvE(10)
273  double complex :: q10,q21,q32,q43,q54,q50,q20,q31,q42,q53,q40
274  double complex :: q51,q30,q41,q52,mm02,mm12,mm22,mm32,mm42,mm52,ff(5)
275  double complex :: mx(0:5,0:5), mxinv(0:5,0:5),mx0k(5,5),mx0kinv(5,5),mxinvP(5,0:3)
276  double complex :: mxinvs,mxinvPs(0:3)
277  double complex :: det,newdet
278  double complex :: Smod, Faux(5), elimminf2_coli,chdet,P1ten
279  integer :: r,n0,n1,n2,n3,n4,n5,np0,np1,np2,np3,k,i,j,kbest,rmaxE,mu,combi,nid(0:4),bin
280 
281 
282  ! determine inverse modified Caley matrix
283  mm02 = elimminf2_coli(masses2(0))
284  mm12 = elimminf2_coli(masses2(1))
285  mm22 = elimminf2_coli(masses2(2))
286  mm32 = elimminf2_coli(masses2(3))
287  mm42 = elimminf2_coli(masses2(4))
288  mm52 = elimminf2_coli(masses2(5))
289  q10 = elimminf2_coli(mominv(1))
290  q21 = elimminf2_coli(mominv(2))
291  q32 = elimminf2_coli(mominv(3))
292  q43 = elimminf2_coli(mominv(4))
293  q54 = elimminf2_coli(mominv(5))
294  q50 = elimminf2_coli(mominv(6))
295  q20 = elimminf2_coli(mominv(7))
296  q31 = elimminf2_coli(mominv(8))
297  q42 = elimminf2_coli(mominv(9))
298  q53 = elimminf2_coli(mominv(10))
299  q40 = elimminf2_coli(mominv(11))
300  q51 = elimminf2_coli(mominv(12))
301  q30 = elimminf2_coli(mominv(13))
302  q41 = elimminf2_coli(mominv(14))
303  q52 = elimminf2_coli(mominv(15))
304 
305  ff(1) = q10+mm02-mm12
306  ff(2) = q20+mm02-mm22
307  ff(3) = q30+mm02-mm32
308  ff(4) = q40+mm02-mm42
309  ff(5) = q50+mm02-mm52
310 
311 
312  mx(0,0) = 2d0*mm02
313  mx(1,0) = q10 - mm12 + mm02
314  mx(2,0) = q20 - mm22 + mm02
315  mx(3,0) = q30 - mm32 + mm02
316  mx(4,0) = q40 - mm42 + mm02
317  mx(5,0) = q50 - mm52 + mm02
318  mx(0,1) = mx(1,0)
319  mx(1,1) = 2d0*q10
320  mx(2,1) = q10+q20-q21
321  mx(3,1) = q10+q30-q31
322  mx(4,1) = q10+q40-q41
323  mx(5,1) = q10+q50-q51
324  mx(0,2) = mx(2,0)
325  mx(1,2) = mx(2,1)
326  mx(2,2) = 2d0*q20
327  mx(3,2) = q20+q30-q32
328  mx(4,2) = q20+q40-q42
329  mx(5,2) = q20+q50-q52
330  mx(0,3) = mx(3,0)
331  mx(1,3) = mx(3,1)
332  mx(2,3) = mx(3,2)
333  mx(3,3) = 2d0*q30
334  mx(4,3) = q30+q40-q43
335  mx(5,3) = q30+q50-q53
336  mx(0,4) = mx(4,0)
337  mx(1,4) = mx(4,1)
338  mx(2,4) = mx(4,2)
339  mx(3,4) = mx(4,3)
340  mx(4,4) = 2d0*q40
341  mx(5,4) = q40+q50-q54
342  mx(0,5) = mx(5,0)
343  mx(1,5) = mx(5,1)
344  mx(2,5) = mx(5,2)
345  mx(3,5) = mx(5,3)
346  mx(4,5) = mx(5,4)
347  mx(5,5) = 2d0*q50
348 
349  call chinv(6,mx,mxinv)
350 
351 
352  ! determine X_(0,5)
353  do j=1,5
354  do i=1,5
355  mx0k(i,j) = mx(i,j-1)
356  end do
357  end do
358 
359  det = chdet(5,mx0k)
360  kbest = 5
361 
362  do j=5,2,-1
363  do i=1,5
364  mx0k(i,j) = mx(i,j)
365  end do
366 
367  newdet = chdet(5,mx0k)
368  if (abs(newdet).gt.abs(det)) then
369  kbest = j-1
370  det = newdet
371  end if
372 
373  end do
374 
375  do i=1,5
376  mx0k(i,1) = mx(i,1)
377  mx0k(i,kbest) = mx(i,0)
378  end do
379 
380  call chinv(5,mx0k,mx0kinv)
381  do i=1,5
382  mx0kinv(kbest,i) = 0d0
383  end do
384 
385  mxinvs = sum(mxinv(0:5,0))
386 
387  ! build rank5 tensors
388  rmaxe = max(rmax-1,0)
389  allocate(ce0uv(0:rmaxe/2,0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe))
390  allocate(ce(0:rmaxe/2,0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe))
391  allocate(ceuv(0:rmaxe/2,0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe))
392  allocate(ceerr(0:rmaxe))
393  allocate(te(0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe,0:5))
394  allocate(te_0aux(0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe))
395  allocate(teuv(0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe,0:5))
396  allocate(teuv_0aux(0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe))
397  allocate(teerr(0:5,0:rmaxe))
398 
399  mominve = submominv(6,0,mominv)
400  masses2e = submasses(6,0,masses2)
401  call e_main_cll(ce(:,:,:,:,:),ce0uv(:,0,:,:,:,:), &
402  mominve(1),mominve(2),mominve(3),mominve(4),mominve(5),mominve(6), &
403  mominve(7),mominve(8),mominve(9),mominve(10),masses2e(0),masses2e(1), &
404  masses2e(2),masses2e(3),masses2e(4),rmaxe,eerr2=ceerr,id_in=1)
405  call calctensore(te_0aux(:,:,:,:),teuv_0aux(:,:,:,:),teerr(0,0:), &
406  ce(:,:,:,:,:),ce0uv(:,0,:,:,:,:),ceerr,submomvec(6,0,momvec),rmaxe)
407 
408  ! shift of integration momentum in TE\{0}
409  te(:,:,:,:,0) = 0d0
410 
411  do r=0,rmaxe
412 
413  do np0=0,r
414  do np1=0,r-np0
415  do np2=0,r-np0-np1
416  do np3=0,r-np0-np1-np2
417 
418  p1ten = (-momvec(0,1))**np0*(-momvec(1,1))**np1* &
419  (-momvec(2,1))**np2*(-momvec(3,1))**np3
420 
421  do n0=0,r-np0-np1-np2-np3
422  do n1=0,r-np0-np1-np2-np3-n0
423  do n2=0,r-np0-np1-np2-np3-n0-n1
424  n3 = r-np0-np1-np2-np3-n0-n1-n2
425 
426  combi = binomtable(n0,n0+np0)*binomtable(n1,n1+np1)* &
427  binomtable(n2,n2+np2)*binomtable(n3,n3+np3)
428 
429  te(np0+n0,np1+n1,np2+n2,np3+n3,0) = te(np0+n0,np1+n1,np2+n2,np3+n3,0) &
430  + combi*te_0aux(n0,n1,n2,n3)*p1ten
431 
432  end do
433  end do
434  end do
435 
436  end do
437  end do
438  end do
439  end do
440 
441  end do
442 
443  ! shift of integration momentum in TEuv\{0} coefficients
444  do n1=1,rmaxe-6
445  do n2=0,rmaxe-n1-6
446  do n3=0,rmaxe-n1-n2-6
447  do n4=0,rmaxe-n1-n2-n3-6
448  do n5=0,rmaxe-n1-n2-n3-n4-6
449  n0 = (rmaxe-n1-n2-n3-n4-n5)/2
450  ce0uv(0:n0,n1,n2,n3,n4,n5) = -ce0uv(0:n0,n1-1,n2,n3,n4,n5)-ce0uv(0:n0,n1-1,n2+1,n3,n4,n5) &
451  -ce0uv(0:n0,n1-1,n2,n3+1,n4,n5)-ce0uv(0:n0,n1-1,n2,n3,n4+1,n5) &
452  -ce0uv(0:n0,n1-1,n2,n3,n4,n5+1)
453  end do
454  end do
455  end do
456  end do
457  end do
458 
459 
460  p1max = maxval(abs(momvec(0:3,1)))
461  do r=1,rmaxe
462  do i=1,r
463  teerr(0,r) = max(teerr(0,i),teerr(0,r-i)*p1max**i)
464  end do
465  end do
466 
467 
468  do i=1,5
469  mominve = submominv(6,i,mominv)
470  masses2e = submasses(6,i,masses2)
471  call e_main_cll(ce,ceuv,mominve(1),mominve(2),mominve(3),mominve(4),mominve(5),mominve(6), &
472  mominve(7),mominve(8),mominve(9),mominve(10),masses2e(0),masses2e(1), &
473  masses2e(2),masses2e(3),masses2e(4),rmaxe,eerr2=ceerr,id_in=2**i)
474  call calctensore(te(:,:,:,:,i),teuv(:,:,:,:,i),teerr(i,0:), &
475  ce,ceuv,ceerr,submomvec(6,i,momvec),rmaxe)
476  end do
477 
478  ! UV divergent tensor
479  ! calculation of UV-divergent parts
480  ! F_(n0,n1,n2,n3) UV-finite for n0<4
481  cfuv(0:min(rmax/2,3),:,:,:,:,:) = 0d0
482  ! PV reduction (5.10)
483  do r=8,rmax
484  do n0=3,r/2
485  do n1=0,r-2*n0
486  do n2=0,r-2*n0-n1
487  do n3=0,r-2*n0-n1-n2
488  do n4=0,r-2*n0-n1-n2-n3
489  n5 = r-2*n0-n1-n2-n3-n4
490 
491  cfuv(n0,n1,n2,n3,n4,n5) = (ce0uv(n0-1,n1,n2,n3,n4,n5) &
492  + 2*masses2(0)*cfuv(n0-1,n1,n2,n3,n4,n5) &
493  + ff(1)*cfuv(n0-1,n1+1,n2,n3,n4,n5) &
494  + ff(2)*cfuv(n0-1,n1,n2+1,n3,n4,n5) &
495  + ff(3)*cfuv(n0-1,n1,n2,n3+1,n4,n5) &
496  + ff(4)*cfuv(n0-1,n1,n2,n3,n4+1,n5) &
497  + ff(5)*cfuv(n0-1,n1,n2,n3,n4,n5+1))/ (2*(r-3))
498 
499  end do
500  end do
501  end do
502  end do
503  end do
504  end do
505 
506 
507  tf = 0d0
508  tf(0,0,0,0) = -mxinv(0,0)*te(0,0,0,0,0)
509  do i=1,5
510  tf(0,0,0,0) = tf(0,0,0,0) + mxinv(i,0)*(te(0,0,0,0,i)-te(0,0,0,0,0))
511  end do
512  tferr(0) = max( abs(mxinvs)*teerr(0,0), &
513  abs(mxinv(1,0))*teerr(1,0) , &
514  abs(mxinv(2,0))*teerr(2,0) , &
515  abs(mxinv(3,0))*teerr(3,0) , &
516  abs(mxinv(4,0))*teerr(4,0) , &
517  abs(mxinv(5,0))*teerr(5,0) )
518 
519  do mu=0,3
520  do i=1,5
521  mxinvp(i,mu)=0d0
522  do j=1,5
523  mxinvp(i,mu) = mxinvp(i,mu) + mx0kinv(j,i)*momvec(mu,j)
524  end do
525  end do
526  mxinvps(mu) = sum(mxinvp(1:5,mu))
527  end do
528 
529 
530 
531 
532  do r=1,rmax
533  do n0=0,r
534  do n1=0,r-n0
535  do n2=0,r-n0-n1
536  n3 = r-n0-n1-n2
537 
538  if (n0.ge.1) then
539  do i=1,5
540  smod = te(n0-1,n1,n2,n3,i) - te(n0-1,n1,n2,n3,0)
541  tf(n0,n1,n2,n3) = tf(n0,n1,n2,n3) + mxinvp(i,0)*smod
542  end do
543  else if (n1.ge.1) then
544  do i=1,5
545  smod = te(n0,n1-1,n2,n3,i) - te(n0,n1-1,n2,n3,0)
546  tf(n0,n1,n2,n3) = tf(n0,n1,n2,n3) + mxinvp(i,1)*smod
547  end do
548  else if (n2.ge.1) then
549  do i=1,5
550  smod = te(n0,n1,n2-1,n3,i) - te(n0,n1,n2-1,n3,0)
551  tf(n0,n1,n2,n3) = tf(n0,n1,n2,n3) + mxinvp(i,2)*smod
552  end do
553  else
554  do i=1,5
555  smod = te(n0,n1,n2,n3-1,i) - te(n0,n1,n2,n3-1,0)
556  tf(n0,n1,n2,n3) = tf(n0,n1,n2,n3) + mxinvp(i,3)*smod
557  end do
558  end if
559 
560  end do
561  end do
562  end do
563 
564  tferr(r) = max(maxval(abs(mxinvps(0:3)))*teerr(0,r-1), &
565  maxval(abs(mxinvp(1,0:3)))*teerr(1,r-1), &
566  maxval(abs(mxinvp(2,0:3)))*teerr(2,r-1), &
567  maxval(abs(mxinvp(3,0:3)))*teerr(3,r-1))
568 
569  end do
570 
571 
572  end subroutine calctensorfrred
573 
574 
575 
576 
577 
578  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
579  ! subroutine CalcTensorFr_list(TF,TFuv,TFerr,MomVec,MomInv,masses2,rmax)
580  !
581  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
582 
583  subroutine calctensorfr_list(TF,TFuv,TFerr,MomVec,MomInv,masses2,rmax)
584 
585  integer, intent(in) :: rmax
586  double complex, intent(in) :: MomVec(0:3,5), masses2(0:5), MomInv(15)
587  double complex, intent(out) :: TF(RtS(rmax)), TFuv(RtS(rmax))
588  double precision, intent(out) :: TFerr(0:rmax)
589  double complex :: TF_aux(0:rmax,0:rmax,0:rmax,0:rmax)
590  double complex :: TFuv_aux(0:rmax,0:rmax,0:rmax,0:rmax)
591  integer :: mu
592 
593  call calctensorfr(tf_aux,tfuv_aux,tferr,momvec,mominv,masses2,rmax)
594 
595  do mu=1,rts(rmax)
596  tf(mu) = tf_aux(lorindtab(0,mu),lorindtab(1,mu),lorindtab(2,mu),lorindtab(3,mu))
597  end do
598 
599  if (calcuv_cll) then
600  do mu=1,rts(rmax)
601  tfuv(mu) = tfuv_aux(lorindtab(0,mu),lorindtab(1,mu),lorindtab(2,mu),lorindtab(3,mu))
602  end do
603  end if
604 
605  end subroutine calctensorfr_list
606 
607 
608 
609 
610 
611  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
612  ! subroutine CalcTensorTNr(TN,TNuv,TNerr,MomVec,MomInv,masses2,rmax,id,CNuv)
613  !
614  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
615 
616  recursive subroutine calctensortnr(TN,TNuv,TNerr,MomVec,MomInv,masses2,N,rmax,id,CNuv)
617 
618  integer, intent(in) :: n,rmax,id
619  double complex, intent(in) :: momvec(0:3,n-1), masses2(0:n-1)
620  double complex, intent(in) :: mominv(binomtable(2,n))
621  double complex, intent(out) :: tn(0:rmax,0:rmax,0:rmax,0:rmax)
622  double complex, intent(out) :: tnuv(0:rmax,0:rmax,0:rmax,0:rmax)
623  double complex, intent(out), optional :: cnuv(binomtable(rmax-2*n+4,max(rmax-n+2,0)),n-2:rmax/2,2*n-4:rmax)
624  double complex :: cnuvaux0(binomtable(max(rmax-2*n+4,0),max(rmax-n+2,0)),n-2:max(rmax/2,n-2),2*n-4:max(rmax,2*n-4))
625  double precision, intent(out) :: tnerr(0:rmax)
626  double complex, allocatable :: tnaux(:,:,:,:),tnuvaux(:,:,:,:)
627  double precision, allocatable :: tnerraux(:),cnerr(:)
628  double complex :: x(binomtable(2,n)+n+4*(n-1))
629  double complex, allocatable :: fct(:),cn(:),cnuvaux1(:),cnuvaux2(:,:,:)
630  integer :: r,n0,n1,n2,n3,n4,i,rank,bino,cnt,cntc
631  logical :: nocalc,wrica
632 
633  if ((use_cache_system).and.(n.ge.tencache)) then
634  if ((ncache.gt.0).and.(ncache.le.ncache_max)) then
635 ! if (use_cache(ncache).ge.N) then
636  do i=1,n-1
637  x((i-1)*4+1:i*4) = momvec(0:3,i)
638  end do
639  bino = binomtable(2,n)
640  x(4*n-3:bino+4*n-4) = mominv
641  x(4*n-3+bino:bino+5*n-4) = masses2
642  rank = rmax
643 
644  if (rmax.ge.2*n-4) then
645  allocate(fct(rts(rmax)+ncoefs(rmax-2*n+4,n)+rmax+1))
646  call readcache(fct,rts(rmax)+ncoefs(rmax-2*n+4,n)+rmax+1,x,bino+5*n-4,mode_cll+10,id,n,rank,nocalc,wrica)
647  else
648  allocate(fct(rts(rmax)+rmax+1))
649  call readcache(fct,rts(rmax)+rmax+1,x,bino+5*n-4,mode_cll+10,id,n,rank,nocalc,wrica)
650  end if
651  ! write(*,*) N,id, nocalc,wrica, masses2(0)
652 
653  if(nocalc)then
654  cnt = 1
655  do r=0,rmax
656  do n0=0,r
657  do n1=0,r-n0
658  do n2=0,r-n0-n1
659  n3 = r-n0-n1-n2
660  tn(n0,n1,n2,n3) = fct(cnt)
661  cnt = cnt+1
662  end do
663  end do
664  end do
665  if(id.eq.0) then
666  do n0=r/2,n-2,-1
667  do i=1,binomtable(r-2*n0,max(n+r-2*n0-2,0))
668  cnuvaux0(i,n0,r) = fct(cnt)
669  cnt = cnt+1
670  end do
671  end do
672  else
673  do n0=r/2,n-2,-1
674  do i=1,binomtable(r-2*n0,max(n+r-2*n0-2,0))
675  cnuv(i,n0,r) = fct(cnt)
676  cnt = cnt+1
677  end do
678  end do
679  end if
680  tnerr(r) = real(fct(cnt))
681  cnt = cnt+1
682  end do
683 ! write(*,*) 'nocalc', N, rmax, cnt-1, fct(4)
684  if(id.eq.0) then
685  if(rmax.ge.2*n-4) then
686  call calctensortnuv_list(tnuv,cnuvaux0(1:binomtable(rmax-2*n+4,max(rmax-n+2,0)), &
687  n-2:rmax/2,2*n-4:rmax),momvec,n,rmax)
688  else
689  tnuv=0d0
690  end if
691  end if
692  return
693  endif
694 
695 
696  if(rank.eq.rmax) then
697 
698  if (n.eq.tenred_cll-1) then
699  allocate(cn(ncoefs(rank,n)))
700  allocate(cnuvaux1(ncoefs(rank,n)))
701  allocate(cnerr(0:rank))
702  allocate(tnuvaux(0:rank,0:rank,0:rank,0:rank))
703  call tn_cll(cn,cnuvaux1,mominv,masses2,n,rank,tnerr2=cnerr,id_in=id)
704  call calctensortn(tn,tnuv,tnerr,cn,cnuvaux1,cnerr,momvec,n,rank)
705  if((id.ne.0).and.(rmax.ge.2*n-4)) then
706  do r=2*n-4,rmax
707  cnt = ncoefs(r-1,n)+1
708  do n0=r/2,n-2,-1
709  do i=1,binomtable(r-2*n0,max(n+r-2*n0-2,0))
710  cnuv(i,n0,r) = cnuvaux1(cnt)
711  cnt = cnt+1
712  end do
713  end do
714  end do
715  end if
716  else
717  if(rank.ge.2*n-4) then
718  if(id.eq.0) then
719  call calctensortnrred(tn,tnerr,momvec,mominv,masses2,n,rank,id,cnuvaux0)
720  call calctensortnuv_list(tnuv,cnuvaux0,momvec,n,rank)
721  else
722  call calctensortnrred(tn,tnerr,momvec,mominv,masses2,n,rank,id,cnuv)
723  end if
724  else
725  call calctensortnrred(tn,tnerr,momvec,mominv,masses2,n,rank,id)
726  if(id.eq.0) tnuv=0d0
727  end if
728  end if
729 
730  if (wrica) then
731  cnt = 0
732  do r=0,rank
733  do n0=0,r
734  do n1=0,r-n0
735  do n2=0,r-n0-n1
736  n3 = r-n0-n1-n2
737  cnt = cnt+1
738  fct(cnt) = tn(n0,n1,n2,n3)
739  end do
740  end do
741  end do
742  if(r.ge.2*n-4) then
743  if(n.eq.tenred_cll-1) then
744  do i=ncoefs(r-2*n+3,n)+1,ncoefs(r-2*n+4,n)
745  cnt = cnt+1
746  fct(cnt) = cnuvaux1(i)
747  end do
748  else
749  if(id.eq.0) then
750  do n0=r/2,n-2,-1
751  do i=1,binomtable(r-2*n0,max(n+r-2*n0-2,0))
752  cnt = cnt+1
753  fct(cnt) = cnuvaux0(i,n0,r)
754  end do
755  end do
756  else
757  do n0=r/2,n-2,-1
758  do i=1,binomtable(r-2*n0,max(n+r-2*n0-2,0))
759  cnt = cnt+1
760  fct(cnt) = cnuv(i,n0,r)
761  end do
762  end do
763  end if
764  end if
765  end if
766 
767  cnt = cnt+1
768  fct(cnt) = tnerr(r)
769 
770  end do
771 ! write(*,*) 'wrica', N, rank, cnt, fct(4)
772 
773  if(rank.ge.2*n-4) then
774  call writecache(fct,rts(rank)+ncoefs(rank-2*n+4,n)+rank+1,id,n,rank)
775  else
776  call writecache(fct,rts(rank)+rank+1,id,n,rank)
777  end if
778 
779  end if
780  return
781 
782 
783  else
784  allocate(tnaux(0:rank,0:rank,0:rank,0:rank))
785  allocate(tnerraux(0:rank))
786 
787  if (n.eq.tenred_cll-1) then
788  allocate(cn(ncoefs(rank,n)))
789  allocate(cnuvaux1(ncoefs(rank,n)))
790  allocate(cnerr(0:rank))
791  allocate(tnuvaux(0:rank,0:rank,0:rank,0:rank))
792  call tn_cll(cn,cnuvaux1,mominv,masses2,n,rank,tnerr2=cnerr,id_in=id)
793  call calctensortn(tnaux,tnuvaux,tnerraux,cn,cnuvaux1,cnerr,momvec,n,rank)
794  if(id.eq.0) then
795  tnuv = tnuvaux(0:rmax,0:rmax,0:rmax,0:rmax)
796  else if(rmax.ge.2*n-4) then
797  do r=2*n-4,rmax
798  cnt = ncoefs(r-1,n)+1
799  do n0=r/2,n-2,-1
800  do i=1,binomtable(r-2*n0,max(n+r-2*n0-2,0))
801  cnuv(i,n0,r) = cnuvaux1(cnt)
802  cnt = cnt+1
803  end do
804  end do
805  end do
806  end if
807  else
808  if(rank.ge.2*n-4) then
809  allocate(cnuvaux2(binomtable(rank-2*n+4,max(rank-n+2,0)),n-2:rank/2,2*n-4:rank))
810  call calctensortnrred(tnaux,tnerraux,momvec,mominv,masses2,n,rank,id,cnuvaux2)
811  if(id.eq.0) then
812  tnuv = 0d0
813  call calctensortnuv_list(tnuv,cnuvaux2(1:binomtable(rmax-2*n+4,max(rmax-n+2,0)), &
814  n-2:rmax/2,2*n-4:rmax),momvec,n,rmax)
815  else
816  cnuv(1:binomtable(rmax-2*n+4,max(rmax-n+2,0)),n-2:rmax/2,2*n-4:rmax) = &
817  cnuvaux2(1:binomtable(rmax-2*n+4,max(rmax-n+2,0)),n-2:rmax/2,2*n-4:rmax)
818  end if
819  else
820  call calctensortnrred(tnaux,tnerraux,momvec,mominv,masses2,n,rank,id)
821  if(id.eq.0) tnuv=0d0
822  end if
823  end if
824 
825  if (wrica) then
826  cnt = 0
827  cntc = 0
828  deallocate(fct)
829  if(rank.ge.2*n-4) then
830  allocate(fct(rts(rank)+ncoefs(rank-2*n+4,n)+rank+1))
831  else
832  allocate(fct(rts(rank)+rank+1))
833  end if
834  do r=0,rank
835  do n0=0,r
836  do n1=0,r-n0
837  do n2=0,r-n0-n1
838  n3 = r-n0-n1-n2
839  cnt = cnt+1
840  fct(cnt) = tnaux(n0,n1,n2,n3)
841  end do
842  end do
843  end do
844  if(r.ge.2*n-4) then
845  if(n.eq.tenred_cll-1) then
846  do i=ncoefs(r-2*n+3,n)+1,ncoefs(r-2*n+4,n)
847  cnt = cnt+1
848  fct(cnt) = cnuvaux1(i)
849  end do
850  else
851  do n0=r/2,n-2,-1
852  do i=1,binomtable(r-2*n0,max(n+r-2*n0-2,0))
853  cnt = cnt+1
854  fct(cnt) = cnuvaux2(i,n0,r)
855  end do
856  end do
857  end if
858  end if
859 
860  cnt = cnt+1
861  fct(cnt) = tnerraux(r)
862 
863  end do
864 
865  if(rank.ge.2*n-4) then
866  call writecache(fct,rts(rank)+ncoefs(rank-2*n+4,n)+rank+1,id,n,rank)
867  else
868  call writecache(fct,rts(rank)+rank+1,id,n,rank)
869  end if
870 
871  end if
872 
873  tn = tnaux(0:rmax,0:rmax,0:rmax,0:rmax)
874  tnerr = tnerraux(0:rmax)
875  return
876 
877  end if
878 
879 ! end if
880  end if
881  end if
882 
883  if (n.le.tenred_cll-1) then
884  allocate(cn(ncoefs(rmax,n)))
885  allocate(cnuvaux1(ncoefs(rmax,n)))
886  allocate(cnerr(0:rmax))
887  call tn_cll(cn,cnuvaux1,mominv,masses2,n,rmax,tnerr2=cnerr,id_in=id)
888  call calctensortn(tn,tnuv,tnerr,cn,cnuvaux1,cnerr,momvec,n,rmax)
889  if((id.ne.0).and.(rmax.ge.2*n-4)) then
890  do r=2*n-4,rmax
891  cnt = ncoefs(r-1,n)+1
892  do n0=r/2,n-2,-1
893  do i=1,binomtable(r-2*n0,max(n+r-2*n0-2,0))
894  cnuv(i,n0,r) = cnuvaux1(cnt)
895  cnt = cnt+1
896  end do
897  end do
898  end do
899  end if
900  else
901  if(rmax.ge.2*n-4) then
902  if(id.eq.0) then
903  call calctensortnrred(tn,tnerr,momvec,mominv,masses2,n,rmax,id,cnuvaux0)
904  call calctensortnuv_list(tnuv,cnuvaux0,momvec,n,rmax)
905  else
906  call calctensortnrred(tn,tnerr,momvec,mominv,masses2,n,rmax,id,cnuv)
907  end if
908  else
909  call calctensortnrred(tn,tnerr,momvec,mominv,masses2,n,rmax,id)
910  if(id.eq.0) tnuv=0d0
911  end if
912  end if
913 
914 
915  end subroutine calctensortnr
916 
917 
918 
919 
920 
921  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
922  ! subroutine CalcTensorTNrRed(TN,TNerr,MomVec,MomInv,masses2,rmax,id,CNuv)
923  !
924  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
925 
926  recursive subroutine calctensortnrred(TN,TNerr,MomVec,MomInv,masses2,N,rmax,id,CNuv)
927 
928  integer, intent(in) :: n,rmax,id
929  double complex, intent(in) :: momvec(0:3,n-1), masses2(0:n-1)
930  double complex, intent(in) :: mominv(binomtable(2,n))
931  double complex, intent(out) :: tn(0:rmax,0:rmax,0:rmax,0:rmax)
932  double complex, intent(out), optional :: cnuv(binomtable(rmax-2*n+4,max(rmax-n+2,0)),n-2:rmax/2,2*n-4:rmax)
933  double precision, intent(out) :: tnerr(0:rmax)
934  double complex, allocatable :: tnm1(:,:,:,:,:),tnm1_0aux(:,:,:,:)
935  double complex, allocatable :: tnm1uv_0aux(:,:,:,:)
936  double complex, allocatable :: cnm1uv_0(:,:,:),cnm1uv_0aux(:,:,:),cnm1uv_i(:,:,:)
937  double precision, allocatable :: tnm1err(:,:)
938  double precision :: p1max
939  double complex :: q10,q21,q32,q43,q54,q50,q20,q31,q42,q53,q40
940  double complex :: q51,q30,q41,q52,mm02,mm12,mm22,mm32,mm42,mm52
941  double complex :: mx(0:5,0:5), mxinv(0:5,0:5),mx0k(5,5),mx0kinv(5,5)
942  double complex :: det,newdet,ff(n-1)
943  double complex :: smod, elimminf2_coli,chdet,p1ten,mxinvp(5,0:3)
944  double complex :: mxinvs,mxinvps(0:3)
945  integer :: r,n0,n1,n2,n3,n4,n5,np0,np1,np2,np3,k,i,j,kbest,rmax_m1,combi,mu,cnt
946  integer :: bin,nid(0:5),r0,rbcd,mia,bino_0,bino_i,cind,mia1,mia2,mia3
947 
948 
949  ! determine inverse modified Caley matrix
950  mm02 = elimminf2_coli(masses2(0))
951  mm12 = elimminf2_coli(masses2(1))
952  mm22 = elimminf2_coli(masses2(2))
953  mm32 = elimminf2_coli(masses2(3))
954  mm42 = elimminf2_coli(masses2(4))
955  mm52 = elimminf2_coli(masses2(5))
956  q10 = elimminf2_coli(mominv(1))
957  q21 = elimminf2_coli(mominv(2))
958  q32 = elimminf2_coli(mominv(3))
959  q43 = elimminf2_coli(mominv(4))
960  q54 = elimminf2_coli(mominv(5))
961  if (n.le.9) then
962  q50 = elimminf2_coli(mominv((n-6)*n+6))
963  else
964  q50 = elimminf2_coli(mominv(4*n+1))
965  end if
966  q20 = elimminf2_coli(mominv(n+1))
967  q31 = elimminf2_coli(mominv(n+2))
968  q42 = elimminf2_coli(mominv(n+3))
969  q53 = elimminf2_coli(mominv(n+4))
970  if (n.le.7) then
971  q40 = elimminf2_coli(mominv((n-5)*n+5))
972  q51 = elimminf2_coli(mominv((n-5)*n+6))
973  else
974  q40 = elimminf2_coli(mominv(3*n+1))
975  q51 = elimminf2_coli(mominv(3*n+2))
976  end if
977  q30 = elimminf2_coli(mominv(2*n+1))
978  q41 = elimminf2_coli(mominv(2*n+2))
979  q52 = elimminf2_coli(mominv(2*n+3))
980 
981  cnt = 1
982  do i=1,n-1
983  ff(i) = elimminf2_coli(mominv(cnt)) + mm02 &
984  - elimminf2_coli(masses2(i))
985  cnt = cnt+n-i
986  end do
987 
988  mx(0,0) = 2d0*mm02
989  mx(1,0) = q10 - mm12 + mm02
990  mx(2,0) = q20 - mm22 + mm02
991  mx(3,0) = q30 - mm32 + mm02
992  mx(4,0) = q40 - mm42 + mm02
993  mx(5,0) = q50 - mm52 + mm02
994  mx(0,1) = mx(1,0)
995  mx(1,1) = 2d0*q10
996  mx(2,1) = q10+q20-q21
997  mx(3,1) = q10+q30-q31
998  mx(4,1) = q10+q40-q41
999  mx(5,1) = q10+q50-q51
1000  mx(0,2) = mx(2,0)
1001  mx(1,2) = mx(2,1)
1002  mx(2,2) = 2d0*q20
1003  mx(3,2) = q20+q30-q32
1004  mx(4,2) = q20+q40-q42
1005  mx(5,2) = q20+q50-q52
1006  mx(0,3) = mx(3,0)
1007  mx(1,3) = mx(3,1)
1008  mx(2,3) = mx(3,2)
1009  mx(3,3) = 2d0*q30
1010  mx(4,3) = q30+q40-q43
1011  mx(5,3) = q30+q50-q53
1012  mx(0,4) = mx(4,0)
1013  mx(1,4) = mx(4,1)
1014  mx(2,4) = mx(4,2)
1015  mx(3,4) = mx(4,3)
1016  mx(4,4) = 2d0*q40
1017  mx(5,4) = q40+q50-q54
1018  mx(0,5) = mx(5,0)
1019  mx(1,5) = mx(5,1)
1020  mx(2,5) = mx(5,2)
1021  mx(3,5) = mx(5,3)
1022  mx(4,5) = mx(5,4)
1023  mx(5,5) = 2d0*q50
1024 
1025  call chinv(6,mx,mxinv)
1026 
1027  ! determine X_(0,5)
1028  do j=1,5
1029  do i=1,5
1030  mx0k(i,j) = mx(i,j-1)
1031  end do
1032  end do
1033 
1034  det = chdet(5,mx0k)
1035  kbest = 5
1036 
1037  do j=5,2,-1
1038  do i=1,5
1039  mx0k(i,j) = mx(i,j)
1040  end do
1041 
1042  newdet = chdet(5,mx0k)
1043  if (abs(newdet).gt.abs(det)) then
1044  kbest = j-1
1045  det = newdet
1046  end if
1047 
1048  end do
1049 
1050  do i=1,5
1051  mx0k(i,1) = mx(i,1)
1052  mx0k(i,kbest) = mx(i,0)
1053  end do
1054 
1055  call chinv(5,mx0k,mx0kinv)
1056  do i=1,5
1057  mx0kinv(kbest,i) = 0d0
1058  end do
1059 
1060  mxinvs = sum(mxinv(0:5,0))
1061 
1062  ! determine binaries for F-coefficients
1063  k=0
1064  bin = 1
1065  do while (k.le.5)
1066  if (mod(id/bin,2).eq.0) then
1067  nid(k) = id+bin
1068  k = k+1
1069  end if
1070  bin = 2*bin
1071  end do
1072 
1073 
1074  rmax_m1 = max(rmax-1,0)
1075  bino_0 = binomtable(rmax_m1,n+rmax_m1-2)
1076  bino_i = binomtable(rmax_m1,n+rmax_m1-3)
1077  allocate(tnm1(0:rmax_m1,0:rmax_m1,0:rmax_m1,0:rmax_m1,0:5))
1078  allocate(tnm1_0aux(0:rmax_m1,0:rmax_m1,0:rmax_m1,0:rmax_m1))
1079  allocate(tnm1uv_0aux(0:rmax_m1,0:rmax_m1,0:rmax_m1,0:rmax_m1))
1080  allocate(tnm1err(0:5,0:rmax_m1))
1081 
1082  if(rmax_m1.ge.2*n-6) then
1083  allocate(cnm1uv_0(bino_0,n-3:rmax/2,2*n-6:rmax))
1084  allocate(cnm1uv_0aux(bino_i,n-3:rmax_m1/2,2*n-6:rmax_m1))
1085  allocate(cnm1uv_i(bino_i,0:rmax_m1/2,0:rmax_m1))
1086  call calctensortnr(tnm1_0aux(:,:,:,:),tnm1uv_0aux,tnm1err(0,0:rmax_m1),submomvec(n,0,momvec),submominv(n,0,mominv), &
1087  submasses(n,0,masses2),n-1,rmax_m1,nid(0),cnm1uv_0aux)
1088  do i=1,5
1089  call calctensortnr(tnm1(:,:,:,:,i),tnm1uv_0aux,tnm1err(i,0:rmax_m1),submomvec(n,i,momvec),submominv(n,i,mominv), &
1090  submasses(n,i,masses2),n-1,rmax_m1,nid(i),cnm1uv_i)
1091  end do
1092  else
1093  call calctensortnr(tnm1_0aux(:,:,:,:),tnm1uv_0aux,tnm1err(0,0:rmax_m1),submomvec(n,0,momvec), &
1094  submominv(n,0,mominv),submasses(n,0,masses2),n-1,rmax_m1,nid(0))
1095  do i=1,5
1096  call calctensortnr(tnm1(:,:,:,:,i),tnm1uv_0aux,tnm1err(i,0:rmax_m1),submomvec(n,i,momvec), &
1097  submominv(n,i,mominv),submasses(n,i,masses2),n-1,rmax_m1,nid(i))
1098  end do
1099  end if
1100 
1101  ! shift of integration momentum in TE\{0}
1102  tnm1(:,:,:,:,0) = 0d0
1103 
1104  do r=0,rmax_m1
1105 
1106  do np0=0,r
1107  do np1=0,r-np0
1108  do np2=0,r-np0-np1
1109  do np3=0,r-np0-np1-np2
1110 
1111  p1ten = (-momvec(0,1))**np0*(-momvec(1,1))**np1* &
1112  (-momvec(2,1))**np2*(-momvec(3,1))**np3
1113 
1114  do n0=0,r-np0-np1-np2-np3
1115  do n1=0,r-np0-np1-np2-np3-n0
1116  do n2=0,r-np0-np1-np2-np3-n0-n1
1117  n3 = r-np0-np1-np2-np3-n0-n1-n2
1118 
1119  combi = binomtable(n0,n0+np0)*binomtable(n1,n1+np1)* &
1120  binomtable(n2,n2+np2)*binomtable(n3,n3+np3)
1121 
1122  tnm1(np0+n0,np1+n1,np2+n2,np3+n3,0) = tnm1(np0+n0,np1+n1,np2+n2,np3+n3,0) &
1123  + combi*tnm1_0aux(n0,n1,n2,n3)*p1ten
1124 
1125  end do
1126  end do
1127  end do
1128 
1129  end do
1130  end do
1131  end do
1132  end do
1133 
1134  end do
1135 
1136  p1max = maxval(abs(momvec(0:3,1)))
1137  do r=1,rmax_m1
1138  do i=1,r
1139  tnm1err(0,r) = max(tnm1err(0,i),tnm1err(0,r-i)*p1max**i)
1140  end do
1141  end do
1142 
1143  tn = 0d0
1144  tn(0,0,0,0) = -mxinv(0,0)*tnm1(0,0,0,0,0)
1145  do i=1,5
1146  tn(0,0,0,0) = tn(0,0,0,0) + mxinv(i,0)*(tnm1(0,0,0,0,i)-tnm1(0,0,0,0,0))
1147  end do
1148  tnerr(0) = max( abs(mxinvs)*tnm1err(0,0), &
1149  abs(mxinv(1,0))*tnm1err(1,0) , &
1150  abs(mxinv(2,0))*tnm1err(2,0) , &
1151  abs(mxinv(3,0))*tnm1err(3,0) , &
1152  abs(mxinv(4,0))*tnm1err(4,0) , &
1153  abs(mxinv(5,0))*tnm1err(5,0) )
1154 
1155  do mu=0,3
1156  do i=1,5
1157  mxinvp(i,mu)=0d0
1158  do j=1,5
1159  mxinvp(i,mu) = mxinvp(i,mu) + mx0kinv(j,i)*momvec(mu,j)
1160  end do
1161  end do
1162  mxinvps(mu) = sum(mxinvp(1:5,mu))
1163  end do
1164 
1165  ! calculation of UV-divergent parts
1166  ! TN UV-finite for n0<N-2
1167  ! TNuv(:,0:min(rmax/2,N-3),:) = 0d0
1168 
1169  ! PV reduction (5.10)
1170  if(rmax.ge.2*n-4) then
1171  cnm1uv_0 = 0d0
1172  do n0=n-3,rmax_m1/2
1173  cnm1uv_0(1,n0,2*n0) = cnm1uv_0aux(1,n0,2*n0)
1174  do r=1,rmax_m1-2*n0
1175  mia1 = binomtable(r-1,n+r-3)+1
1176  mia2 = binomtable(r,n+r-3)
1177  mia3 = binomtable(r,n+r-2)
1178  cnm1uv_0(mia1:mia3,n0,r+2*n0) = cnm1uv_0aux(1:mia2,n0,r+2*n0)
1179  end do
1180  end do
1181  do r=2*n-6,rmax_m1
1182  do n0=n-3,(r-1)/2
1183  do i=binomtable(r-2*n0-1,n+r-2*n0-3),1,-1
1184  cnm1uv_0(i,n0,r) = -cnm1uv_0(i,n0,r-1)
1185  do j=2,n-1
1186  cnm1uv_0(i,n0,r) = cnm1uv_0(i,n0,r) &
1187  - cnm1uv_0(addtocind(j,i,r-2*n0-1,n-1),n0,r)
1188  end do
1189  end do
1190  end do
1191  end do
1192  cnuv = 0d0
1193  do r=2*n-4,rmax
1194  do n0=n-2,r/2
1195  do cind=1,binomtable(r-2*n0,n+r-2*n0-2)
1196 
1197  cnuv(cind,n0,r) = cnm1uv_0(cind,n0-1,r-2)
1198  if(n0.ge.n-1) then
1199  cnuv(cind,n0,r) = cnuv(cind,n0,r) + 2*mm02*cnuv(cind,n0-1,r-2)
1200  do i=1,n-1
1201  cnuv(cind,n0,r) = cnuv(cind,n0,r) &
1202  + ff(i)*cnuv(addtocind(i,cind,r-2*n0,n-1),n0-1,r-1)
1203  end do
1204  end if
1205  cnuv(cind,n0,r) = cnuv(cind,n0,r)/(2*(r+3-n))
1206 
1207  end do
1208  end do
1209  end do
1210  end if
1211 
1212  do r=1,rmax
1213  do n0=0,r
1214  do n1=0,r-n0
1215  do n2=0,r-n0-n1
1216  n3 = r-n0-n1-n2
1217 
1218  if (n0.ge.1) then
1219  do i=1,5
1220  smod = tnm1(n0-1,n1,n2,n3,i) - tnm1(n0-1,n1,n2,n3,0)
1221  tn(n0,n1,n2,n3) = tn(n0,n1,n2,n3) + mxinvp(i,0)*smod
1222  end do
1223  else if (n1.ge.1) then
1224  do i=1,5
1225  smod = tnm1(n0,n1-1,n2,n3,i) - tnm1(n0,n1-1,n2,n3,0)
1226  tn(n0,n1,n2,n3) = tn(n0,n1,n2,n3) + mxinvp(i,1)*smod
1227  end do
1228  else if (n2.ge.1) then
1229  do i=1,5
1230  smod = tnm1(n0,n1,n2-1,n3,i) - tnm1(n0,n1,n2-1,n3,0)
1231  tn(n0,n1,n2,n3) = tn(n0,n1,n2,n3) + mxinvp(i,2)*smod
1232  end do
1233  else
1234  do i=1,5
1235  smod = tnm1(n0,n1,n2,n3-1,i) - tnm1(n0,n1,n2,n3-1,0)
1236  tn(n0,n1,n2,n3) = tn(n0,n1,n2,n3) + mxinvp(i,3)*smod
1237  end do
1238  end if
1239 
1240  end do
1241  end do
1242  end do
1243 
1244  tnerr(r) = max(maxval(abs(mxinvps(0:3)))*tnm1err(0,r-1), &
1245  maxval(abs(mxinvp(1,0:3)))*tnm1err(1,r-1), &
1246  maxval(abs(mxinvp(2,0:3)))*tnm1err(2,r-1), &
1247  maxval(abs(mxinvp(3,0:3)))*tnm1err(3,r-1))
1248 
1249  end do
1250 
1251 
1252  end subroutine calctensortnrred
1253 
1254 
1255 
1256 
1257 
1258  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1259  ! subroutine CalcTensorTNr_list(TN,TNuv,TNerr,MomVec,MomInv,masses2,rmax)
1260  !
1261  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1262 
1263  subroutine calctensortnr_list(TN,TNuv,TNerr,MomVec,MomInv,masses2,N,rmax)
1264 
1265  integer, intent(in) :: N,rmax
1266  double complex, intent(in) :: MomVec(0:3,N-1), masses2(0:N-1)
1267  double complex, intent(in) :: MomInv(BinomTable(2,N))
1268  double complex, intent(out) :: TN(RtS(rmax)),TNuv(RtS(rmax))
1269  double precision, intent(out) :: TNerr(0:rmax)
1270  double complex :: TN_aux(0:rmax,0:rmax,0:rmax,0:rmax)
1271  double complex :: TNuv_aux(0:rmax,0:rmax,0:rmax,0:rmax)
1272  integer :: mu
1273 
1274  call calctensortnr(tn_aux,tnuv_aux,tnerr,momvec,mominv,masses2,n,rmax,0)
1275 
1276  do mu=1,rts(rmax)
1277  tn(mu) = tn_aux(lorindtab(0,mu),lorindtab(1,mu),lorindtab(2,mu),lorindtab(3,mu))
1278  end do
1279 
1280  if (calcuv_cll) then
1281  do mu=1,rts(rmax)
1282  tnuv(mu) = tnuv_aux(lorindtab(0,mu),lorindtab(1,mu),lorindtab(2,mu),lorindtab(3,mu))
1283  end do
1284  end if
1285 
1286  end subroutine calctensortnr_list
1287 
1288 
1289 
1290 
1291 
1292  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1293  ! subroutine CheckTensors(TN,MomInv,N,rmax)
1294  !
1295  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1296 
1297  subroutine checktensors(TN,MomVec,MomInv,masses2,N,rmax,writeout,acc)
1298 
1299  integer, intent(in) :: N,rmax,writeout
1300  double precision, intent(in) :: acc
1301  double complex, intent(in) :: TN(0:rmax,0:rmax,0:rmax,0:rmax)
1302  double complex, intent(in) :: MomVec(0:3,N-1), MomInv(BinomTable(2,N)), masses2(0:N-1)
1303  double complex :: CoefsN(BinomTable(rmax,N+rmax-2),0:rmax/2,0:rmax)
1304  double complex :: CoefsNuv(BinomTable(rmax,N+rmax-2),0:rmax/2,0:rmax)
1305  double complex :: check1, check2, fac, Gram(N-1,N-1)
1306  integer, allocatable :: pinds(:),loinds(:)
1307  integer :: struc1(0:N-1), struc2(0:N-1)
1308  integer :: r,n0,i,j,mu,mu0,sgn,tinds(0:3),cnt,nn0,ii
1309  double precision :: TNerr(0:rmax),TNerr2(0:rmax)
1310 
1311  call calctn(coefsn,coefsnuv,mominv,masses2,n,rmax,0,tnerr,tnerr2)
1312 
1313 
1314  ! determine Gram matrix
1315  cnt = 1
1316  do i=1,n-1
1317  gram(i,i) = mominv(cnt)
1318  cnt = cnt+n-i
1319  end do
1320 
1321  cnt = 2
1322  do i=1,n-1
1323  do j=i+1,n-1
1324  gram(i,j) = (gram(i,i)+gram(j,j)-mominv(cnt))/2d0
1325  gram(j,i) = gram(i,j)
1326  cnt = cnt+1
1327  end do
1328  cnt = cnt+1
1329  end do
1330 
1331 
1332  do r=1,rmax
1333 
1334  if (allocated(loinds)) then
1335  deallocate(loinds)
1336  end if
1337  allocate(loinds(r))
1338 
1339  do n0=r/2,0,-1
1340 
1341  if (allocated(pinds)) then
1342  deallocate(pinds)
1343  end if
1344  allocate(pinds(r-2*n0))
1345 
1346  do i=1,binomtable(r-2*n0,r-2*n0+n-2)
1347  if (r.gt.2*n0) then
1348  pinds = indcombiseq(1:r-2*n0,i,r-2*n0,n-1)
1349  end if
1350 
1351  struc1(1:n-1) = calccindarr(n-1,r-2*n0,i)
1352  ! contract TN(r) with g(n0)*MomVec(pinds)
1353  ! --> check1
1354  check1 = 0d0
1355 
1356  loinds = 0
1357  do mu=0,4**r-1
1358 
1359  mu0 = mu
1360  tinds = 0
1361  do j=1,r
1362  loinds(j) = mod(mu0,4)
1363  tinds(loinds(j)) = tinds(loinds(j))+1
1364  mu0 = mu0/4
1365  end do
1366 
1367  ! sign of contraction from Minkowski metric
1368  sgn=(-1)**(r-tinds(0))
1369  do j=1,n0
1370  select case (loinds(2*j)-loinds(2*j-1))
1371  case (0)
1372  if (loinds(2*j).ne.0) then
1373  sgn = -sgn
1374  end if
1375  case default
1376  sgn = 0
1377  end select
1378  end do
1379 
1380  if (sgn.ne.0) then
1381 
1382  fac = 1d0
1383  do j=1,r-2*n0
1384  fac = fac*momvec(loinds(2*n0+j),pinds(j))
1385  end do
1386 
1387  check1 = check1 + sgn*fac*tn(tinds(0),tinds(1),tinds(2),tinds(3))
1388 
1389  end if
1390 
1391  end do
1392 
1393 
1394  ! calculate contraction directly from coefficients
1395  ! --> check2
1396  check2 = 0d0
1397 
1398  struc1(0) = n0
1399  struc1(1:n-1) = calccindarr(n-1,r-2*n0,i)
1400 
1401  do nn0=r/2,0,-1
1402  do ii=1,binomtable(r-2*nn0,r-2*nn0+n-2)
1403 
1404  struc2(0) = nn0
1405  struc2(1:n-1) = calccindarr(n-1,r-2*nn0,ii)
1406 
1407  check2 = check2 + coefsn(ii,nn0,r)* &
1408  contractlostruc(n-1,struc1,struc2,gram)
1409 
1410  end do
1411  end do
1412 
1413  select case (writeout)
1414 
1415  case(1)
1416  if (abs(check1-check2)/abs(check2).gt.acc) then
1417  write(ninfout_cll,*) struc1
1418  write(ninfout_cll,*) 'check1:', check1
1419  write(ninfout_cll,*) 'check2:', check2
1420  end if
1421 
1422  case(2)
1423  write(ninfout_cll,*) struc1
1424  write(ninfout_cll,*) 'check1:', check1
1425  write(ninfout_cll,*) 'check2:', check2
1426 
1427  end select
1428 
1429  end do
1430  end do
1431  end do
1432 
1433  end subroutine checktensors
1434 
1435 
1436 
1437 
1438 
1439 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1440 ! ! function ContractLoStruc(Nm1,struc1,struc2,Gram)
1441 ! !
1442 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1443 !
1444 ! recursive function ContractLoStruc(Nm1,struc1,struc2,Gram) result(res)
1445 !
1446 ! integer, intent(in) :: Nm1
1447 ! integer, intent(in) :: struc1(0:Nm1), Struc2(0:Nm1)
1448 ! integer :: struc1aux(0:Nm1), struc2aux(0:Nm1), struc2aux2(0:Nm1)
1449 ! double complex, intent(in) :: Gram(Nm1,Nm1)
1450 ! double complex :: res
1451 ! integer :: i,j,con,sum1,sum2,fac
1452 !
1453 ! sum1 = 2*struc1(0)
1454 ! sum2 = 2*struc2(0)
1455 ! do i=1,Nm1
1456 ! sum1 = sum1 + struc1(i)
1457 ! sum2 = sum2 + struc2(i)
1458 ! end do
1459 !
1460 ! if (sum1.ne.sum2) then
1461 ! write(*,*) 'error in ContractLoStruc:'
1462 ! write(*,*) 'struc1 and struc2 must be of equal rank!'
1463 ! stop
1464 ! end if
1465 !
1466 ! if (sum1.eq.0) then
1467 ! res = 1d0
1468 ! return
1469 ! end if
1470 !
1471 !
1472 ! do i=0,Nm1
1473 ! if (struc1(i).ge.1) then
1474 ! con = i
1475 ! end if
1476 ! end do
1477 !
1478 ! res = 0d0
1479 ! struc1aux = struc1
1480 ! struc1aux(con) = struc1aux(con)-1
1481 !
1482 ! if (con.ge.1) then
1483 !
1484 ! ! contract p_con from T1 with g from T2
1485 ! if (struc2(0).ge.1) then
1486 ! struc2aux = struc2
1487 ! struc2aux(0) = struc2aux(0)-1
1488 ! struc2aux(con) = struc2aux(con)+1
1489 ! ! go on contracting recursively
1490 ! ! (factor struc2aux(con) because of symmetrization wrt. g and pi)
1491 ! res = res + struc2aux(con)*ContractLoStruc(Nm1,struc1aux,struc2aux,Gram)
1492 ! end if
1493 !
1494 ! ! contract p_con from T1 with all the pi from T2
1495 ! do i=1,Nm1
1496 ! if (struc2(i).ge.1) then
1497 ! struc2aux = struc2
1498 ! struc2aux(i) = struc2aux(i)-1
1499 ! ! go on contracting recursively
1500 ! res = res + Gram(i,con)*ContractLoStruc(Nm1,struc1aux,struc2aux,Gram)
1501 ! end if
1502 ! end do
1503 !
1504 !
1505 ! else
1506 !
1507 ! ! contract g from T1 with g from T2
1508 ! if (struc2(0).ge.1) then
1509 ! struc2aux = struc2
1510 ! struc2aux(0) = struc2aux(0)-1
1511 ! ! full contraction g^{mu,nu}.g_{mu,nu}
1512 ! ! tensor in D=4 dimensions g*g=4
1513 ! fac = 4
1514 ! do i=0,Nm1
1515 ! ! partial contration g^{mu,nu}.(pi_mu g_{nu,rho})
1516 ! ! or g^{mu,nu}.(g_{mu,rho}g_{nu,sigma})
1517 ! ! factor 2 for mu <--> nu
1518 ! fac = fac + 2*struc2aux(i)
1519 ! end do
1520 ! ! go on contracting recursively
1521 ! res = res + fac*ContractLoStruc(Nm1,struc1aux,struc2aux,Gram)
1522 ! end if
1523 !
1524 ! ! contract g^{mu,nu} from T1 with pi,pj from T2
1525 ! do i=1,Nm1
1526 ! if (struc2(i).ge.1) then
1527 ! struc2aux = struc2
1528 ! struc2aux(i) = struc2aux(i)-1
1529 ! do j=1,Nm1
1530 ! if (struc2aux(j).ge.1) then
1531 ! struc2aux2 = struc2aux
1532 ! struc2aux2(j) = struc2aux2(j)-1
1533 ! ! go on contracting recursively
1534 ! res = res + Gram(i,j)*ContractLoStruc(Nm1,struc1aux,struc2aux2,Gram)
1535 ! end if
1536 ! end do
1537 ! end if
1538 ! end do
1539 !
1540 ! end if
1541 !
1542 !
1543 ! end function ContractLoStruc
1544 
1545 
1546 
1547 
1548 
1549  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1550  ! subroutine CompareTensors
1551  !
1552  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1553 
1554  subroutine comparetensors(T1,T2,rmax,writeout,acc)
1555 
1556  integer, intent(in) :: rmax, writeout
1557  double complex, intent(in) :: T1(0:rmax,0:rmax,0:rmax,0:rmax)
1558  double complex, intent(in) :: T2(0:rmax,0:rmax,0:rmax,0:rmax)
1559  double precision, intent(in) :: acc
1560  integer :: r,n0,n1,n2,n3
1561 
1562  do r=0,rmax
1563  do n0=0,r
1564  do n1=0,r-n0
1565  do n2=0,r-n0-n1
1566  n3 = r-n0-n1-n2
1567 
1568  if (writeout.eq.1) then
1569  if (abs(t1(n0,n1,n2,n3)-t2(n0,n1,n2,n3)).gt. &
1570  acc*abs(t1(n0,n1,n2,n3)+t2(n0,n1,n2,n3))) then
1571  write(ninfout_cll,*) n0,n1,n2,n3
1572  write(ninfout_cll,*) t1(n0,n1,n2,n3), t2(n0,n1,n2,n3)
1573  end if
1574  end if
1575 
1576  end do
1577  end do
1578  end do
1579  end do
1580 
1581  end subroutine comparetensors
1582 
1583 
1584 end module tensorreduction
1585 
submasses
double complex function, dimension(0:n-2) submasses(N, k, masses)
Definition: democache.f90:468
tensorreduction::checktensors
subroutine checktensors(TN, MomVec, MomInv, masses2, N, rmax, writeout, acc)
Definition: TensorReduction.F90:1298
collier_coefs::e_main_cll
subroutine e_main_cll(E, Euv, p10, p21, p32, p43, p40, p20, p31, p42, p30, p41, m02, m12, m22, m32, m42, rmax, Eerr, id_in, Eerr2)
Definition: collier_coefs.F90:1780
collier_coefs
Definition: collier_coefs.F90:28
tensorreduction::comparetensors
subroutine comparetensors(T1, T2, rmax, writeout, acc)
Definition: TensorReduction.F90:1555
buildtensors::calctensortnuv_list
subroutine calctensortnuv_list(TNuv, CoefsNuv, MomVec, N, rmax)
Definition: BuildTensors.F90:1586
submominv
double complex function, dimension((n-1) *(n-2)/2) submominv(N, k, MomInv)
Definition: democache.f90:495
buildtensors
Definition: BuildTensors.F90:25
buildtensors::calctensorfuv
subroutine calctensorfuv(TFuv, CoefsFuv, MomVec, rmax)
Definition: BuildTensors.F90:1153
buildtensors::calctensore
subroutine calctensore(TE, TEuv, TEerr, CoefsE, CoefsEuv, CoefsEerr, MomVec, rmax)
Definition: BuildTensors.F90:804
collier_coefs::tn_cll
Definition: collier_coefs.F90:76
tensorreduction::calctensorfrred
subroutine calctensorfrred(TF, CFuv, TFerr, MomVec, MomInv, masses2, rmax)
Definition: TensorReduction.F90:260
tensorreduction::calctensortnr
recursive subroutine calctensortnr(TN, TNuv, TNerr, MomVec, MomInv, masses2, N, rmax, id, CNuv)
Definition: TensorReduction.F90:617
tensorreduction::calctensortnrred
recursive subroutine calctensortnrred(TN, TNerr, MomVec, MomInv, masses2, N, rmax, id, CNuv)
Definition: TensorReduction.F90:927
tensorreduction
Definition: TensorReduction.F90:25
tensorreduction::calctensorfr_list
subroutine calctensorfr_list(TF, TFuv, TFerr, MomVec, MomInv, masses2, rmax)
Definition: TensorReduction.F90:584
tensorreduction::calctensorfr
subroutine calctensorfr(TF, TFuv, TFerr, MomVec, MomInv, masses2, rmax)
Definition: TensorReduction.F90:45
buildtensors::calctensortn
subroutine calctensortn(TN, TNuv, TNerr, CoefsN, CoefsNuv, CoefsNerr, MomVec, N, rmax)
Definition: BuildTensors.F90:1546
tensorreduction::calctensortnr_list
subroutine calctensortnr_list(TN, TNuv, TNerr, MomVec, MomInv, masses2, N, rmax)
Definition: TensorReduction.F90:1264