JHUGen MELA  v2.4.1
Matrix element calculations as used in JHUGen. MELA is an important tool that was used for the Higgs boson discovery and for precise measurements of its structure and interactions. Please see the website https://spin.pha.jhu.edu/ and papers cited there for more details, and kindly cite those papers when using this code.
reductionEFG.F90
Go to the documentation of this file.
1 !!
2 !! File reductionEFG.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 reductionEFG *
14 ! * by Lars Hofer *
15 ! *************************
16 !
17 ! functions and subroutines:
18 ! CalcE, CalcF, CalcG
19 !
20 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
21 
22 
23 
24 
25 module reductionefg
26 
27  use reductiond
28 
29  implicit none
30 
31 
32 
33 contains
34 
35 
36  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
37  ! subroutine CalcE(E,Euv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,Eerr)
38  !
39  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
40 
41  subroutine calce(E,Euv,p10,p21,p32,p43,p40,p20,p31,p42,p30,p41, &
42  m02,m12,m22,m32,m42,rmax,id,Eerr,Eerr2)
43 
44  integer, intent(in) :: rmax,id
45  double complex, intent(in) :: p10,p21,p32,p43,p40,p20,p31,p42,p30,p41
46  double complex, intent(in) :: m02,m12,m22,m32,m42
47  double complex, intent(out) :: E(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
48  double complex, intent(out) :: Euv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
49  double precision :: Eerr(0:rmax),Eerr2(0:rmax)
50  double complex, allocatable :: Eaux(:,:,:,:,:), Euvaux(:,:,:,:,:), fct(:)
51  double precision, allocatable :: Eerraux(:),Eerr2aux(:)
52  double complex :: x(15)
53  integer :: rank,switch,cnt,n0,n1,n2,n3,n4,r
54  logical :: nocalc,wrica
55 
56 
57  if ((use_cache_system).and.(tencache.gt.5)) then
58  if ((ncache.gt.0).and.(ncache.le.ncache_max)) then
59 ! if (use_cache(ncache).ge.5) then
60  x(1)=p10
61  x(2)=p21
62  x(3)=p32
63  x(4)=p43
64  x(5)=p40
65  x(6)=p20
66  x(7)=p31
67  x(8)=p42
68  x(9)=p30
69  x(10)=p41
70  x(11)=m02
71  x(12)=m12
72  x(13)=m22
73  x(14)=m32
74  x(15)=m42
75 
76  rank = rmax
77 
78  if (rmax.ge.6) then
79  allocate(fct(ncoefs(rmax,5)+ncoefs(rmax-6,5)+2*(rmax+1)))
80  call readcache(fct,ncoefs(rmax,5)+ncoefs(rmax-6,5)+2*(rmax+1),x,15,1,id,5,rank,nocalc,wrica)
81  else
82  allocate(fct(ncoefs(rmax,5)+2*(rmax+1)))
83  call readcache(fct,ncoefs(rmax,5)+2*(rmax+1),x,15,1,id,5,rank,nocalc,wrica)
84  end if
85 
86  if(nocalc)then
87  cnt = 0
88  euv(0:min(rmax/2,2),:,:,:,:) = 0d0
89  do r=0,rmax
90  do n0=0,r/2
91  do n1=0,r-2*n0
92  do n2=0,r-2*n0-n1
93  do n3=0,r-2*n0-n1-n2
94  n4 = r-2*n0-n1-n2-n3
95 
96  cnt = cnt+1
97  e(n0,n1,n2,n3,n4) = fct(cnt)
98 
99  end do
100  end do
101  end do
102  end do
103  do n0=3,r/2
104  do n1=0,r-2*n0
105  do n2=0,r-2*n0-n1
106  do n3=0,r-2*n0-n1-n2
107  n4 = r-2*n0-n1-n2-n3
108 
109  cnt = cnt+1
110  euv(n0,n1,n2,n3,n4) = fct(cnt)
111 
112  end do
113  end do
114  end do
115  end do
116  cnt = cnt+1
117  eerr(r) = real(fct(cnt))
118  cnt = cnt+1
119  eerr2(r) = real(fct(cnt))
120  end do
121  return
122  endif
123 
124 
125  if(rank.eq.rmax) then
126 
127  call calcered(e,euv,p10,p21,p32,p43,p40,p20,p31,p42,p30,p41, &
128  m02,m12,m22,m32,m42,rank,id,eerr,eerr2)
129 
130  if (wrica) then
131  cnt = 0
132  do r=0,rank
133  do n0=0,r/2
134  do n1=0,r-2*n0
135  do n2=0,r-2*n0-n1
136  do n3=0,r-2*n0-n1-n2
137  n4 = r-2*n0-n1-n2-n3
138 
139  cnt = cnt+1
140  fct(cnt) = e(n0,n1,n2,n3,n4)
141 
142  end do
143  end do
144  end do
145  end do
146  do n0=3,r/2
147  do n1=0,r-2*n0
148  do n2=0,r-2*n0-n1
149  do n3=0,r-2*n0-n1-n2
150  n4 = r-2*n0-n1-n2-n3
151 
152  cnt = cnt+1
153  fct(cnt) = euv(n0,n1,n2,n3,n4)
154 
155  end do
156  end do
157  end do
158  end do
159  cnt = cnt+1
160  fct(cnt) = eerr(r)
161  cnt = cnt+1
162  fct(cnt) = eerr2(r)
163  end do
164 
165  if (rank.ge.6) then
166  call writecache(fct,ncoefs(rank,5)+ncoefs(rank-6,5)+2*(rank+1),id,5,rank)
167  else
168  call writecache(fct,ncoefs(rank,5)+2*(rank+1),id,5,rank)
169  end if
170 
171  end if
172 
173  return
174 
175 
176  else
177  allocate(eaux(0:rank/2,0:rank,0:rank,0:rank,0:rank))
178  allocate(euvaux(0:rank/2,0:rank,0:rank,0:rank,0:rank))
179  allocate(eerraux(0:rank))
180  allocate(eerr2aux(0:rank))
181 
182  call calcered(eaux,euvaux,p10,p21,p32,p43,p40,p20,p31,p42,p30,p41, &
183  m02,m12,m22,m32,m42,rank,id,eerraux,eerr2aux)
184 
185  if (wrica) then
186  cnt = 0
187  if (rank.ge.6) then
188  deallocate(fct)
189  allocate(fct(ncoefs(rank,5)+ncoefs(rank-6,5)+2*(rank+1)))
190  else
191  deallocate(fct)
192  allocate(fct(ncoefs(rank,5)+2*(rank+1)))
193  end if
194  do r=0,rank
195  do n0=0,r/2
196  do n1=0,r-2*n0
197  do n2=0,r-2*n0-n1
198  do n3=0,r-2*n0-n1-n2
199  n4 = r-2*n0-n1-n2-n3
200 
201  cnt = cnt+1
202  fct(cnt) = eaux(n0,n1,n2,n3,n4)
203 
204  end do
205  end do
206  end do
207  end do
208  do n0=3,r/2
209  do n1=0,r-2*n0
210  do n2=0,r-2*n0-n1
211  do n3=0,r-2*n0-n1-n2
212  n4 = r-2*n0-n1-n2-n3
213 
214  cnt = cnt+1
215  fct(cnt) = euvaux(n0,n1,n2,n3,n4)
216 
217  end do
218  end do
219  end do
220  end do
221  cnt = cnt+1
222  fct(cnt) = eerraux(r)
223  cnt = cnt+1
224  fct(cnt) = eerr2aux(r)
225  end do
226 
227  if (rank.ge.6) then
228  call writecache(fct,ncoefs(rank,5)+ncoefs(rank-6,5)+2*(rank+1),id,5,rank)
229  else
230  call writecache(fct,ncoefs(rank,5)+2*(rank+1),id,5,rank)
231  end if
232 
233  end if
234 
235  e = eaux(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
236  euv = euvaux(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
237  eerr = eerraux(0:rmax)
238  eerr2 = eerr2aux(0:rmax)
239  return
240 
241  end if
242 
243 ! end if
244  end if
245  end if
246 
247 
248  call calcered(e,euv,p10,p21,p32,p43,p40,p20,p31,p42,p30,p41, &
249  m02,m12,m22,m32,m42,rmax,id,eerr,eerr2)
250 
251 
252  end subroutine calce
253 
254 
255 
256 
257 
258  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
259  ! subroutine CalcEred(E,Euv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,Eerr,Eerr2)
260  !
261  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
262 
263  subroutine calcered(E,Euv,p10,p21,p32,p43,p40,p20,p31,p42,p30,p41, &
264  m02,m12,m22,m32,m42,rmax,id,Eerr,Eerr2)
265 
266  integer, intent(in) :: rmax,id
267  double complex, intent(in) :: p10,p21,p32,p43,p40,p20,p31,p42,p30,p41
268  double complex, intent(in) :: m02,m12,m22,m32,m42
269  double precision, intent(out) :: Eerr(0:rmax),Eerr2(0:rmax)
270  double complex :: q10,q21,q32,q43,q40,q20,q31,q42,q30,q41
271  double complex :: mm02,mm12,mm22,mm32,mm42
272  double complex :: mx(0:4,0:4), mxinv(0:4,0:4), f(4),mxinvs(0:4)
273  double complex :: zmxinv(4,0:4),Z(4,4),zmxinvs(4)
274  double precision :: maxZ,maxzmxinv(0:4),maxzmxinvs
275  double complex, intent(out) :: E(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
276  double complex, intent(out) :: Euv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
277  double complex :: Euv_aux(0:(rmax+1)/2,0:(rmax+1),0:(rmax+1),0:(rmax+1),0:(rmax+1))
278  double complex, allocatable :: D_0(:,:,:,:,:), Duv_0(:,:,:,:,:)
279  double complex, allocatable :: D_i(:,:,:,:,:), Duv_i(:,:,:,:,:)
280 ! double complex :: D_0(0:rmaxD,0:rmaxD,0:rmaxD,0:rmaxD,0:rmaxD)
281 ! double complex :: Duv_0(0:rmaxD,0:rmaxD,0:rmaxD,0:rmaxD,0:rmaxD)
282 ! double complex :: D_i(0:rmaxD,0:rmaxD,0:rmaxD,0:rmaxD,4)
283 ! double complex :: Duv_i(0:rmaxD,0:rmaxD,0:rmaxD,0:rmaxD,4)
284  double complex :: S(0:4), S2(4,4), Eaux(4), Eaux2, Eaux3, elimminf2_coli
285  double precision, allocatable :: Derr(:,:),Derr2(:,:)
286 ! double precision :: Derr(0:4,0:rmaxD)
287  integer :: r,n0,n1,n2,n3,n4,i,n,k,nid(0:4),rid,r0,bin,rBCD,rmaxD
288 
289 ! write(*,*) 'CalcEred in',p10,p21,p32,p43,p40,p20,p31,p42,p30,p41, &
290 ! m02,m12,m22,m32,m42,rmax,id
291 
292  rmaxd = max(rmax-1,0)
293 ! rmaxD = rmax
294  allocate(d_0(0:rmaxd,0:rmaxd,0:rmaxd,0:rmaxd,0:rmaxd))
295  allocate(duv_0(0:rmaxd,0:rmaxd,0:rmaxd,0:rmaxd,0:rmaxd))
296  allocate(d_i(0:rmaxd,0:rmaxd,0:rmaxd,0:rmaxd,4))
297  allocate(duv_i(0:rmaxd,0:rmaxd,0:rmaxd,0:rmaxd,4))
298  allocate(derr(0:4,0:rmaxd))
299  allocate(derr2(0:4,0:rmaxd))
300 
301  r0 = 0
302  ! write(*,*) 'CalcE', id, rmax_cp(id)
303  ! determine binaries for B-coefficients
304  k=0
305  bin = 1
306  do while (k.le.4)
307  if (mod(id/bin,2).eq.0) then
308  nid(k) = id+bin
309  k = k+1
310  end if
311  bin = 2*bin
312  end do
313 
314  call calcd(d_0(:,0,:,:,:),duv_0(:,0,:,:,:),p21,p32,p43,p41,p31,p42, &
315  m12,m22,m32,m42,rmaxd,nid(0),derr(0,0:rmaxd),derr2(0,0:rmaxd))
316  call calcd(d_i(:,:,:,:,1),duv_i(:,:,:,:,1),p20,p32,p43,p40,p30,p42, &
317  m02,m22,m32,m42,rmaxd,nid(1),derr(1,0:rmaxd),derr2(1,0:rmaxd))
318  call calcd(d_i(:,:,:,:,2),duv_i(:,:,:,:,2),p10,p31,p43,p40,p30,p41, &
319  m02,m12,m32,m42,rmaxd,nid(2),derr(2,0:rmaxd),derr2(2,0:rmaxd))
320  call calcd(d_i(:,:,:,:,3),duv_i(:,:,:,:,3),p10,p21,p42,p40,p20,p41, &
321  m02,m12,m22,m42,rmaxd,nid(3),derr(3,0:rmaxd),derr2(3,0:rmaxd))
322  call calcd(d_i(:,:,:,:,4),duv_i(:,:,:,:,4),p10,p21,p32,p30,p20,p31, &
323  m02,m12,m22,m32,rmaxd,nid(4),derr(4,0:rmaxd),derr2(4,0:rmaxd))
324 
325  ! shift of integration momentum in C\{0}
326  do n1=1,rmaxd
327  do n2=0,rmaxd-n1
328  do n3=0,rmaxd-n1-n2
329  do n4=0,rmaxd-n1-n2-n3
330  n0 = (rmaxd-n1-n2-n3-n4)
331  duv_0(0:n0,n1,n2,n3,n4) = -duv_0(0:n0,n1-1,n2,n3,n4)-duv_0(0:n0,n1-1,n2+1,n3,n4) &
332  -duv_0(0:n0,n1-1,n2,n3+1,n4)-duv_0(0:n0,n1-1,n2,n3,n4+1)
333  d_0(0:n0,n1,n2,n3,n4) = -d_0(0:n0,n1-1,n2,n3,n4)-d_0(0:n0,n1-1,n2+1,n3,n4) &
334  -d_0(0:n0,n1-1,n2,n3+1,n4)-d_0(0:n0,n1-1,n2,n3,n4+1)
335  end do
336  end do
337  end do
338  end do
339 
340 
341  ! determine inverse modified Caley matrix
342  mm02 = elimminf2_coli(m02)
343  mm12 = elimminf2_coli(m12)
344  mm22 = elimminf2_coli(m22)
345  mm32 = elimminf2_coli(m32)
346  mm42 = elimminf2_coli(m42)
347  q10 = elimminf2_coli(p10)
348  q21 = elimminf2_coli(p21)
349  q32 = elimminf2_coli(p32)
350  q43 = elimminf2_coli(p43)
351  q40 = elimminf2_coli(p40)
352  q31 = elimminf2_coli(p31)
353  q42 = elimminf2_coli(p42)
354  q30 = elimminf2_coli(p30)
355  q41 = elimminf2_coli(p41)
356  q20 = elimminf2_coli(p20)
357 
358  f(1) = q10+mm02-mm12
359  f(2) = q20+mm02-mm22
360  f(3) = q30+mm02-mm32
361  f(4) = q40+mm02-mm42
362 
363  mx(0,0) = 2d0*mm02
364  mx(1,0) = q10 - mm12 + mm02
365  mx(2,0) = q20 - mm22 + mm02
366  mx(3,0) = q30 - mm32 + mm02
367  mx(4,0) = q40 - mm42 + mm02
368  mx(0,1) = mx(1,0)
369  mx(0,2) = mx(2,0)
370  mx(0,3) = mx(3,0)
371  mx(0,4) = mx(4,0)
372  mx(1,1) = 2d0*q10
373  mx(2,2) = 2d0*q20
374  mx(3,3) = 2d0*q30
375  mx(4,4) = 2d0*q40
376  mx(1,2) = q10+q20-q21
377  mx(1,3) = q10+q30-q31
378  mx(1,4) = q10+q40-q41
379  mx(2,3) = q20+q30-q32
380  mx(2,4) = q20+q40-q42
381  mx(3,4) = q30+q40-q43
382  mx(2,1) = mx(1,2)
383  mx(3,1) = mx(1,3)
384  mx(4,1) = mx(1,4)
385  mx(3,2) = mx(2,3)
386  mx(4,2) = mx(2,4)
387  mx(4,3) = mx(3,4)
388 
389  call chinv(5,mx,mxinv)
390 
391  do i=0,4
392  mxinvs(i) = sum(mxinv(i,0:4))
393  end do
394 
395  ! for alternative error estimate
396  z(1:4,1:4) = mx(1:4,1:4)
397 
398  zmxinv = matmul(z,mxinv(1:4,0:4))
399 
400  do i=0,4
401  maxzmxinv(i) = maxval(abs(zmxinv(1:4,i)))
402  end do
403 
404  do i=1,4
405  zmxinvs(i) = sum(zmxinv(i,0:4))
406  end do
407 
408  maxzmxinvs = maxval(abs(zmxinvs(1:4)))
409 
410  maxz = maxval(abs(z))
411 
412 
413  ! calculation of UV-divergent parts
414  ! E_(n0,n1,n2,n3) UV-finite for n0<3
415  euv_aux(0:min((rmax+1)/2,2),:,:,:,:) = 0d0
416  euv(0:min(rmax/2,2),:,:,:,:) = 0d0
417 
418  ! PV reduction (5.10)
419  do r=max(r0,6),rmax+1
420  do n0=3,r/2
421  do n1=0,r-2*n0
422  do n2=0,r-2*n0-n1
423  do n3=0,r-2*n0-n1-n2
424  n4 = r-2*n0-n1-n2-n3
425 
426  euv_aux(n0,n1,n2,n3,n4) = (duv_0(n0-1,n1,n2,n3,n4) &
427  + 2*m02*euv_aux(n0-1,n1,n2,n3,n4) &
428  + f(1)*euv_aux(n0-1,n1+1,n2,n3,n4) &
429  + f(2)*euv_aux(n0-1,n1,n2+1,n3,n4) &
430  + f(3)*euv_aux(n0-1,n1,n2,n3+1,n4) &
431  + f(4)*euv_aux(n0-1,n1,n2,n3,n4+1)) / (2*(r-2))
432  if (r.le.rmax) then
433  euv(n0,n1,n2,n3,n4) = euv_aux(n0,n1,n2,n3,n4)
434  end if
435 
436  end do
437  end do
438  end do
439  end do
440  end do
441 
442 ! write(*,*) 'CalcEred E0', -mxinv(0,0) -mxinv(1,0) -mxinv(2,0) -mxinv(3,0)-mxinv(4,0)
443 ! write(*,*) 'CalcEred E0', -mxinvs(0)
444 ! write(*,*) 'CalcEred E0', -mxinvs(0)+mxinv(0,0),mxinv(0,0)
445 
446  ! scalar coefficient
447  if (r0.eq.0) then
448  e = 0d0
449 ! E(0,0,0,0,0) = -mxinv(0,0)*D_0(0,0,0,0,0)
450 ! do k=1,4
451 ! E(0,0,0,0,0) = E(0,0,0,0,0) &
452 ! + mxinv(k,0)*(D_i(0,0,0,0,k)-D_0(0,0,0,0,0))
453 ! end do
454  e(0,0,0,0,0) = -mxinvs(0)*d_0(0,0,0,0,0)
455  do k=1,4
456  e(0,0,0,0,0) = e(0,0,0,0,0) + mxinv(k,0)*d_i(0,0,0,0,k)
457  end do
458  end if
459 
460 ! Eerr(0) = max( maxval(abs(mxinv(0:4,0)))*Derr(0,0), &
461  eerr(0) = max( abs(mxinvs(0))*derr(0,0), &
462  abs(mxinv(1,0))*derr(1,0) , &
463  abs(mxinv(2,0))*derr(2,0) , &
464  abs(mxinv(3,0))*derr(3,0) , &
465  abs(mxinv(4,0))*derr(4,0) )
466  eerr2(0) = max( abs(mxinvs(0))*derr2(0,0), &
467  abs(mxinv(1,0))*derr2(1,0) , &
468  abs(mxinv(2,0))*derr2(2,0) , &
469  abs(mxinv(3,0))*derr2(3,0) , &
470  abs(mxinv(4,0))*derr2(4,0) )
471 
472 ! do k=1,4
473 ! write(*,*) 'CalcEred En', -mxinv(0,k) -mxinv(1,k) -mxinv(2,k) -mxinv(3,k)-mxinv(4,k)
474 ! write(*,*) 'CalcEred En', -mxinvs(k)
475 ! end do
476 
477  ! formula (6.12) and (6.13)
478  do r=r0,rmax
479  do n0=0,r/2
480  do n1=0,r-2*n0
481  do n2=0,r-2*n0-n1
482  do n3=0,r-2*n0-n1-n2
483  n4 = r-2*n0-n1-n2-n3
484 
485  if (n0.gt.0.or.r.le.rmax-1) then
486  do n=0,4
487 ! S(n) = -D_0(n0,n1,n2,n3,n4)
488  s(n) = 0d0
489  end do
490  endif
491 
492  if (n1.eq.0) then
493  if (n0.gt.0.or.r.le.rmax-1) then
494  s(1) = s(1) + d_i(n0,n2,n3,n4,1)
495  s2(:,1) = 0d0
496  end if
497  else
498  if (r.le.rmax-1) then
499  s2(:,1) = -n1*d_0(n0+1,n1-1,n2,n3,n4)
500  if (n2.eq.0) then
501  s2(2,1) = s2(2,1) + n1*d_i(n0+1,n1-1,n3,n4,2)
502  end if
503  if (n3.eq.0) then
504  s2(3,1) = s2(3,1) + n1*d_i(n0+1,n1-1,n2,n4,3)
505  end if
506  if (n4.eq.0) then
507  s2(4,1) = s2(4,1) + n1*d_i(n0+1,n1-1,n2,n3,4)
508  end if
509  end if
510  end if
511 
512  if (n2.eq.0) then
513  if (n0.gt.0.or.r.le.rmax-1) then
514  s(2) = s(2) + d_i(n0,n1,n3,n4,2)
515  s2(:,2) = 0d0
516  end if
517  else
518  if (r.le.rmax-1) then
519  s2(:,2) = -n2*d_0(n0+1,n1,n2-1,n3,n4)
520  if (n1.eq.0) then
521  s2(1,2) = s2(1,2) + n2*d_i(n0+1,n2-1,n3,n4,1)
522  end if
523  if (n3.eq.0) then
524  s2(3,2) = s2(3,2) + n2*d_i(n0+1,n1,n2-1,n4,3)
525  end if
526  if (n4.eq.0) then
527  s2(4,2) = s2(4,2) + n2*d_i(n0+1,n1,n2-1,n3,4)
528  end if
529  end if
530  end if
531 
532  if (n3.eq.0) then
533  if (n0.gt.0.or.r.le.rmax-1) then
534  s(3) = s(3) + d_i(n0,n1,n2,n4,3)
535  s2(:,3) = 0d0
536  end if
537  else
538  if (r.le.rmax-1) then
539  s2(:,3) = -n3*d_0(n0+1,n1,n2,n3-1,n4)
540  if (n1.eq.0) then
541  s2(1,3) = s2(1,3) + n3*d_i(n0+1,n2,n3-1,n4,1)
542  end if
543  if (n2.eq.0) then
544  s2(2,3) = s2(2,3) + n3*d_i(n0+1,n1,n3-1,n4,2)
545  end if
546  if (n4.eq.0) then
547  s2(4,3) = s2(4,3) + n3*d_i(n0+1,n1,n2,n3-1,4)
548  end if
549  end if
550  end if
551 
552  if (n4.eq.0) then
553  if (n0.gt.0.or.r.le.rmax-1) then
554  s(4) = s(4) + d_i(n0,n1,n2,n3,4)
555  s2(:,4) = 0d0
556  end if
557  else
558  if (r.le.rmax-1) then
559  s2(:,4) = -n4*d_0(n0+1,n1,n2,n3,n4-1)
560  if (n1.eq.0) then
561  s2(1,4) = s2(1,4) + n4*d_i(n0+1,n2,n3,n4-1,1)
562  end if
563  if (n2.eq.0) then
564  s2(2,4) = s2(2,4) + n4*d_i(n0+1,n1,n3,n4-1,2)
565  end if
566  if (n3.eq.0) then
567  s2(3,4) = s2(3,4) + n4*d_i(n0+1,n1,n2,n4-1,3)
568  end if
569  end if
570  end if
571 
572  if (r.le.rmax-1) then
573 
574  ! rational term
575  s(0) = s(0) - 4d0*euv_aux(n0+1,n1,n2,n3,n4)
576 
577  do k=1,4
578  eaux(k) = mxinv(0,k)*s(0)+mxinv(1,k)*s(1)+mxinv(2,k)*s(2) &
579  + mxinv(3,k)*s(3)+mxinv(4,k)*s(4) &
580  - mxinvs(k) * d_0(n0,n1,n2,n3,n4) &
581  + ((mxinv(1,k)*mxinv(2,0)-mxinv(2,k)*mxinv(1,0))*(s2(1,2)-s2(2,1)) &
582  + (mxinv(1,k)*mxinv(3,0)-mxinv(3,k)*mxinv(1,0))*(s2(1,3)-s2(3,1)) &
583  + (mxinv(1,k)*mxinv(4,0)-mxinv(4,k)*mxinv(1,0))*(s2(1,4)-s2(4,1)) &
584  + (mxinv(2,k)*mxinv(3,0)-mxinv(3,k)*mxinv(2,0))*(s2(2,3)-s2(3,2)) &
585  + (mxinv(2,k)*mxinv(4,0)-mxinv(4,k)*mxinv(2,0))*(s2(2,4)-s2(4,2)) &
586  + (mxinv(3,k)*mxinv(4,0)-mxinv(4,k)*mxinv(3,0))*(s2(3,4)-s2(4,3)))*2
587  end do
588 
589  e(n0,n1+1,n2,n3,n4) = e(n0,n1+1,n2,n3,n4) + (n1+1)*eaux(1)/(r+1)
590  e(n0,n1,n2+1,n3,n4) = e(n0,n1,n2+1,n3,n4) + (n2+1)*eaux(2)/(r+1)
591  e(n0,n1,n2,n3+1,n4) = e(n0,n1,n2,n3+1,n4) + (n3+1)*eaux(3)/(r+1)
592  e(n0,n1,n2,n3,n4+1) = e(n0,n1,n2,n3,n4+1) + (n4+1)*eaux(4)/(r+1)
593 
594  end if
595 
596  if (n0.ge.1) then
597  eaux2 = mxinv(1,0)*s(1)+mxinv(2,0)*s(2) &
598  + mxinv(3,0)*s(3)+mxinv(4,0)*s(4) &
599  - (mxinvs(0)-mxinv(0,0)) * d_0(n0,n1,n2,n3,n4)
600 
601  e(n0,n1,n2,n3,n4) = e(n0,n1,n2,n3,n4) + 2*n0*eaux2/r
602  end if
603 
604  end do
605  end do
606  end do
607  end do
608 
609  if (r.le.rmax-1) then
610 ! Eerr(r+1) = max( maxval(abs(mxinv(0:4,1:4)))*Derr(0,r), &
611  eerr(r+1) = max( maxval(abs(mxinvs(1:4)))*derr(0,r), &
612  maxval(abs(mxinv(1,1:4)))*derr(1,r) , &
613  maxval(abs(mxinv(2,1:4)))*derr(2,r) , &
614  maxval(abs(mxinv(3,1:4)))*derr(3,r) , &
615  maxval(abs(mxinv(4,1:4)))*derr(4,r) )
616  eerr2(r+1) = max( abs(maxzmxinvs)*derr2(0,r), &
617  abs(maxzmxinv(1))*derr2(1,r) , &
618  abs(maxzmxinv(2))*derr2(2,r) , &
619  abs(maxzmxinv(3))*derr2(3,r) , &
620  abs(maxzmxinv(4))*derr2(4,r) ) /maxz
621 
622 ! write(*,*) 'CalcEred s ', maxval(abs(mxinvs(1:4)))
623 ! write(*,*) 'CalcEred 1 ', maxval(abs(mxinv(1,1:4)))
624 ! write(*,*) 'CalcEred 2 ', maxval(abs(mxinv(2,1:4)))
625 ! write(*,*) 'CalcEred 3 ', maxval(abs(mxinv(3,1:4)))
626 ! write(*,*) 'CalcEred 4 ', maxval(abs(mxinv(4,1:4)))
627 ! write(*,*) 'CalcEred err ',r,Eerr(r+1)
628 ! if (r.gt.0) write(*,*) 'CalcEred err ',1,Eerr(2)
629 ! write(*,*) 'CalcEred err ',Derr(0:4,r)
630 ! write(*,*) 'CalcEred err2',r,Eerr2(r+1)
631 ! write(*,*) 'CalcEred err2',Derr2(0:4,r)
632 
633  end if
634 
635 
636  end do
637 
638  ! write(*,*) 'E', id, E(0,0,0,0,2)
639 
640 
641  end subroutine calcered
642 
643 
644 
645 
646 
647  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
648  ! subroutine CalcF(F,Fuv,p10,p21,p32,p43,p54,p50,p20,p31,p42,p53,p40, &
649  ! p51,p30,p41,p52,m02,m12,m22,m32,m42,m52, &
650  ! rmax,id,Ferr,Ferr2)
651  !
652  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
653 
654  subroutine calcf(F,Fuv,p10,p21,p32,p43,p54,p50,p20,p31,p42,p53,p40, &
655  p51,p30,p41,p52,m02,m12,m22,m32,m42,m52, &
656  rmax,id,Ferr,Ferr2)
657 
658  integer, intent(in) :: rmax,id
659  double complex, intent(in) :: p10,p21,p32,p43,p54,p50,p20,p31,p42,p53,p40
660  double complex, intent(in) :: p51,p30,p41,p52,m02,m12,m22,m32,m42,m52
661  double complex, intent(out) :: F(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
662  double complex, intent(out) :: Fuv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
663  double precision, intent(out) :: Ferr(0:rmax),Ferr2(0:rmax)
664  double complex, allocatable :: Faux(:,:,:,:,:,:), Fuvaux(:,:,:,:,:,:), fct(:)
665  double precision, allocatable :: Ferraux(:),Ferr2aux(:)
666  double complex :: x(21)
667  integer :: rank,switch,cnt,n0,n1,n2,n3,n4,n5,r
668  logical :: nocalc,wrica
669 
670 ! write(*,*) 'CalcF in'
671 
672  if ((use_cache_system).and.(tencache.gt.6)) then
673  if ((ncache.gt.0).and.(ncache.le.ncache_max)) then
674 ! if (use_cache(ncache).ge.6) then
675  x(1)=p10
676  x(2)=p21
677  x(3)=p32
678  x(4)=p43
679  x(5)=p54
680  x(6)=p50
681  x(7)=p20
682  x(8)=p31
683  x(9)=p42
684  x(10)=p53
685  x(11)=p40
686  x(12)=p51
687  x(13)=p30
688  x(14)=p41
689  x(15)=p52
690  x(16)=m02
691  x(17)=m12
692  x(18)=m22
693  x(19)=m32
694  x(20)=m42
695  x(21)=m52
696 
697  rank = rmax
698 
699  if (rmax.ge.8) then
700  allocate(fct(ncoefs(rmax,6)+ncoefs(rmax-8,6)+2*(rmax+1)))
701  call readcache(fct,ncoefs(rmax,6)+ncoefs(rmax-8,6)+2*(rmax+1),x,21,1,id,6,rank,nocalc,wrica)
702  else
703  allocate(fct(ncoefs(rmax,6)+2*(rmax+1)))
704  call readcache(fct,ncoefs(rmax,6)+2*(rmax+1),x,21,1,id,6,rank,nocalc,wrica)
705  end if
706 
707  if(nocalc) then
708  cnt = 0
709  fuv(0:min(rmax/2,3),:,:,:,:,:) = 0d0
710  do r=0,rmax
711  do n0=0,r/2
712  do n1=0,r-2*n0
713  do n2=0,r-2*n0-n1
714  do n3=0,r-2*n0-n1-n2
715  do n4=0,r-2*n0-n1-n2-n3
716  n5 = r-2*n0-n1-n2-n3-n4
717 
718  cnt = cnt+1
719  f(n0,n1,n2,n3,n4,n5) = fct(cnt)
720 
721  end do
722  end do
723  end do
724  end do
725  end do
726  do n0=4,r/2
727  do n1=0,r-2*n0
728  do n2=0,r-2*n0-n1
729  do n3=0,r-2*n0-n1-n2
730  do n4=0,r-2*n0-n1-n2-n3
731  n5 = r-2*n0-n1-n2-n3-n4
732 
733  cnt = cnt+1
734  fuv(n0,n1,n2,n3,n4,n5) = fct(cnt)
735 
736  end do
737  end do
738  end do
739  end do
740  end do
741  cnt = cnt+1
742  ferr(r) = real(fct(cnt))
743  cnt = cnt+1
744  ferr2(r) = real(fct(cnt))
745  end do
746  return
747  endif
748 
749 
750  if(rank.eq.rmax) then
751 
752  call calcfred(f,fuv,p10,p21,p32,p43,p54,p50,p20,p31,p42,p53,p40, &
753  p51,p30,p41,p52,m02,m12,m22,m32,m42,m52,rank,id,ferr,ferr2)
754 
755  if (wrica) then
756  cnt = 0
757  do r=0,rank
758  do n0=0,r/2
759  do n1=0,r-2*n0
760  do n2=0,r-2*n0-n1
761  do n3=0,r-2*n0-n1-n2
762  do n4=0,r-2*n0-n1-n2-n3
763  n5 = r-2*n0-n1-n2-n3-n4
764 
765  cnt = cnt+1
766  fct(cnt) = f(n0,n1,n2,n3,n4,n5)
767 
768  end do
769  end do
770  end do
771  end do
772  end do
773  do n0=4,r/2
774  do n1=0,r-2*n0
775  do n2=0,r-2*n0-n1
776  do n3=0,r-2*n0-n1-n2
777  do n4=0,r-2*n0-n1-n2-n3
778  n5 = r-2*n0-n1-n2-n3-n4
779 
780  cnt = cnt+1
781  fct(cnt) = fuv(n0,n1,n2,n3,n4,n5)
782 
783  end do
784  end do
785  end do
786  end do
787  end do
788  cnt = cnt+1
789  fct(cnt) = ferr(r)
790  cnt = cnt+1
791  fct(cnt) = ferr2(r)
792  end do
793 
794  if (rank.ge.8) then
795  call writecache(fct,ncoefs(rank,6)+ncoefs(rank-8,6)+2*(rank+1),id,6,rank)
796  else
797  call writecache(fct,ncoefs(rank,6)+2*(rank+1),id,6,rank)
798  end if
799 
800  end if
801 
802  return
803 
804 
805  else
806  allocate(faux(0:rank/2,0:rank,0:rank,0:rank,0:rank,0:rank))
807  allocate(fuvaux(0:rank/2,0:rank,0:rank,0:rank,0:rank,0:rank))
808  allocate(ferraux(0:rank))
809  allocate(ferr2aux(0:rank))
810 
811  call calcfred(faux,fuvaux,p10,p21,p32,p43,p54,p50,p20,p31,p42,p53,p40, &
812  p51,p30,p41,p52,m02,m12,m22,m32,m42,m52,rank,id,ferraux,ferr2aux)
813 
814 
815  if (wrica) then
816  cnt = 0
817  deallocate(fct)
818  if (rank.ge.8) then
819  allocate(fct(ncoefs(rank,6)+ncoefs(rank-8,6)+2*(rank+1)))
820  else
821  allocate(fct(ncoefs(rank,6)+2*(rank+1)))
822  end if
823  do r=0,rank
824  do n0=0,r/2
825  do n1=0,r-2*n0
826  do n2=0,r-2*n0-n1
827  do n3=0,r-2*n0-n1-n2
828  do n4=0,r-2*n0-n1-n2-n3
829  n5 = r-2*n0-n1-n2-n3-n4
830 
831  cnt = cnt+1
832  fct(cnt) = faux(n0,n1,n2,n3,n4,n5)
833 
834  end do
835  end do
836  end do
837  end do
838  end do
839  do n0=4,r/2
840  do n1=0,r-2*n0
841  do n2=0,r-2*n0-n1
842  do n3=0,r-2*n0-n1-n2
843  do n4=0,r-2*n0-n1-n2-n3
844  n5 = r-2*n0-n1-n2-n3-n4
845 
846  cnt = cnt+1
847  fct(cnt) = fuvaux(n0,n1,n2,n3,n4,n5)
848 
849  end do
850  end do
851  end do
852  end do
853  end do
854  cnt = cnt+1
855  fct(cnt) = ferraux(r)
856  cnt = cnt+1
857  fct(cnt) = ferr2aux(r)
858  end do
859 
860  if (rank.ge.8) then
861  call writecache(fct,ncoefs(rank,6)+ncoefs(rank-8,6)+2*(rank+1),id,6,rank)
862  else
863  call writecache(fct,ncoefs(rank,6)+2*(rank+1),id,6,rank)
864  end if
865 
866  end if
867 
868  f = faux(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
869  fuv = fuvaux(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
870  ferr = ferraux(0:rmax)
871  ferr2 = ferr2aux(0:rmax)
872  return
873 
874  end if
875 
876 ! end if
877  end if
878  end if
879 
880 
881  call calcfred(f,fuv,p10,p21,p32,p43,p54,p50,p20,p31,p42,p53,p40, &
882  p51,p30,p41,p52,m02,m12,m22,m32,m42,m52,rmax,id,ferr,ferr2)
883 
884 
885  end subroutine calcf
886 
887 
888 
889 
890 
891  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
892  ! subroutine CalcFred(F,Fuv,p10,p21,p32,p43,p54,p50,p20,p31,p42,p53,p40, &
893  ! p51,p30,p41,p52,m02,m12,m22,m32,m42,m52, &
894  ! rmax,id,Ferr,Ferr2)
895  !
896  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
897 
898  subroutine calcfred(F,Fuv,p10,p21,p32,p43,p54,p50,p20,p31,p42,p53,p40, &
899  p51,p30,p41,p52,m02,m12,m22,m32,m42,m52, &
900  rmax,id,Ferr,Ferr2)
901 
902  integer, intent(in) :: rmax,id
903  double complex, intent(in) :: p10,p21,p32,p43,p54,p50,p20,p31,p42,p53,p40
904  double complex, intent(in) :: p51,p30,p41,p52,m02,m12,m22,m32,m42,m52
905  double precision, intent(out) :: Ferr(0:rmax),Ferr2(0:rmax)
906  double complex, intent(out) :: F(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
907  double complex, intent(out) :: Fuv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
908  double complex :: q10,q21,q32,q43,q54,q50,q20,q31,q42,q53,q40
909  double complex :: q51,q30,q41,q52,mm02,mm12,mm22,mm32,mm42,mm52
910  double complex :: mx(0:5,0:5), mxinv(0:5,0:5),mx0k(5,5),mx0kinv(5,5),ff(5)
911  double complex :: det,newdet,mxinvs,mx0kinvs(5)
912  double complex :: zmx0kinv(5,5),Z(5,5),zmx0kinvs(5)
913  double precision :: maxZ,maxzmx0kinv(5),maxzmx0kinvs
914  double complex, allocatable :: E_0(:,:,:,:,:,:), E_i(:,:,:,:,:,:)
915  double complex, allocatable :: Euv_0(:,:,:,:,:,:), Euv_i(:,:,:,:,:,:)
916  double complex :: S(5), Faux(5), elimminf2_coli,chdet,gram(1:5,1:5),gramdet
917  double precision :: Eerr(0:5,0:rmax),Eerr2(0:5,0:rmax)
918  integer :: r,n0,n1,n2,n3,n4,n5,n,k,i,j,nid(0:5),r0,bin,kbest,rmaxE,rBCD
919  logical :: errorwriteflag
920  character(len=*),parameter :: fmt10 = "(A17,'(',d25.18,' , ',d25.18,' )')"
921 
922 ! write(*,*) 'CalcFred in'
923 
924  r0 = 0
925 
926  ! allocation of C functions
927  rmaxe = max(rmax-1,0)
928  allocate(e_0(0:rmaxe/2,0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe))
929  allocate(e_i(0:rmaxe/2,0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe,5))
930  allocate(euv_0(0:rmaxe/2,0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe))
931  allocate(euv_i(0:rmaxe/2,0:rmaxe,0:rmaxe,0:rmaxe,0:rmaxe,5))
932 
933 
934  ! determine binaries for B-coefficients
935  k=0
936  bin = 1
937  do while (k.le.5)
938  if (mod(id/bin,2).eq.0) then
939  nid(k) = id+bin
940  k = k+1
941  end if
942  bin = 2*bin
943  end do
944 
945 
946  call calce(e_0(:,0,:,:,:,:),euv_0(:,0,:,:,:,:),p21,p32,p43,p54,p51,p31,p42,p53,p41,p52, &
947  m12,m22,m32,m42,m52,rmaxe,nid(0),eerr=eerr(0,0:rmaxe),eerr2=eerr2(0,0:rmaxe))
948  call calce(e_i(:,:,:,:,:,1),euv_i(:,:,:,:,:,1),p20,p32,p43,p54,p50,p30, &
949  p42,p53,p40,p52,m02,m22,m32,m42,m52,rmaxe,nid(1),eerr=eerr(1,0:rmaxe),eerr2=eerr2(1,0:rmaxe))
950  call calce(e_i(:,:,:,:,:,2),euv_i(:,:,:,:,:,2),p10,p31,p43,p54,p50,p30, &
951  p41,p53,p40,p51,m02,m12,m32,m42,m52,rmaxe,nid(2),eerr=eerr(2,0:rmaxe),eerr2=eerr2(2,0:rmaxe))
952  call calce(e_i(:,:,:,:,:,3),euv_i(:,:,:,:,:,3),p10,p21,p42,p54,p50,p20, &
953  p41,p52,p40,p51,m02,m12,m22,m42,m52,rmaxe,nid(3),eerr=eerr(3,0:rmaxe),eerr2=eerr2(3,0:rmaxe))
954  call calce(e_i(:,:,:,:,:,4),euv_i(:,:,:,:,:,4),p10,p21,p32,p53,p50,p20, &
955  p31,p52,p30,p51,m02,m12,m22,m32,m52,rmaxe,nid(4),eerr=eerr(4,0:rmaxe),eerr2=eerr2(4,0:rmaxe))
956  call calce(e_i(:,:,:,:,:,5),euv_i(:,:,:,:,:,5),p10,p21,p32,p43,p40,p20, &
957  p31,p42,p30,p41,m02,m12,m22,m32,m42,rmaxe,nid(5),eerr=eerr(5,0:rmaxe),eerr2=eerr2(5,0:rmaxe))
958 
959  ! shift of integration momentum in E\{0}
960  do n1=1,rmaxe
961  do n2=0,rmaxe-n1
962  do n3=0,rmaxe-n1-n2
963  do n4=0,rmaxe-n1-n2-n3
964  do n5=0,rmaxe-n1-n2-n3-n4
965  n0 = (rmaxe-n1-n2-n3-n4-n5)/2
966  e_0(0:n0,n1,n2,n3,n4,n5) = -e_0(0:n0,n1-1,n2,n3,n4,n5)-e_0(0:n0,n1-1,n2+1,n3,n4,n5) &
967  -e_0(0:n0,n1-1,n2,n3+1,n4,n5)-e_0(0:n0,n1-1,n2,n3,n4+1,n5) &
968  -e_0(0:n0,n1-1,n2,n3,n4,n5+1)
969  end do
970  end do
971  end do
972  end do
973  end do
974  do n1=1,rmaxe-6
975  do n2=0,rmaxe-n1-6
976  do n3=0,rmaxe-n1-n2-6
977  do n4=0,rmaxe-n1-n2-n3-6
978  do n5=0,rmaxe-n1-n2-n3-n4-6
979  n0 = (rmaxe-n1-n2-n3-n4-n5)/2
980  euv_0(0:n0,n1,n2,n3,n4,n5) = -euv_0(0:n0,n1-1,n2,n3,n4,n5)-euv_0(0:n0,n1-1,n2+1,n3,n4,n5) &
981  -euv_0(0:n0,n1-1,n2,n3+1,n4,n5)-euv_0(0:n0,n1-1,n2,n3,n4+1,n5) &
982  -euv_0(0:n0,n1-1,n2,n3,n4,n5+1)
983  end do
984  end do
985  end do
986  end do
987  end do
988 
989 
990 
991  ! determine inverse modified Caley matrix
992  mm02 = elimminf2_coli(m02)
993  mm12 = elimminf2_coli(m12)
994  mm22 = elimminf2_coli(m22)
995  mm32 = elimminf2_coli(m32)
996  mm42 = elimminf2_coli(m42)
997  mm52 = elimminf2_coli(m52)
998  q10 = elimminf2_coli(p10)
999  q21 = elimminf2_coli(p21)
1000  q32 = elimminf2_coli(p32)
1001  q43 = elimminf2_coli(p43)
1002  q54 = elimminf2_coli(p54)
1003  q50 = elimminf2_coli(p50)
1004  q20 = elimminf2_coli(p20)
1005  q31 = elimminf2_coli(p31)
1006  q42 = elimminf2_coli(p42)
1007  q53 = elimminf2_coli(p53)
1008  q40 = elimminf2_coli(p40)
1009  q51 = elimminf2_coli(p51)
1010  q30 = elimminf2_coli(p30)
1011  q41 = elimminf2_coli(p41)
1012  q52 = elimminf2_coli(p52)
1013 
1014  ff(1) = q10+mm02-mm12
1015  ff(2) = q20+mm02-mm22
1016  ff(3) = q30+mm02-mm32
1017  ff(4) = q40+mm02-mm42
1018  ff(5) = q50+mm02-mm52
1019 
1020  mx(0,0) = 2d0*mm02
1021  mx(1,0) = q10 - mm12 + mm02
1022  mx(2,0) = q20 - mm22 + mm02
1023  mx(3,0) = q30 - mm32 + mm02
1024  mx(4,0) = q40 - mm42 + mm02
1025  mx(5,0) = q50 - mm52 + mm02
1026  mx(0,1) = mx(1,0)
1027  mx(1,1) = 2d0*q10
1028  mx(2,1) = q10+q20-q21
1029  mx(3,1) = q10+q30-q31
1030  mx(4,1) = q10+q40-q41
1031  mx(5,1) = q10+q50-q51
1032  mx(0,2) = mx(2,0)
1033  mx(1,2) = mx(2,1)
1034  mx(2,2) = 2d0*q20
1035  mx(3,2) = q20+q30-q32
1036  mx(4,2) = q20+q40-q42
1037  mx(5,2) = q20+q50-q52
1038  mx(0,3) = mx(3,0)
1039  mx(1,3) = mx(3,1)
1040  mx(2,3) = mx(3,2)
1041  mx(3,3) = 2d0*q30
1042  mx(4,3) = q30+q40-q43
1043  mx(5,3) = q30+q50-q53
1044  mx(0,4) = mx(4,0)
1045  mx(1,4) = mx(4,1)
1046  mx(2,4) = mx(4,2)
1047  mx(3,4) = mx(4,3)
1048  mx(4,4) = 2d0*q40
1049  mx(5,4) = q40+q50-q54
1050  mx(0,5) = mx(5,0)
1051  mx(1,5) = mx(5,1)
1052  mx(2,5) = mx(5,2)
1053  mx(3,5) = mx(5,3)
1054  mx(4,5) = mx(5,4)
1055  mx(5,5) = 2d0*q50
1056 
1057  call chinv(6,mx,mxinv)
1058 
1059  ! determine X_(0,5)
1060  do j=1,5
1061  do i=1,5
1062  mx0k(i,j) = mx(i,j-1)
1063  end do
1064  end do
1065 
1066  ! determine best reduction matrix
1067  det = chdet(5,mx0k)
1068  kbest = 5
1069 
1070 ! write(*,*) 'Fred det',5,det
1071 
1072  do j=5,2,-1
1073  do i=1,5
1074  mx0k(i,j) = mx(i,j)
1075  end do
1076 
1077  newdet = chdet(5,mx0k)
1078 
1079 ! write(*,*) 'Fred det',j-1,newdet
1080 
1081  if (abs(newdet).gt.abs(det)) then
1082  kbest = j-1
1083  det = newdet
1084  end if
1085 
1086  end do
1087 
1088 ! write(*,*) 'Fred kbest',kbest
1089 
1090  do i=1,5
1091  mx0k(i,1) = mx(i,1)
1092  mx0k(i,kbest) = mx(i,0)
1093  end do
1094 
1095  call chinv(5,mx0k,mx0kinv)
1096  do i=1,5
1097  mx0kinv(kbest,i) = 0d0
1098  end do
1099 
1100  mxinvs = sum(mxinv(0:5,0))
1101  do i=1,5
1102  mx0kinvs(i) = sum(mx0kinv(i,1:5))
1103  end do
1104 
1105  ! for alternative error estimate
1106  z(1:5,1:5) = mx(1:5,1:5)
1107 
1108  zmx0kinv = matmul(z,mx0kinv)
1109 
1110  do i=1,5
1111  maxzmx0kinv(i) = maxval(abs(zmx0kinv(1:5,i)))
1112  zmx0kinvs(i) = sum(zmx0kinv(i,1:5))
1113  end do
1114 
1115  maxzmx0kinvs = maxval(abs(zmx0kinvs(1:5)))
1116 
1117  maxz = maxval(abs(z))
1118 
1119 
1120 
1121  ! calculation of UV-divergent parts
1122  ! F_(n0,n1,n2,n3) UV-finite for n0<4
1123  fuv(0:min(rmax/2,3),:,:,:,:,:) = 0d0
1124 
1125  ! PV reduction (5.10)
1126  do r=max(r0,8),rmax
1127  do n0=3,r/2
1128  do n1=0,r-2*n0
1129  do n2=0,r-2*n0-n1
1130  do n3=0,r-2*n0-n1-n2
1131  do n4=0,r-2*n0-n1-n2-n3
1132  n5 = r-2*n0-n1-n2-n3-n4
1133 
1134  fuv(n0,n1,n2,n3,n4,n5) = (euv_0(n0-1,n1,n2,n3,n4,n5) &
1135  + 2*m02*fuv(n0-1,n1,n2,n3,n4,n5) &
1136  + ff(1)*fuv(n0-1,n1+1,n2,n3,n4,n5) &
1137  + ff(2)*fuv(n0-1,n1,n2+1,n3,n4,n5) &
1138  + ff(3)*fuv(n0-1,n1,n2,n3+1,n4,n5) &
1139  + ff(4)*fuv(n0-1,n1,n2,n3,n4+1,n5) &
1140  + ff(5)*fuv(n0-1,n1,n2,n3,n4,n5+1))/ (2*(r-3))
1141 
1142  end do
1143  end do
1144  end do
1145  end do
1146  end do
1147  end do
1148 
1149  ! scalar coefficient
1150  if (r0.eq.0) then
1151  f = 0d0
1152  f(0,0,0,0,0,0) = -mxinv(0,0)*e_0(0,0,0,0,0,0)
1153  do k=1,5
1154  f(0,0,0,0,0,0) = f(0,0,0,0,0,0) &
1155  + mxinv(k,0)*(e_i(0,0,0,0,0,k)-e_0(0,0,0,0,0,0))
1156  end do
1157  end if
1158 
1159 ! write(*,*) 'CalcFred: Eerr(0)',Eerr(0,0:rmaxE)
1160 ! write(*,*) 'CalcFred: Eerr(1)',Eerr(1,0:rmaxE)
1161 ! write(*,*) 'CalcFred: Eerr(2)',Eerr(2,0:rmaxE)
1162 ! write(*,*) 'CalcFred: Eerr(3)',Eerr(3,0:rmaxE)
1163 ! write(*,*) 'CalcFred: Eerr(4)',Eerr(4,0:rmaxE)
1164 ! write(*,*) 'CalcFred: Eerr(5)',Eerr(5,0:rmaxE)
1165 ! write(*,*) 'CalcFred: mxinv(0)',mxinv(0,0)
1166 ! write(*,*) 'CalcFred: mxinv(1)',mxinv(1,0)
1167 ! write(*,*) 'CalcFred: mxinv(2)',mxinv(2,0)
1168 ! write(*,*) 'CalcFred: mxinv(3)',mxinv(3,0)
1169 ! write(*,*) 'CalcFred: mxinv(4)',mxinv(4,0)
1170 ! write(*,*) 'CalcFred: mxinv(5)',mxinv(5,0)
1171 
1172  ferr(0) = max( abs(mxinvs)*eerr(0,0), &
1173  abs(mxinv(1,0))*eerr(1,0) , &
1174  abs(mxinv(2,0))*eerr(2,0) , &
1175  abs(mxinv(3,0))*eerr(3,0) , &
1176  abs(mxinv(4,0))*eerr(4,0) , &
1177  abs(mxinv(5,0))*eerr(5,0) )
1178  ferr2(0) = max( abs(mxinvs)*eerr2(0,0), &
1179  abs(mxinv(1,0))*eerr2(1,0) , &
1180  abs(mxinv(2,0))*eerr2(2,0) , &
1181  abs(mxinv(3,0))*eerr2(3,0) , &
1182  abs(mxinv(4,0))*eerr2(4,0) , &
1183  abs(mxinv(5,0))*eerr2(5,0) )
1184 
1185 ! write(*,*) 'CalcFred: Ferr(0)',Ferr(0)
1186 
1187  ! formula (7.13)
1188  do r=r0,rmax-1
1189  do n0=0,r/2
1190  do n1=0,r-2*n0
1191  do n2=0,r-2*n0-n1
1192  do n3=0,r-2*n0-n1-n2
1193  do n4=0,r-2*n0-n1-n2-n3
1194  n5 = r-2*n0-n1-n2-n3-n4
1195 
1196  do n=1,5
1197  s(n) = -e_0(n0,n1,n2,n3,n4,n5)
1198  end do
1199 
1200  if (n1.eq.0) then
1201  s(1) = s(1) + e_i(n0,n2,n3,n4,n5,1)
1202  end if
1203  if (n2.eq.0) then
1204  s(2) = s(2) + e_i(n0,n1,n3,n4,n5,2)
1205  end if
1206  if (n3.eq.0) then
1207  s(3) = s(3) + e_i(n0,n1,n2,n4,n5,3)
1208  end if
1209  if (n4.eq.0) then
1210  s(4) = s(4) + e_i(n0,n1,n2,n3,n5,4)
1211  end if
1212  if (n5.eq.0) then
1213  s(5) = s(5) + e_i(n0,n1,n2,n3,n4,5)
1214  end if
1215 
1216  do k=1,5
1217  faux(k) = mx0kinv(k,1)*s(1)+mx0kinv(k,2)*s(2) &
1218  + mx0kinv(k,3)*s(3)+mx0kinv(k,4)*s(4)+mx0kinv(k,5)*s(5)
1219  end do
1220 
1221  f(n0,n1+1,n2,n3,n4,n5) = f(n0,n1+1,n2,n3,n4,n5) + (n1+1)*faux(1)/(r+1)
1222  f(n0,n1,n2+1,n3,n4,n5) = f(n0,n1,n2+1,n3,n4,n5) + (n2+1)*faux(2)/(r+1)
1223  f(n0,n1,n2,n3+1,n4,n5) = f(n0,n1,n2,n3+1,n4,n5) + (n3+1)*faux(3)/(r+1)
1224  f(n0,n1,n2,n3,n4+1,n5) = f(n0,n1,n2,n3,n4+1,n5) + (n4+1)*faux(4)/(r+1)
1225  f(n0,n1,n2,n3,n4,n5+1) = f(n0,n1,n2,n3,n4,n5+1) + (n5+1)*faux(5)/(r+1)
1226 
1227  end do
1228  end do
1229  end do
1230  end do
1231  end do
1232 
1233  if (r.le.rmax-1) then
1234  ferr(r+1) = max( maxval(abs(mx0kinvs(1:5)))*eerr(0,r), &
1235  maxval(abs(mx0kinv(1:5,1)))*eerr(1,r) , &
1236  maxval(abs(mx0kinv(1:5,2)))*eerr(2,r) , &
1237  maxval(abs(mx0kinv(1:5,3)))*eerr(3,r) , &
1238  maxval(abs(mx0kinv(1:5,4)))*eerr(4,r) , &
1239  maxval(abs(mx0kinv(1:5,5)))*eerr(5,r) )
1240  ferr2(r+1) = max( abs(maxzmx0kinvs)*eerr2(0,r), &
1241  abs(maxzmx0kinv(1))*eerr2(1,r) , &
1242  abs(maxzmx0kinv(2))*eerr2(2,r) , &
1243  abs(maxzmx0kinv(3))*eerr2(3,r) , &
1244  abs(maxzmx0kinv(4))*eerr2(4,r) , &
1245  abs(maxzmx0kinv(5))*eerr2(5,r) )/maxz
1246 
1247 ! write(*,*) 'CalcFred Ferr',r, Ferr(r+1)
1248 ! write(*,*) 'CalcFred Ferr',Eerr(0:5,r)
1249 ! write(*,*) 'CalcFred Ferr2',r, Ferr2(r+1)
1250 ! write(*,*) 'CalcFred Ferr2',Eerr2(0:5,r)
1251 ! write(*,*) 'CalcFred Ferr',maxval(abs(mx0kinv(1:5,1:5))),maxval(abs(mx0kinvs(1:5)))
1252 ! write(*,*) 'CalcFred Ferr',maxval(abs(mx0kinvs(1:5)))*Eerr(0,r)
1253 ! write(*,*) 'CalcFred Ferr',maxval(abs(mx0kinv(1:5,1:5)))*Eerr(0,r), &
1254 ! maxval(abs(mx0kinv(1:5,1)))*Eerr(1,r) , &
1255 ! maxval(abs(mx0kinv(1:5,2)))*Eerr(2,r) , &
1256 ! maxval(abs(mx0kinv(1:5,3)))*Eerr(3,r) , &
1257 ! maxval(abs(mx0kinv(1:5,4)))*Eerr(4,r) , &
1258 ! maxval(abs(mx0kinv(1:5,5)))*Eerr(5,r)
1259 ! write(*,*) 'CalcFred Ferr',Eerr(0,r), &
1260 ! Eerr(1,r) , &
1261 ! Eerr(2,r) , &
1262 ! Eerr(3,r) , &
1263 ! Eerr(4,r) , &
1264 ! Eerr(5,r)
1265 ! write(*,*) 'CalcFred Ferr',maxval(abs(mxinv(1:5,1:5))), &
1266 ! maxval(abs(mx0kinv(1:5,1))), &
1267 ! maxval(abs(mx0kinv(1:5,2))), &
1268 ! maxval(abs(mx0kinv(1:5,3))), &
1269 ! maxval(abs(mx0kinv(1:5,4))), &
1270 ! maxval(abs(mx0kinv(1:5,5)))
1271 
1272  if (mode_coli.lt.1) then
1273  gram= mx(1:5,1:5)
1274 
1275  gramdet= chdet(5,gram)
1276 
1277 ! write(*,*) 'CalcFred gram=',det
1278 
1279  if (max(abs(f(0,0,0,0,0,0)),abs(f(0,1,0,0,0,0)),abs(f(0,0,1,0,0,0)), &
1280  abs(f(0,0,0,1,0,0)),abs(f(0,0,0,0,1,0)),abs(f(0,0,0,0,0,1)) &
1281  )*abs(gramdet/det).gt. ferr(r+1)) then
1282 
1283 
1284 ! write(*,*) 'CalcFred Ferr=',r+1,Ferr(r+1)
1285 
1286  ferr(r+1)=max(ferr(r+1),max( &
1287  abs(f(0,0,0,0,0,0)),abs(f(0,1,0,0,0,0)), &
1288  abs(f(0,0,1,0,0,0)), abs(f(0,0,0,1,0,0)), &
1289  abs(f(0,0,0,0,1,0)),abs(f(0,0,0,0,0,1)) &
1290  )*abs(gramdet/det) )
1291 
1292  ferr2(r+1)=max(ferr2(r+1),max( &
1293  abs(f(0,0,0,0,0,0)),abs(f(0,1,0,0,0,0)), &
1294  abs(f(0,0,1,0,0,0)), abs(f(0,0,0,1,0,0)), &
1295  abs(f(0,0,0,0,1,0)),abs(f(0,0,0,0,0,1)) &
1296  )*abs(gramdet/det) )
1297 
1298 ! write(*,*) 'CalcFred ',abs(F(0,0,0,0,0,0)),abs(F(0,1,0,0,0,0)), &
1299 ! abs(F(0,0,1,0,0,0)), abs(F(0,0,0,1,0,0)), &
1300 ! abs(F(0,0,0,0,1,0)),abs(F(0,0,0,0,0,1)),abs(gramdet/det)
1301 ! write(*,*) 'CalcFred Ferr=',r+1,Ferr(r+1)
1302 
1303  if (abs(gramdet/det).gt.reqacc_coli) then
1304  call seterrflag_coli(-6)
1305  call errout_coli('CalcFred', &
1306  'input momenta inconsistent! (not 4-dimensional)', &
1307  errorwriteflag)
1308  if (errorwriteflag) then
1309  write(nerrout_coli,fmt10) ' CalcFred: q10 = ',q10
1310  write(nerrout_coli,fmt10) ' CalcFred: q21 = ',q21
1311  write(nerrout_coli,fmt10) ' CalcFred: q32 = ',q32
1312  write(nerrout_coli,fmt10) ' CalcFred: q43 = ',q43
1313  write(nerrout_coli,fmt10) ' CalcFred: q54 = ',q54
1314  write(nerrout_coli,fmt10) ' CalcFred: q50 = ',q50
1315  write(nerrout_coli,fmt10) ' CalcFred: q20 = ',q10
1316  write(nerrout_coli,fmt10) ' CalcFred: q31 = ',q31
1317  write(nerrout_coli,fmt10) ' CalcFred: q42 = ',q42
1318  write(nerrout_coli,fmt10) ' CalcFred: q53 = ',q53
1319  write(nerrout_coli,fmt10) ' CalcFred: q40 = ',q40
1320  write(nerrout_coli,fmt10) ' CalcFred: q51 = ',q51
1321  write(nerrout_coli,fmt10) ' CalcFred: q30 = ',q30
1322  write(nerrout_coli,fmt10) ' CalcFred: q41 = ',q41
1323  write(nerrout_coli,fmt10) ' CalcFred: q52 = ',q52
1324  write(nerrout_coli,fmt10) ' CalcFred: mm02 = ',mm02
1325  write(nerrout_coli,fmt10) ' CalcFred: mm12 = ',mm12
1326  write(nerrout_coli,fmt10) ' CalcFred: mm22 = ',mm22
1327  write(nerrout_coli,fmt10) ' CalcFred: mm32 = ',mm32
1328  write(nerrout_coli,fmt10) ' CalcFred: mm42 = ',mm42
1329  write(nerrout_coli,fmt10) ' CalcFred: mm52 = ',mm52
1330  write(nerrout_coli,fmt10) ' CalcFred: gram = ',gramdet/det
1331  end if
1332  end if
1333 
1334  end if
1335  end if
1336 
1337 ! write(*,*) 'CalcFred Ferr',r, Ferr(r+1)
1338 ! write(*,*) 'CalcFred Ferr',Eerr(0:5,r)
1339 ! write(*,*) 'CalcFred Ferr2',r, Ferr2(r+1)
1340 ! write(*,*) 'CalcFred Ferr2',Eerr2(0:5,r)
1341 
1342  end if
1343 
1344  end do
1345 
1346  end subroutine calcfred
1347 
1348 
1349 
1350 
1351  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1352  ! subroutine CalcG(G,Guv,p10,p21,p32,p43,p54,p65,p60,p20,p31,p42,p53, &
1353  ! p64,p50,p61,p30,p41,p52,p63,p40,p51,p62, &
1354  ! m02,m12,m22,m32,m42,m52,m62,rmax,id,Gerr,Gerr2)
1355  !
1356  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1357 
1358  subroutine calcg(G,Guv,p10,p21,p32,p43,p54,p65,p60,p20,p31,p42,p53, &
1359  p64,p50,p61,p30,p41,p52,p63,p40,p51,p62, &
1360  m02,m12,m22,m32,m42,m52,m62,rmax,id,Gerr,Gerr2)
1361 
1362  integer, intent(in) :: rmax,id
1363  double complex, intent(in) :: p10,p21,p32,p43,p54,p65,p60,p20,p31,p42,p53
1364  double complex, intent(in) :: p64,p50,p61,p30,p41,p52,p63,p40,p51,p62
1365  double complex, intent(in) :: m02,m12,m22,m32,m42,m52,m62
1366  double complex, intent(out) :: G(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
1367  double complex, intent(out) :: Guv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
1368  double precision, intent(out) :: Gerr(0:rmax),Gerr2(0:rmax)
1369  double complex, allocatable :: Gaux(:,:,:,:,:,:,:), Guvaux(:,:,:,:,:,:,:), fct(:)
1370  double precision, allocatable :: Gerraux(:),Gerr2aux(:)
1371  double complex :: x(28)
1372  integer :: rank,switch,cnt,n0,n1,n2,n3,n4,n5,n6,r,r0
1373  logical :: nocalc,wrica,noten
1374 
1375  r0 = 0
1376 
1377  if ((use_cache_system).and.(tencache.gt.7)) then
1378  if ((ncache.gt.0).and.(ncache.le.ncache_max)) then
1379 ! if (use_cache(ncache).ge.7) then
1380  x(1) = p10
1381  x(2) = p21
1382  x(3) = p32
1383  x(4) = p43
1384  x(5) = p54
1385  x(6) = p65
1386  x(7) = p60
1387  x(8) = p20
1388  x(9) = p31
1389  x(10) = p42
1390  x(11) = p53
1391  x(12) = p64
1392  x(13) = p50
1393  x(14) = p61
1394  x(15) = p30
1395  x(16) = p41
1396  x(17) = p52
1397  x(18) = p63
1398  x(19) = p40
1399  x(20) = p51
1400  x(21) = p62
1401  x(22) = m02
1402  x(23) = m12
1403  x(24) = m22
1404  x(25) = m32
1405  x(26) = m42
1406  x(27) = m52
1407  x(28) = m62
1408 
1409  rank = rmax
1410 
1411  if (rmax.ge.10) then
1412  allocate(fct(ncoefs(rmax,7)+ncoefs(rmax-10,7)+2*(rmax+1)))
1413  call readcache(fct,ncoefs(rmax,7)+ncoefs(rmax-10,7)+2*(rmax+1),x,28,1,id,7,rank,nocalc,wrica)
1414  else
1415  allocate(fct(ncoefs(rmax,7)+2*(rmax+1)))
1416  call readcache(fct,ncoefs(rmax,7)+2*(rmax+1),x,28,1,id,7,rank,nocalc,wrica)
1417  end if
1418 
1419  if(nocalc) then
1420  cnt = 0
1421  guv(0:min(rmax/2,4),:,:,:,:,:,:) = 0d0
1422  do r=0,rmax
1423  do n0=0,r/2
1424  do n1=0,r-2*n0
1425  do n2=0,r-2*n0-n1
1426  do n3=0,r-2*n0-n1-n2
1427  do n4=0,r-2*n0-n1-n2-n3
1428  do n5=0,r-2*n0-n1-n2-n3-n4
1429  n6 = r-2*n0-n1-n2-n3-n4-n5
1430 
1431  cnt = cnt+1
1432  g(n0,n1,n2,n3,n4,n5,n6) = fct(cnt)
1433 
1434  end do
1435  end do
1436  end do
1437  end do
1438  end do
1439  end do
1440  do n0=5,r/2
1441  do n1=0,r-2*n0
1442  do n2=0,r-2*n0-n1
1443  do n3=0,r-2*n0-n1-n2
1444  do n4=0,r-2*n0-n1-n2-n3
1445  do n5=0,r-2*n0-n1-n2-n3-n4
1446  n6 = r-2*n0-n1-n2-n3-n4-n5
1447 
1448  cnt = cnt+1
1449  guv(n0,n1,n2,n3,n4,n5,n6) = fct(cnt)
1450 
1451  end do
1452  end do
1453  end do
1454  end do
1455  end do
1456  end do
1457  cnt = cnt+1
1458  gerr(r) = real(fct(cnt))
1459  cnt = cnt+1
1460  gerr2(r) = real(fct(cnt))
1461  end do
1462  return
1463  endif
1464 
1465 
1466  if(rank.eq.rmax) then
1467 
1468  call calcgred(g,guv,p10,p21,p32,p43,p54,p65,p60,p20,p31,p42,p53, &
1469  p64,p50,p61,p30,p41,p52,p63,p40,p51,p62, &
1470  m02,m12,m22,m32,m42,m52,m62,rank,id,gerr,gerr2)
1471 
1472  if (wrica) then
1473  cnt = 0
1474  do r=0,rank
1475  do n0=0,r/2
1476  do n1=0,r-2*n0
1477  do n2=0,r-2*n0-n1
1478  do n3=0,r-2*n0-n1-n2
1479  do n4=0,r-2*n0-n1-n2-n3
1480  do n5=0,r-2*n0-n1-n2-n3-n4
1481  n6 = r-2*n0-n1-n2-n3-n4-n5
1482 
1483  cnt = cnt+1
1484  fct(cnt) = g(n0,n1,n2,n3,n4,n5,n6)
1485 
1486  end do
1487  end do
1488  end do
1489  end do
1490  end do
1491  end do
1492  do n0=5,r/2
1493  do n1=0,r-2*n0
1494  do n2=0,r-2*n0-n1
1495  do n3=0,r-2*n0-n1-n2
1496  do n4=0,r-2*n0-n1-n2-n3
1497  do n5=0,r-2*n0-n1-n2-n3-n4
1498  n6 = r-2*n0-n1-n2-n3-n4-n5
1499 
1500  cnt = cnt+1
1501  fct(cnt) = guv(n0,n1,n2,n3,n4,n5,n6)
1502 
1503  end do
1504  end do
1505  end do
1506  end do
1507  end do
1508  end do
1509  cnt = cnt+1
1510  fct(cnt) = gerr(r)
1511  cnt = cnt+1
1512  fct(cnt) = gerr2(r)
1513  end do
1514 
1515  if (rmax.ge.10) then
1516  call writecache(fct,ncoefs(rank,7)+ncoefs(rank-10,7)+2*(rank+1),id,7,rank)
1517  else
1518  call writecache(fct,ncoefs(rank,7)+2*(rank+1),id,7,rank)
1519  end if
1520 
1521  end if
1522 
1523  return
1524 
1525 
1526  else
1527  allocate(gaux(0:rank/2,0:rank,0:rank,0:rank,0:rank,0:rank,0:rank))
1528  allocate(guvaux(0:rank/2,0:rank,0:rank,0:rank,0:rank,0:rank,0:rank))
1529  allocate(gerraux(0:rank))
1530  allocate(gerr2aux(0:rank))
1531 
1532  call calcgred(gaux,guvaux,p10,p21,p32,p43,p54,p65,p60,p20,p31,p42,p53, &
1533  p64,p50,p61,p30,p41,p52,p63,p40,p51,p62, &
1534  m02,m12,m22,m32,m42,m52,m62,rank,id,gerraux,gerr2aux)
1535 
1536  if (wrica) then
1537  cnt = 0
1538  if (rmax.ge.10) then
1539  deallocate(fct)
1540  allocate(fct(ncoefs(rank,7)+ncoefs(rank-10,7)+2*(rank+1)))
1541  else
1542  deallocate(fct)
1543  allocate(fct(ncoefs(rank,7)+2*(rank+1)))
1544  end if
1545  do r=0,rank
1546  do n0=0,r/2
1547  do n1=0,r-2*n0
1548  do n2=0,r-2*n0-n1
1549  do n3=0,r-2*n0-n1-n2
1550  do n4=0,r-2*n0-n1-n2-n3
1551  do n5=0,r-2*n0-n1-n2-n3-n4
1552  n6 = r-2*n0-n1-n2-n3-n4-n5
1553 
1554  cnt = cnt+1
1555  fct(cnt) = gaux(n0,n1,n2,n3,n4,n5,n6)
1556 
1557  end do
1558  end do
1559  end do
1560  end do
1561  end do
1562  end do
1563  do n0=5,r/2
1564  do n1=0,r-2*n0
1565  do n2=0,r-2*n0-n1
1566  do n3=0,r-2*n0-n1-n2
1567  do n4=0,r-2*n0-n1-n2-n3
1568  do n5=0,r-2*n0-n1-n2-n3-n4
1569  n6 = r-2*n0-n1-n2-n3-n4-n5
1570 
1571  cnt = cnt+1
1572  fct(cnt) = guvaux(n0,n1,n2,n3,n4,n5,n6)
1573 
1574  end do
1575  end do
1576  end do
1577  end do
1578  end do
1579  end do
1580  cnt = cnt+1
1581  fct(cnt) = gerraux(r)
1582  cnt = cnt+1
1583  fct(cnt) = gerr2aux(r)
1584  end do
1585 
1586  if (rmax.ge.10) then
1587  call writecache(fct,ncoefs(rank,7)+ncoefs(rank-10,7)+2*(rank+1),id,7,rank)
1588  else
1589  call writecache(fct,ncoefs(rank,7)+2*(rank+1),id,7,rank)
1590  end if
1591 
1592  end if
1593 
1594  g = gaux(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
1595  guv = guvaux(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
1596  gerr = gerraux(0:rmax)
1597  gerr2 = gerr2aux(0:rmax)
1598  return
1599 
1600  end if
1601 
1602 ! end if
1603  end if
1604  end if
1605 
1606 
1607  call calcgred(g,guv,p10,p21,p32,p43,p54,p65,p60,p20,p31,p42,p53, &
1608  p64,p50,p61,p30,p41,p52,p63,p40,p51,p62, &
1609  m02,m12,m22,m32,m42,m52,m62,rmax,id,gerr,gerr2)
1610 
1611 
1612  end subroutine calcg
1613 
1614 
1615 
1616 
1617 
1618  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1619  ! subroutine CalcGred(G,Guv,p10,p21,p32,p43,p54,p65,p60,p20,p31,p42,p53, &
1620  ! p64,p50,p61,p30,p41,p52,p63,p40,p51,p62, &
1621  ! m02,m12,m22,m32,m42,m52,m62,rmax,id,Gerr,Gerr2)
1622  !
1623  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1624 
1625  subroutine calcgred(G,Guv,p10,p21,p32,p43,p54,p65,p60,p20,p31,p42,p53, &
1626  p64,p50,p61,p30,p41,p52,p63,p40,p51,p62, &
1627  m02,m12,m22,m32,m42,m52,m62,rmax,id,Gerr,Gerr2)
1628 
1629  integer, intent(in) :: rmax,id
1630  double complex, intent(in) :: p10,p21,p32,p43,p54,p65,p60,p20,p31,p42,p53
1631  double complex, intent(in) :: p64,p50,p61,p30,p41,p52,p63,p40,p51,p62
1632  double complex, intent(in) :: m02,m12,m22,m32,m42,m52,m62
1633  double precision, intent(out) :: Gerr(0:rmax),Gerr2(0:rmax)
1634  double complex :: q10,q21,q32,q43,q54,q20,q31,q42,q53,q50,q30,q41,q52,q40,q51,q60
1635  double complex :: mm02,mm12,mm22,mm32,mm42,mm52,mm62
1636  double complex :: mx(0:5,0:5), mxinv(0:5,0:5),mx0k(5,5),mx0kinv(5,5), f(6)
1637  double complex :: det,newdet,mxinvs,mx0kinvs(5)
1638  double complex :: zmx0kinv(5,5),Z(5,5),zmx0kinvs(5)
1639  double precision :: maxZ,maxzmx0kinv(5),maxzmx0kinvs
1640  double complex, intent(out) :: G(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
1641  double complex, intent(out) :: Guv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax,0:rmax)
1642  double complex, allocatable :: F_0(:,:,:,:,:,:,:), F_i(:,:,:,:,:,:,:)
1643  double complex, allocatable :: Fuv_0(:,:,:,:,:,:,:), Fuv_i(:,:,:,:,:,:,:)
1644  double complex :: S(5), Gaux(5), elimminf2_coli,chdet,gram(1:5,1:5),gramdet
1645  double precision :: Ferr(0:5,0:rmax), Ferr2(0:5,0:rmax)
1646  integer :: r,n0,n1,n2,n3,n4,n5,n6,n,k,i,j,nid(0:5),r0,bin,kbest,rmaxF,rBCD,up
1647  logical :: errorwriteflag
1648  character(len=*),parameter :: fmt10 = "(A17,'(',d25.18,' , ',d25.18,' )')"
1649 
1650  r0=0
1651 
1652  ! allocation of F functions
1653  rmaxf = max(rmax-1,0)
1654  allocate(f_0(0:rmaxf/2,0:rmaxf,0:rmaxf,0:rmaxf,0:rmaxf,0:rmaxf,0:rmaxf))
1655  allocate(f_i(0:rmaxf/2,0:rmaxf,0:rmaxf,0:rmaxf,0:rmaxf,0:rmaxf,5))
1656  allocate(fuv_0(0:rmaxf/2,0:rmaxf,0:rmaxf,0:rmaxf,0:rmaxf,0:rmaxf,0:rmaxf))
1657  allocate(fuv_i(0:rmaxf/2,0:rmaxf,0:rmaxf,0:rmaxf,0:rmaxf,0:rmaxf,5))
1658 
1659 
1660  ! determine binaries for F-coefficients
1661  k=0
1662  bin = 1
1663  do while (k.le.5)
1664  if (mod(id/bin,2).eq.0) then
1665  nid(k) = id+bin
1666  k = k+1
1667  end if
1668  bin = 2*bin
1669  end do
1670 
1671  call calcf(f_0(:,0,:,:,:,:,:),fuv_0(:,0,:,:,:,:,:),p21,p32,p43,p54,p65,p61, &
1672  p31,p42,p53,p64,p51,p62,p41,p52,p63,m12,m22,m32,m42,m52,m62, &
1673  rmaxf,nid(0),ferr=ferr(0,0:rmax),ferr2=ferr2(0,0:rmax))
1674  call calcf(f_i(:,:,:,:,:,:,1),fuv_i(:,:,:,:,:,:,1),p20,p32,p43,p54,p65,p60, &
1675  p30,p42,p53,p64,p50,p62,p40,p52,p63,m02,m22,m32,m42,m52,m62, &
1676  rmaxf,nid(1),ferr=ferr(1,0:rmax),ferr2=ferr2(1,0:rmax))
1677  call calcf(f_i(:,:,:,:,:,:,2),fuv_i(:,:,:,:,:,:,2),p10,p31,p43,p54,p65,p60, &
1678  p30,p41,p53,p64,p50,p61,p40,p51,p63,m02,m12,m32,m42,m52,m62, &
1679  rmaxf,nid(2),ferr=ferr(2,0:rmax),ferr2=ferr(2,0:rmax))
1680  call calcf(f_i(:,:,:,:,:,:,3),fuv_i(:,:,:,:,:,:,3),p10,p21,p42,p54,p65,p60, &
1681  p20,p41,p52,p64,p50,p61,p40,p51,p62,m02,m12,m22,m42,m52,m62, &
1682  rmaxf,nid(3),ferr=ferr(3,0:rmax),ferr2=ferr2(3,0:rmax))
1683  call calcf(f_i(:,:,:,:,:,:,4),fuv_i(:,:,:,:,:,:,4),p10,p21,p32,p53,p65,p60, &
1684  p20,p31,p52,p63,p50,p61,p30,p51,p62,m02,m12,m22,m32,m52,m62, &
1685  rmaxf,nid(4),ferr=ferr(4,0:rmax),ferr2=ferr2(4,0:rmax))
1686  call calcf(f_i(:,:,:,:,:,:,5),fuv_i(:,:,:,:,:,:,5),p10,p21,p32,p43,p64,p60, &
1687  p20,p31,p42,p63,p40,p61,p30,p41,p62,m02,m12,m22,m32,m42,m62, &
1688  rmaxf,nid(5),ferr=ferr(5,0:rmax),ferr2=ferr2(5,0:rmax))
1689 
1690  ! shift of integration momentum in F\{0}
1691  do n1=1,rmaxf
1692  do n2=0,rmaxf-n1
1693  do n3=0,rmaxf-n1-n2
1694  do n4=0,rmaxf-n1-n2-n3
1695  do n5=0,rmaxf-n1-n2-n3-n4
1696  do n6=0,rmaxf-n1-n2-n3-n4-n5
1697  n0 = (rmaxf-n1-n2-n3-n4-n5-n6)/2
1698  fuv_0(0:n0,n1,n2,n3,n4,n5,n6) = -fuv_0(0:n0,n1-1,n2,n3,n4,n5,n6) &
1699  -fuv_0(0:n0,n1-1,n2+1,n3,n4,n5,n6) &
1700  -fuv_0(0:n0,n1-1,n2,n3+1,n4,n5,n6) &
1701  -fuv_0(0:n0,n1-1,n2,n3,n4+1,n5,n6) &
1702  -fuv_0(0:n0,n1-1,n2,n3,n4,n5+1,n6) &
1703  -fuv_0(0:n0,n1-1,n2,n3,n4,n5,n6+1)
1704  f_0(0:n0,n1,n2,n3,n4,n5,n6) = -f_0(0:n0,n1-1,n2,n3,n4,n5,n6) &
1705  -f_0(0:n0,n1-1,n2+1,n3,n4,n5,n6) &
1706  -f_0(0:n0,n1-1,n2,n3+1,n4,n5,n6) &
1707  -f_0(0:n0,n1-1,n2,n3,n4+1,n5,n6) &
1708  -f_0(0:n0,n1-1,n2,n3,n4,n5+1,n6) &
1709  -f_0(0:n0,n1-1,n2,n3,n4,n5,n6+1)
1710  end do
1711  end do
1712  end do
1713  end do
1714  end do
1715  end do
1716 
1717 
1718  ! determine inverse modified Caley matrix
1719  mm02 = elimminf2_coli(m02)
1720  mm12 = elimminf2_coli(m12)
1721  mm22 = elimminf2_coli(m22)
1722  mm32 = elimminf2_coli(m32)
1723  mm42 = elimminf2_coli(m42)
1724  mm52 = elimminf2_coli(m52)
1725  mm62 = elimminf2_coli(m62)
1726  q10 = elimminf2_coli(p10)
1727  q21 = elimminf2_coli(p21)
1728  q32 = elimminf2_coli(p32)
1729  q43 = elimminf2_coli(p43)
1730  q54 = elimminf2_coli(p54)
1731  q50 = elimminf2_coli(p50)
1732  q20 = elimminf2_coli(p20)
1733  q31 = elimminf2_coli(p31)
1734  q42 = elimminf2_coli(p42)
1735  q53 = elimminf2_coli(p53)
1736  q40 = elimminf2_coli(p40)
1737  q51 = elimminf2_coli(p51)
1738  q30 = elimminf2_coli(p30)
1739  q41 = elimminf2_coli(p41)
1740  q52 = elimminf2_coli(p52)
1741  q60 = elimminf2_coli(p60)
1742 
1743  f(1) = q10+mm02-mm12
1744  f(2) = q20+mm02-mm22
1745  f(3) = q30+mm02-mm32
1746  f(4) = q40+mm02-mm42
1747  f(5) = q50+mm02-mm52
1748  f(6) = q60+mm02-mm62
1749 
1750  mx(0,0) = 2d0*mm02
1751  mx(1,0) = q10 - mm12 + mm02
1752  mx(2,0) = q20 - mm22 + mm02
1753  mx(3,0) = q30 - mm32 + mm02
1754  mx(4,0) = q40 - mm42 + mm02
1755  mx(5,0) = q50 - mm52 + mm02
1756  mx(0,1) = mx(1,0)
1757  mx(1,1) = 2d0*q10
1758  mx(2,1) = q10+q20-q21
1759  mx(3,1) = q10+q30-q31
1760  mx(4,1) = q10+q40-q41
1761  mx(5,1) = q10+q50-q51
1762  mx(0,2) = mx(2,0)
1763  mx(1,2) = mx(2,1)
1764  mx(2,2) = 2d0*q20
1765  mx(3,2) = q20+q30-q32
1766  mx(4,2) = q20+q40-q42
1767  mx(5,2) = q20+q50-q52
1768  mx(0,3) = mx(3,0)
1769  mx(1,3) = mx(3,1)
1770  mx(2,3) = mx(3,2)
1771  mx(3,3) = 2d0*q30
1772  mx(4,3) = q30+q40-q43
1773  mx(5,3) = q30+q50-q53
1774  mx(0,4) = mx(4,0)
1775  mx(1,4) = mx(4,1)
1776  mx(2,4) = mx(4,2)
1777  mx(3,4) = mx(4,3)
1778  mx(4,4) = 2d0*q40
1779  mx(5,4) = q40+q50-q54
1780  mx(0,5) = mx(5,0)
1781  mx(1,5) = mx(5,1)
1782  mx(2,5) = mx(5,2)
1783  mx(3,5) = mx(5,3)
1784  mx(4,5) = mx(5,4)
1785  mx(5,5) = 2d0*q50
1786 
1787  call chinv(6,mx,mxinv)
1788 
1789  ! determine X_(0,5)
1790  do j=1,5
1791  do i=1,5
1792  mx0k(i,j) = mx(i,j-1)
1793  end do
1794  end do
1795 
1796  det = chdet(5,mx0k)
1797  kbest = 5
1798 
1799  do j=5,2,-1
1800  do i=1,5
1801  mx0k(i,j) = mx(i,j)
1802  end do
1803 
1804  newdet = chdet(5,mx0k)
1805  if (abs(newdet).gt.abs(det)) then
1806  kbest = j-1
1807  det = newdet
1808  end if
1809 
1810  end do
1811 
1812  do i=1,5
1813  mx0k(i,1) = mx(i,1)
1814  mx0k(i,kbest) = mx(i,0)
1815  end do
1816 
1817  call chinv(5,mx0k,mx0kinv)
1818  do i=1,5
1819  mx0kinv(kbest,i) = 0d0
1820  end do
1821 
1822  mxinvs = sum(mxinv(0:5,0))
1823  do i=1,5
1824  mx0kinvs(i) = sum(mx0kinv(i,1:5))
1825  end do
1826 
1827  ! for alternative error estimate
1828  z(1:5,1:5) = mx(1:5,1:5)
1829 
1830  zmx0kinv = matmul(z,mx0kinv)
1831 
1832  do i=1,5
1833  maxzmx0kinv(i) = maxval(abs(zmx0kinv(1:5,i)))
1834  zmx0kinvs(i) = sum(zmx0kinv(i,1:5))
1835  end do
1836 
1837  maxzmx0kinvs = maxval(abs(zmx0kinvs(1:5)))
1838 
1839  maxz = maxval(abs(z))
1840 
1841 
1842  ! calculation of UV-divergent parts
1843  ! G_(n0,n1,n2,n3) UV-finite for n0<5
1844  guv(0:min(rmax/2,4),:,:,:,:,:,:) = 0d0
1845 
1846  ! PV reduction (5.10)
1847  do r=max(r0,10),rmax
1848  do n0=5,r/2
1849  do n1=0,r-2*n0
1850  do n2=0,r-2*n0-n1
1851  do n3=0,r-2*n0-n1-n2
1852  do n4=0,r-2*n0-n1-n2-n3
1853  do n5=0,r-2*n0-n1-n2-n3-n4
1854  n6 = r-2*n0-n1-n2-n3-n4-n5
1855 
1856  guv(n0,n1,n2,n3,n4,n5,n6) = (fuv_0(n0-1,n1,n2,n3,n4,n5,n6) &
1857  + 2*m02*guv(n0-1,n1,n2,n3,n4,n5,n6) &
1858  + f(1)*guv(n0-1,n1+1,n2,n3,n4,n5,n6) &
1859  + f(2)*guv(n0-1,n1,n2+1,n3,n4,n5,n6) &
1860  + f(3)*guv(n0-1,n1,n2,n3+1,n4,n5,n6) &
1861  + f(4)*guv(n0-1,n1,n2,n3,n4+1,n5,n6) &
1862  + f(5)*guv(n0-1,n1,n2,n3,n4,n5+1,n6) &
1863  + f(6)*guv(n0-1,n1,n2,n3,n4,n5,n6+1))/ (2*(r-4))
1864 
1865  end do
1866  end do
1867  end do
1868  end do
1869  end do
1870  end do
1871  end do
1872 
1873 
1874  ! scalar coefficient
1875  if (r0.eq.0) then
1876  g = 0d0
1877  g(0,0,0,0,0,0,0) = -mxinv(0,0)*f_0(0,0,0,0,0,0,0)
1878  do k=1,5
1879  g(0,0,0,0,0,0,0) = g(0,0,0,0,0,0,0) &
1880  + mxinv(k,0)*(f_i(0,0,0,0,0,0,k)-f_0(0,0,0,0,0,0,0))
1881  end do
1882  end if
1883 
1884 ! Gerr(0) = max( maxval(abs(mxinv(0:5,0)))*Ferr(0,0), &
1885  gerr(0) = max( abs(mxinvs)*ferr(0,0), &
1886  abs(mxinv(1,0))*ferr(1,0) , &
1887  abs(mxinv(2,0))*ferr(2,0) , &
1888  abs(mxinv(3,0))*ferr(3,0) , &
1889  abs(mxinv(4,0))*ferr(4,0) , &
1890  abs(mxinv(5,0))*ferr(5,0) )
1891  gerr2(0) = max(abs(mxinvs)*ferr2(0,0), &
1892  abs(mxinv(1,0))*ferr2(1,0) , &
1893  abs(mxinv(2,0))*ferr2(2,0) , &
1894  abs(mxinv(3,0))*ferr2(3,0) , &
1895  abs(mxinv(4,0))*ferr2(4,0) , &
1896  abs(mxinv(5,0))*ferr2(5,0) )
1897 
1898  ! formula (7.13) extended to N=7
1899  do r=r0,rmax-1
1900  do n0=0,max((r-1)/2,0)
1901  do n1=0,r-2*n0
1902  do n2=0,r-2*n0-n1
1903  do n3=0,r-2*n0-n1-n2
1904  do n4=0,r-2*n0-n1-n2-n3
1905  do n5=0,r-2*n0-n1-n2-n3-n4
1906  n6 = r-2*n0-n1-n2-n3-n4-n5
1907 
1908  do n=1,5
1909  s(n) = -f_0(n0,n1,n2,n3,n4,n5,n6)
1910  end do
1911 
1912  if (n1.eq.0) then
1913  s(1) = s(1) + f_i(n0,n2,n3,n4,n5,n6,1)
1914  end if
1915  if (n2.eq.0) then
1916  s(2) = s(2) + f_i(n0,n1,n3,n4,n5,n6,2)
1917  end if
1918  if (n3.eq.0) then
1919  s(3) = s(3) + f_i(n0,n1,n2,n4,n5,n6,3)
1920  end if
1921  if (n4.eq.0) then
1922  s(4) = s(4) + f_i(n0,n1,n2,n3,n5,n6,4)
1923  end if
1924  if (n5.eq.0) then
1925  s(5) = s(5) + f_i(n0,n1,n2,n3,n4,n6,5)
1926  end if
1927 
1928  do k=1,5
1929  gaux(k) = mx0kinv(k,1)*s(1)+mx0kinv(k,2)*s(2) &
1930  + mx0kinv(k,3)*s(3)+mx0kinv(k,4)*s(4)+mx0kinv(k,5)*s(5)
1931  end do
1932 
1933  g(n0,n1+1,n2,n3,n4,n5,n6) = g(n0,n1+1,n2,n3,n4,n5,n6) + (n1+1)*gaux(1)/(r+1)
1934  g(n0,n1,n2+1,n3,n4,n5,n6) = g(n0,n1,n2+1,n3,n4,n5,n6) + (n2+1)*gaux(2)/(r+1)
1935  g(n0,n1,n2,n3+1,n4,n5,n6) = g(n0,n1,n2,n3+1,n4,n5,n6) + (n3+1)*gaux(3)/(r+1)
1936  g(n0,n1,n2,n3,n4+1,n5,n6) = g(n0,n1,n2,n3,n4+1,n5,n6) + (n4+1)*gaux(4)/(r+1)
1937  g(n0,n1,n2,n3,n4,n5+1,n6) = g(n0,n1,n2,n3,n4,n5+1,n6) + (n5+1)*gaux(5)/(r+1)
1938 
1939  end do
1940  end do
1941  end do
1942  end do
1943  end do
1944  end do
1945 
1946  if (r.le.rmax-1) then
1947 ! Gerr(r+1) = max( maxval(abs(mx0kinv(1:5,1:5)))*Ferr(0,r), &
1948  gerr(r+1) = max( maxval(abs(mx0kinvs(1:5)))*ferr(0,r), &
1949  maxval(abs(mx0kinv(1:5,1)))*ferr(1,r) , &
1950  maxval(abs(mx0kinv(1:5,2)))*ferr(2,r) , &
1951  maxval(abs(mx0kinv(1:5,3)))*ferr(3,r) , &
1952  maxval(abs(mx0kinv(1:5,4)))*ferr(4,r) , &
1953  maxval(abs(mx0kinv(1:5,5)))*ferr(5,r) )
1954  gerr2(r+1) = max( abs(maxzmx0kinvs)*ferr2(0,r), &
1955  abs(maxzmx0kinv(1))*ferr2(1,r) , &
1956  abs(maxzmx0kinv(2))*ferr2(2,r) , &
1957  abs(maxzmx0kinv(3))*ferr2(3,r) , &
1958  abs(maxzmx0kinv(4))*ferr2(4,r) , &
1959  abs(maxzmx0kinv(5))*ferr2(5,r) ) /maxz
1960  end if
1961 
1962  if (mode_coli.lt.1) then
1963  gram= mx(1:5,1:5)
1964 
1965  gramdet= chdet(5,gram)
1966 
1967  if (max(abs(g(0,0,0,0,0,0,0)),abs(g(0,1,0,0,0,0,0)),abs(g(0,0,1,0,0,0,0)), &
1968  abs(g(0,0,0,1,0,0,0)),abs(g(0,0,0,0,1,0,0)),abs(g(0,0,0,0,0,1,0)), &
1969  abs(g(0,0,0,0,0,0,1)))*abs(gramdet/det).gt. gerr(r+1)) then
1970 
1971 
1972 ! write(*,*) 'CalcGred Gerr=',r+1,Gerr(r+1)
1973 
1974  gerr(r+1)=max(gerr(r+1),max( &
1975  abs(g(0,0,0,0,0,0,0)),abs(g(0,1,0,0,0,0,0)), &
1976  abs(g(0,0,1,0,0,0,0)), abs(g(0,0,0,1,0,0,0)), &
1977  abs(g(0,0,0,0,1,0,0)),abs(g(0,0,0,0,0,1,0)), &
1978  abs(g(0,0,0,0,0,0,1)))*abs(gramdet/det) )
1979 
1980  gerr2(r+1)=max(gerr2(r+1),max( &
1981  abs(g(0,0,0,0,0,0,0)),abs(g(0,1,0,0,0,0,0)), &
1982  abs(g(0,0,1,0,0,0,0)), abs(g(0,0,0,1,0,0,0)), &
1983  abs(g(0,0,0,0,1,0,0)),abs(g(0,0,0,0,0,1,0)), &
1984  abs(g(0,0,0,0,0,0,1)))*abs(gramdet/det) )
1985 
1986 ! write(*,*) 'CalcGred ',abs(G(0,0,0,0,0,0)),abs(G(0,1,0,0,0,0)), &
1987 ! abs(G(0,0,1,0,0,0)), abs(G(0,0,0,1,0,0)), &
1988 ! abs(G(0,0,0,0,1,0)),abs(G(0,0,0,0,0,1)),abs(gramdet/det)
1989 ! write(*,*) 'CalcGred Gerr=',r+1,Gerr(r+1)
1990 
1991  if (abs(gramdet/det).gt.reqacc_coli) then
1992  call seterrflag_coli(-6)
1993  call errout_coli('CalcGred', &
1994  'input momenta inconsistent! (not 4-dimensional)', &
1995  errorwriteflag)
1996  if (errorwriteflag) then
1997  write(nerrout_coli,fmt10) ' CalcGred: q10 = ',q10
1998  write(nerrout_coli,fmt10) ' CalcGred: q21 = ',q21
1999  write(nerrout_coli,fmt10) ' CalcGred: q32 = ',q32
2000  write(nerrout_coli,fmt10) ' CalcGred: q43 = ',q43
2001  write(nerrout_coli,fmt10) ' CalcGred: q54 = ',q54
2002  write(nerrout_coli,fmt10) ' CalcGred: q50 = ',q50
2003  write(nerrout_coli,fmt10) ' CalcGred: q20 = ',q10
2004  write(nerrout_coli,fmt10) ' CalcGred: q31 = ',q31
2005  write(nerrout_coli,fmt10) ' CalcGred: q42 = ',q42
2006  write(nerrout_coli,fmt10) ' CalcGred: q53 = ',q53
2007  write(nerrout_coli,fmt10) ' CalcGred: q40 = ',q40
2008  write(nerrout_coli,fmt10) ' CalcGred: q51 = ',q51
2009  write(nerrout_coli,fmt10) ' CalcGred: q30 = ',q30
2010  write(nerrout_coli,fmt10) ' CalcGred: q41 = ',q41
2011  write(nerrout_coli,fmt10) ' CalcGred: q52 = ',q52
2012  write(nerrout_coli,fmt10) ' CalcGred: mm02 = ',mm02
2013  write(nerrout_coli,fmt10) ' CalcGred: mm12 = ',mm12
2014  write(nerrout_coli,fmt10) ' CalcGred: mm22 = ',mm22
2015  write(nerrout_coli,fmt10) ' CalcGred: mm32 = ',mm32
2016  write(nerrout_coli,fmt10) ' CalcGred: mm42 = ',mm42
2017  write(nerrout_coli,fmt10) ' CalcGred: mm52 = ',mm52
2018  write(nerrout_coli,fmt10) ' CalcGred: gram = ',gramdet/det
2019  end if
2020  end if
2021 
2022  end if
2023  end if
2024 
2025  end do
2026 
2027  end subroutine calcgred
2028 
2029 
2030 
2031 end module reductionefg
reductiond
Definition: reductionD.F90:71
reductionefg::calcg
subroutine calcg(G, Guv, p10, p21, p32, p43, p54, p65, p60, p20, p31, p42, p53, p64, p50, p61, p30, p41, p52, p63, p40, p51, p62, m02, m12, m22, m32, m42, m52, m62, rmax, id, Gerr, Gerr2)
Definition: reductionEFG.F90:1361
reductionefg::calcfred
subroutine calcfred(F, Fuv, p10, p21, p32, p43, p54, p50, p20, p31, p42, p53, p40, p51, p30, p41, p52, m02, m12, m22, m32, m42, m52, rmax, id, Ferr, Ferr2)
Definition: reductionEFG.F90:901
reductionefg::calcered
subroutine calcered(E, Euv, p10, p21, p32, p43, p40, p20, p31, p42, p30, p41, m02, m12, m22, m32, m42, rmax, id, Eerr, Eerr2)
Definition: reductionEFG.F90:265
reductionefg::calcf
subroutine calcf(F, Fuv, p10, p21, p32, p43, p54, p50, p20, p31, p42, p53, p40, p51, p30, p41, p52, m02, m12, m22, m32, m42, m52, rmax, id, Ferr, Ferr2)
Definition: reductionEFG.F90:657
reductionefg::calce
subroutine calce(E, Euv, p10, p21, p32, p43, p40, p20, p31, p42, p30, p41, m02, m12, m22, m32, m42, rmax, id, Eerr, Eerr2)
Definition: reductionEFG.F90:43
reductiond::calcd
subroutine calcd(D, Duv, p10, p21, p32, p30, p20, p31, m02, m12, m22, m32, rmax, id, Derr1, Derr2)
Definition: reductionD.F90:94
reductionefg::calcgred
subroutine calcgred(G, Guv, p10, p21, p32, p43, p54, p65, p60, p20, p31, p42, p53, p64, p50, p61, p30, p41, p52, p63, p40, p51, p62, m02, m12, m22, m32, m42, m52, m62, rmax, id, Gerr, Gerr2)
Definition: reductionEFG.F90:1628
reductionefg
Definition: reductionEFG.F90:25
Generated by   doxygen 1.8.17