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)
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
57 if ((use_cache_system).and.(tencache.gt.5))
then
58 if ((ncache.gt.0).and.(ncache.le.ncache_max))
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)
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)
88 euv(0:min(rmax/2,2),:,:,:,:) = 0d0
97 e(n0,n1,n2,n3,n4) = fct(cnt)
110 euv(n0,n1,n2,n3,n4) = fct(cnt)
117 eerr(r) = real(fct(cnt))
119 eerr2(r) = real(fct(cnt))
125 if(rank.eq.rmax)
then
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)
140 fct(cnt) = e(n0,n1,n2,n3,n4)
153 fct(cnt) = euv(n0,n1,n2,n3,n4)
166 call writecache(fct,ncoefs(rank,5)+ncoefs(rank-6,5)+2*(rank+1),id,5,rank)
168 call writecache(fct,ncoefs(rank,5)+2*(rank+1),id,5,rank)
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))
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)
189 allocate(fct(ncoefs(rank,5)+ncoefs(rank-6,5)+2*(rank+1)))
192 allocate(fct(ncoefs(rank,5)+2*(rank+1)))
202 fct(cnt) = eaux(n0,n1,n2,n3,n4)
215 fct(cnt) = euvaux(n0,n1,n2,n3,n4)
222 fct(cnt) = eerraux(r)
224 fct(cnt) = eerr2aux(r)
228 call writecache(fct,ncoefs(rank,5)+ncoefs(rank-6,5)+2*(rank+1),id,5,rank)
230 call writecache(fct,ncoefs(rank,5)+2*(rank+1),id,5,rank)
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)
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)
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)
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(:,:,:,:,:)
284 double complex :: S(0:4), S2(4,4), Eaux(4), Eaux2, Eaux3, elimminf2_coli
285 double precision,
allocatable :: Derr(:,:),Derr2(:,:)
287 integer :: r,n0,n1,n2,n3,n4,i,n,k,nid(0:4),rid,r0,bin,rBCD,rmaxD
292 rmaxd = max(rmax-1,0)
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))
307 if (mod(id/bin,2).eq.0)
then
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))
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)
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)
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
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
389 call chinv(5,mx,mxinv)
392 mxinvs(i) = sum(mxinv(i,0:4))
396 z(1:4,1:4) = mx(1:4,1:4)
398 zmxinv = matmul(z,mxinv(1:4,0:4))
401 maxzmxinv(i) = maxval(abs(zmxinv(1:4,i)))
405 zmxinvs(i) = sum(zmxinv(i,0:4))
408 maxzmxinvs = maxval(abs(zmxinvs(1:4)))
410 maxz = maxval(abs(z))
415 euv_aux(0:min((rmax+1)/2,2),:,:,:,:) = 0d0
416 euv(0:min(rmax/2,2),:,:,:,:) = 0d0
419 do r=max(r0,6),rmax+1
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))
433 euv(n0,n1,n2,n3,n4) = euv_aux(n0,n1,n2,n3,n4)
454 e(0,0,0,0,0) = -mxinvs(0)*d_0(0,0,0,0,0)
456 e(0,0,0,0,0) = e(0,0,0,0,0) + mxinv(k,0)*d_i(0,0,0,0,k)
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) )
485 if (n0.gt.0.or.r.le.rmax-1)
then
493 if (n0.gt.0.or.r.le.rmax-1)
then
494 s(1) = s(1) + d_i(n0,n2,n3,n4,1)
498 if (r.le.rmax-1)
then
499 s2(:,1) = -n1*d_0(n0+1,n1-1,n2,n3,n4)
501 s2(2,1) = s2(2,1) + n1*d_i(n0+1,n1-1,n3,n4,2)
504 s2(3,1) = s2(3,1) + n1*d_i(n0+1,n1-1,n2,n4,3)
507 s2(4,1) = s2(4,1) + n1*d_i(n0+1,n1-1,n2,n3,4)
513 if (n0.gt.0.or.r.le.rmax-1)
then
514 s(2) = s(2) + d_i(n0,n1,n3,n4,2)
518 if (r.le.rmax-1)
then
519 s2(:,2) = -n2*d_0(n0+1,n1,n2-1,n3,n4)
521 s2(1,2) = s2(1,2) + n2*d_i(n0+1,n2-1,n3,n4,1)
524 s2(3,2) = s2(3,2) + n2*d_i(n0+1,n1,n2-1,n4,3)
527 s2(4,2) = s2(4,2) + n2*d_i(n0+1,n1,n2-1,n3,4)
533 if (n0.gt.0.or.r.le.rmax-1)
then
534 s(3) = s(3) + d_i(n0,n1,n2,n4,3)
538 if (r.le.rmax-1)
then
539 s2(:,3) = -n3*d_0(n0+1,n1,n2,n3-1,n4)
541 s2(1,3) = s2(1,3) + n3*d_i(n0+1,n2,n3-1,n4,1)
544 s2(2,3) = s2(2,3) + n3*d_i(n0+1,n1,n3-1,n4,2)
547 s2(4,3) = s2(4,3) + n3*d_i(n0+1,n1,n2,n3-1,4)
553 if (n0.gt.0.or.r.le.rmax-1)
then
554 s(4) = s(4) + d_i(n0,n1,n2,n3,4)
558 if (r.le.rmax-1)
then
559 s2(:,4) = -n4*d_0(n0+1,n1,n2,n3,n4-1)
561 s2(1,4) = s2(1,4) + n4*d_i(n0+1,n2,n3,n4-1,1)
564 s2(2,4) = s2(2,4) + n4*d_i(n0+1,n1,n3,n4-1,2)
567 s2(3,4) = s2(3,4) + n4*d_i(n0+1,n1,n2,n4-1,3)
572 if (r.le.rmax-1)
then
575 s(0) = s(0) - 4d0*euv_aux(n0+1,n1,n2,n3,n4)
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
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)
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)
601 e(n0,n1,n2,n3,n4) = e(n0,n1,n2,n3,n4) + 2*n0*eaux2/r
609 if (r.le.rmax-1)
then
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
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, &
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
672 if ((use_cache_system).and.(tencache.gt.6))
then
673 if ((ncache.gt.0).and.(ncache.le.ncache_max))
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)
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)
709 fuv(0:min(rmax/2,3),:,:,:,:,:) = 0d0
715 do n4=0,r-2*n0-n1-n2-n3
716 n5 = r-2*n0-n1-n2-n3-n4
719 f(n0,n1,n2,n3,n4,n5) = fct(cnt)
730 do n4=0,r-2*n0-n1-n2-n3
731 n5 = r-2*n0-n1-n2-n3-n4
734 fuv(n0,n1,n2,n3,n4,n5) = fct(cnt)
742 ferr(r) = real(fct(cnt))
744 ferr2(r) = real(fct(cnt))
750 if(rank.eq.rmax)
then
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)
762 do n4=0,r-2*n0-n1-n2-n3
763 n5 = r-2*n0-n1-n2-n3-n4
766 fct(cnt) = f(n0,n1,n2,n3,n4,n5)
777 do n4=0,r-2*n0-n1-n2-n3
778 n5 = r-2*n0-n1-n2-n3-n4
781 fct(cnt) = fuv(n0,n1,n2,n3,n4,n5)
795 call writecache(fct,ncoefs(rank,6)+ncoefs(rank-8,6)+2*(rank+1),id,6,rank)
797 call writecache(fct,ncoefs(rank,6)+2*(rank+1),id,6,rank)
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))
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)
819 allocate(fct(ncoefs(rank,6)+ncoefs(rank-8,6)+2*(rank+1)))
821 allocate(fct(ncoefs(rank,6)+2*(rank+1)))
828 do n4=0,r-2*n0-n1-n2-n3
829 n5 = r-2*n0-n1-n2-n3-n4
832 fct(cnt) = faux(n0,n1,n2,n3,n4,n5)
843 do n4=0,r-2*n0-n1-n2-n3
844 n5 = r-2*n0-n1-n2-n3-n4
847 fct(cnt) = fuvaux(n0,n1,n2,n3,n4,n5)
855 fct(cnt) = ferraux(r)
857 fct(cnt) = ferr2aux(r)
861 call writecache(fct,ncoefs(rank,6)+ncoefs(rank-8,6)+2*(rank+1),id,6,rank)
863 call writecache(fct,ncoefs(rank,6)+2*(rank+1),id,6,rank)
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)
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)
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, &
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,' )')"
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))
938 if (mod(id/bin,2).eq.0)
then
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))
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)
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)
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)
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
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
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
1035 mx(3,2) = q20+q30-q32
1036 mx(4,2) = q20+q40-q42
1037 mx(5,2) = q20+q50-q52
1042 mx(4,3) = q30+q40-q43
1043 mx(5,3) = q30+q50-q53
1049 mx(5,4) = q40+q50-q54
1057 call chinv(6,mx,mxinv)
1062 mx0k(i,j) = mx(i,j-1)
1077 newdet = chdet(5,mx0k)
1081 if (abs(newdet).gt.abs(det))
then
1092 mx0k(i,kbest) = mx(i,0)
1095 call chinv(5,mx0k,mx0kinv)
1097 mx0kinv(kbest,i) = 0d0
1100 mxinvs = sum(mxinv(0:5,0))
1102 mx0kinvs(i) = sum(mx0kinv(i,1:5))
1106 z(1:5,1:5) = mx(1:5,1:5)
1108 zmx0kinv = matmul(z,mx0kinv)
1111 maxzmx0kinv(i) = maxval(abs(zmx0kinv(1:5,i)))
1112 zmx0kinvs(i) = sum(zmx0kinv(i,1:5))
1115 maxzmx0kinvs = maxval(abs(zmx0kinvs(1:5)))
1117 maxz = maxval(abs(z))
1123 fuv(0:min(rmax/2,3),:,:,:,:,:) = 0d0
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
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))
1152 f(0,0,0,0,0,0) = -mxinv(0,0)*e_0(0,0,0,0,0,0)
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))
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) )
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
1197 s(n) = -e_0(n0,n1,n2,n3,n4,n5)
1201 s(1) = s(1) + e_i(n0,n2,n3,n4,n5,1)
1204 s(2) = s(2) + e_i(n0,n1,n3,n4,n5,2)
1207 s(3) = s(3) + e_i(n0,n1,n2,n4,n5,3)
1210 s(4) = s(4) + e_i(n0,n1,n2,n3,n5,4)
1213 s(5) = s(5) + e_i(n0,n1,n2,n3,n4,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)
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)
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
1272 if (mode_coli.lt.1)
then
1275 gramdet= chdet(5,gram)
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
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) )
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) )
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)', &
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
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)
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
1377 if ((use_cache_system).and.(tencache.gt.7))
then
1378 if ((ncache.gt.0).and.(ncache.le.ncache_max))
then
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)
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)
1421 guv(0:min(rmax/2,4),:,:,:,:,:,:) = 0d0
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
1432 g(n0,n1,n2,n3,n4,n5,n6) = fct(cnt)
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
1449 guv(n0,n1,n2,n3,n4,n5,n6) = fct(cnt)
1458 gerr(r) = real(fct(cnt))
1460 gerr2(r) = real(fct(cnt))
1466 if(rank.eq.rmax)
then
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)
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
1484 fct(cnt) = g(n0,n1,n2,n3,n4,n5,n6)
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
1501 fct(cnt) = guv(n0,n1,n2,n3,n4,n5,n6)
1515 if (rmax.ge.10)
then
1516 call writecache(fct,ncoefs(rank,7)+ncoefs(rank-10,7)+2*(rank+1),id,7,rank)
1518 call writecache(fct,ncoefs(rank,7)+2*(rank+1),id,7,rank)
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))
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)
1538 if (rmax.ge.10)
then
1540 allocate(fct(ncoefs(rank,7)+ncoefs(rank-10,7)+2*(rank+1)))
1543 allocate(fct(ncoefs(rank,7)+2*(rank+1)))
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
1555 fct(cnt) = gaux(n0,n1,n2,n3,n4,n5,n6)
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
1572 fct(cnt) = guvaux(n0,n1,n2,n3,n4,n5,n6)
1581 fct(cnt) = gerraux(r)
1583 fct(cnt) = gerr2aux(r)
1586 if (rmax.ge.10)
then
1587 call writecache(fct,ncoefs(rank,7)+ncoefs(rank-10,7)+2*(rank+1),id,7,rank)
1589 call writecache(fct,ncoefs(rank,7)+2*(rank+1),id,7,rank)
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)
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)
1612 end subroutine calcg
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)
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,' )')"
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))
1664 if (mod(id/bin,2).eq.0)
then
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))
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)
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)
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
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
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
1765 mx(3,2) = q20+q30-q32
1766 mx(4,2) = q20+q40-q42
1767 mx(5,2) = q20+q50-q52
1772 mx(4,3) = q30+q40-q43
1773 mx(5,3) = q30+q50-q53
1779 mx(5,4) = q40+q50-q54
1787 call chinv(6,mx,mxinv)
1792 mx0k(i,j) = mx(i,j-1)
1804 newdet = chdet(5,mx0k)
1805 if (abs(newdet).gt.abs(det))
then
1814 mx0k(i,kbest) = mx(i,0)
1817 call chinv(5,mx0k,mx0kinv)
1819 mx0kinv(kbest,i) = 0d0
1822 mxinvs = sum(mxinv(0:5,0))
1824 mx0kinvs(i) = sum(mx0kinv(i,1:5))
1828 z(1:5,1:5) = mx(1:5,1:5)
1830 zmx0kinv = matmul(z,mx0kinv)
1833 maxzmx0kinv(i) = maxval(abs(zmx0kinv(1:5,i)))
1834 zmx0kinvs(i) = sum(zmx0kinv(i,1:5))
1837 maxzmx0kinvs = maxval(abs(zmx0kinvs(1:5)))
1839 maxz = maxval(abs(z))
1844 guv(0:min(rmax/2,4),:,:,:,:,:,:) = 0d0
1847 do r=max(r0,10),rmax
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
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))
1877 g(0,0,0,0,0,0,0) = -mxinv(0,0)*f_0(0,0,0,0,0,0,0)
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))
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) )
1900 do n0=0,max((r-1)/2,0)
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
1909 s(n) = -f_0(n0,n1,n2,n3,n4,n5,n6)
1913 s(1) = s(1) + f_i(n0,n2,n3,n4,n5,n6,1)
1916 s(2) = s(2) + f_i(n0,n1,n3,n4,n5,n6,2)
1919 s(3) = s(3) + f_i(n0,n1,n2,n4,n5,n6,3)
1922 s(4) = s(4) + f_i(n0,n1,n2,n3,n5,n6,4)
1925 s(5) = s(5) + f_i(n0,n1,n2,n3,n4,n6,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)
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)
1946 if (r.le.rmax-1)
then
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
1962 if (mode_coli.lt.1)
then
1965 gramdet= chdet(5,gram)
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
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) )
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) )
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)', &
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