22 #define ALWAYSPV ! default
25 #define Cutrloop ! default
92 subroutine calcd(D,Duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32, &
95 integer,
intent(in) :: rmax,id
96 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
97 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
98 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
99 double precision,
intent(out) :: Derr1(0:rmax),Derr2(0:rmax)
100 double complex,
allocatable :: Daux(:,:,:,:), Duvaux(:,:,:,:), fct(:)
101 double precision,
allocatable :: Derr1aux(:),Derr2aux(:)
102 double complex :: x(10)
103 integer :: rank,switch,cnt,n0,n1,n2,n3,r
104 logical :: nocalc,wrica
109 if (use_cache_system)
then
110 if ((ncache.gt.0).and.(ncache.le.ncache_max))
then
126 allocate(fct(ncoefsg(rmax,4)-ncoefs(rmax-2,4)+ncoefs(rmax-3,4)+2*(rmax+1)))
127 call readcache(fct,ncoefsg(rmax,4)-ncoefs(rmax-2,4)+ncoefs(rmax-3,4)+2*(rmax+1),x,10,1,id,4,rank,nocalc,wrica)
128 else if(rmax.eq.2)
then
129 allocate(fct(ncoefsg(rmax,4)-1+2*(rmax+1)))
130 call readcache(fct,ncoefsg(rmax,4)-1+2*(rmax+1),x,10,1,id,4,rank,nocalc,wrica)
132 allocate(fct(ncoefsg(rmax,4)+2*(rmax+1)))
133 call readcache(fct,ncoefsg(rmax,4)+2*(rmax+1),x,10,1,id,4,rank,nocalc,wrica)
138 duv(0:min(rmax/2,1),:,:,:) = 0d0
144 d(0,n1,n2,n3) = fct(cnt)
153 d(n0,n1,n2,n3) = fct(cnt)
164 duv(n0,n1,n2,n3) = fct(cnt)
169 derr1(r) = real(fct(cnt))
171 derr2(r) = real(fct(cnt))
178 if(rank.eq.rmax)
then
180 call calcdred(d,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rank,id,derr1,derr2)
190 fct(cnt) = d(0,n1,n2,n3)
198 fct(cnt) = d(n0,n1,n2,n3)
207 fct(cnt) = duv(n0,n1,n2,n3)
218 call writecache(fct,ncoefsg(rank,4)-ncoefs(rank-2,4)+ncoefs(rank-3,4)+2*(rank+1),id,4,rank)
219 else if(rank.eq.2)
then
220 call writecache(fct,ncoefsg(rank,4)-1+2*(rank+1),id,4,rank)
222 call writecache(fct,ncoefsg(rank,4)+2*(rank+1),id,4,rank)
231 allocate(daux(0:rank,0:rank,0:rank,0:rank))
232 allocate(duvaux(0:rank,0:rank,0:rank,0:rank))
233 allocate(derr1aux(0:rank))
234 allocate(derr2aux(0:rank))
236 call calcdred(daux,duvaux,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rank,id,derr1aux,derr2aux)
242 allocate(fct(ncoefsg(rank,4)-ncoefs(rank-2,4)+ncoefs(rank-3,4)+2*(rank+1)))
243 else if(rank.eq.2)
then
244 allocate(fct(ncoefsg(rank,4)-1+2*(rank+1)))
246 allocate(fct(ncoefsg(rank,4)+2*(rank+1)))
256 fct(cnt) = daux(n0,n1,n2,n3)
267 fct(cnt) = duvaux(n0,n1,n2,n3)
273 fct(cnt) = derr1aux(r)
275 fct(cnt) = derr2aux(r)
279 call writecache(fct,ncoefsg(rank,4)-ncoefs(rank-2,4)+ncoefs(rank-3,4)+2*(rank+1),id,4,rank)
280 else if(rank.eq.2)
then
281 call writecache(fct,ncoefsg(rank,4)-1+2*(rank+1),id,4,rank)
283 call writecache(fct,ncoefsg(rank,4)+2*(rank+1),id,4,rank)
288 d = daux(0:rmax,0:rmax,0:rmax,0:rmax)
289 duv = duvaux(0:rmax,0:rmax,0:rmax,0:rmax)
290 derr1 = derr1aux(0:rmax)
291 derr2 = derr2aux(0:rmax)
305 call calcdred(d,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,derr1,derr2)
319 subroutine calcdred(D,Duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,Derr1,Derr2)
323 integer,
intent(in) :: rmax,id
324 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
327 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
328 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
329 double complex :: elimminf2_coli
330 double precision,
intent(out) :: Derr1(0:rmax),Derr2(0:rmax)
331 double precision :: D0est,Dtyp
333 double complex :: D0_coli
336 double complex :: chdet
337 integer :: rmaxC,r,rid,n0,n1,n2,n3,g,gy,gp,gr,gm,gpf,i,iexp
338 integer :: bin,k,nid(0:3)
339 logical :: use_pv,use_pv2,use_g,use_gy,use_gp,use_gr,use_gm,use_gpf
341 integer :: r_alt,Drmethod(0:rmax),DrCalc(0:rmax),DCalc
342 double complex,
allocatable :: C_i(:,:,:,:), Cuv_i(:,:,:,:)
343 double complex :: D_alt(0:rmax,0:rmax,0:rmax,0:rmax)
344 double complex :: Duv_alt(0:rmax,0:rmax,0:rmax,0:rmax)
345 double precision :: Derr(0:rmax),Derr_alt(0:rmax),Derr1_alt(0:rmax),Derr2_alt(0:rmax)
346 integer :: Drmethod_alt(0:rmax)
348 double precision :: err_pv(0:rmax),err_pv2(0:rmax),err_g(0:rmax),err_gy(0:rmax),err_gp(0:rmax)
349 double precision :: err_gr(0:rmax),err_gm(0:rmax),err_gpf(0:rmax)
350 double precision :: h_pv,w_pv,v_pv,z_pv,h_pv2,w_pv2,v_pv2,z_pv2,hw_pv2
351 double precision :: u_g,z_g,err_g_C,err_g_Cr,err_g_exp
352 double precision :: u_gm,z_gm,err_gm_C,err_gm_Cr,err_gm_exp
353 double precision :: v1_gy,b_gy,err_gy_C,err_gy_Cr,err_gy_exp
354 double precision :: v_gp,v1_gp,z_gp,err_gp_C,err_gp_Cr,err_gp_exp
355 double precision :: v1_gpf,b_gpf,err_gpf_C,err_gpf_Cr,err_gpf_exp
356 double precision :: x_gr,y_gr,y1_gr,a_gr,err_gr_C,err_gr_Cr,err_gr_exp
357 double precision :: err_C0,Cerr_i(0:rmax_C,0:3),err_C(0:rmax_C),err_D0,acc_D,errfac(0:3),err_req_D,err_inf,Cerr2_i(0:rmax_C,0:3)
358 double precision :: checkest,norm,Dscale2
359 logical :: lerr_D0,errorwriteflag
361 character(len=*),
parameter :: fmt1 =
"(A7,'dcmplx(',d25.18,' , ',d25.18,' )')"
362 character(len=*),
parameter :: fmt10 =
"(A17,'(',d25.18,' , ',d25.18,' )')"
364 #ifdef CritPointsCOLI
365 integer,
parameter :: MaxCritPointD=50
367 integer,
parameter :: MaxCritPointD=0
369 integer,
save :: CritPointCntD
372 data critpointcntd /0/
375 write(*,*)
'CalcDred in',rmax,id,acc_req_d
376 write(*,*)
'CalcDred acc_req',acc_req_d,reqacc_coli
377 write(*,*)
'CalcDred in',p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
388 rmaxc = max(rmax-1,3)
389 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,0:3))
390 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,0:3))
396 if (mod(id/bin,2).eq.0)
then
404 call calcc(c_i(:,:,:,0),cuv_i(:,:,:,0),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(0:rmaxc,0),cerr2_i(0:rmaxc,0))
405 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(0:rmaxc,1),cerr2_i(0:rmaxc,1))
406 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(0:rmaxc,2),cerr2_i(0:rmaxc,2))
407 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(0:rmaxc,3),cerr2_i(0:rmaxc,3))
410 write(*,*)
'CalcDred Cerr 0 =',cerr_i(0:rmaxc,0)
411 write(*,*)
'CalcDred Cerr 1 =',cerr_i(0:rmaxc,1)
412 write(*,*)
'CalcDred Cerr 2 =',cerr_i(0:rmaxc,2)
413 write(*,*)
'CalcDred Cerr 3 =',cerr_i(0:rmaxc,3)
414 if (abs(c_i(0,0,0,0)).ne.0d0) &
415 write(*,*)
'CalcDred Cacc 0 =',cerr_i(0:rmaxc,0)/abs(c_i(0,0,0,0))
416 if (abs(c_i(0,0,0,1)).ne.0d0) &
417 write(*,*)
'CalcDred Cacc 1 =',cerr_i(0:rmaxc,1)/abs(c_i(0,0,0,1))
418 if (abs(c_i(0,0,0,2)).ne.0d0) &
419 write(*,*)
'CalcDred Cacc 2 =',cerr_i(0:rmaxc,2)/abs(c_i(0,0,0,2))
420 if (abs(c_i(0,0,0,3)).ne.0d0) &
421 write(*,*)
'CalcDred Cacc 3 =',cerr_i(0:rmaxc,3)/abs(c_i(0,0,0,3))
435 cerr_i(r,i)=cerr_i(r-1,i)*errfac(i)
440 err_c(r)=maxval(cerr_i(r,0:3))
446 write(*,*)
'CalcDred err_C0:',err_c0
447 write(*,*)
'CalcDred Cerr 0 =',cerr_i(0:rmax_c,0)
448 write(*,*)
'CalcDred Cerr 1 =',cerr_i(0:rmax_c,1)
449 write(*,*)
'CalcDred Cerr 2 =',cerr_i(0:rmax_c,2)
450 write(*,*)
'CalcDred Cerr 3 =',cerr_i(0:rmax_c,3)
451 write(*,*)
'CalcDred Cerr =', err_c(0:rmax_c)
464 mm02 = elimminf2_coli(m02)
465 mm12 = elimminf2_coli(m12)
466 mm22 = elimminf2_coli(m22)
467 mm32 = elimminf2_coli(m32)
468 q10 = elimminf2_coli(p10)
469 q21 = elimminf2_coli(p21)
470 q32 = elimminf2_coli(p32)
471 q30 = elimminf2_coli(p30)
472 q31 = elimminf2_coli(p31)
473 q20 = elimminf2_coli(p20)
496 maxz = maxval(abs(
z))
499 if (
detz.ne.0d0)
then
513 zadj(1,1) = (
z(3,3)*
z(2,2)-
z(2,3)*
z(2,3))
514 zadj(2,1) = (
z(1,3)*
z(2,3)-
z(3,3)*
z(1,2))
515 zadj(3,1) = (
z(1,2)*
z(2,3)-
z(2,2)*
z(1,3))
517 zadj(2,2) = (
z(1,1)*
z(3,3)-
z(1,3)*
z(1,3))
518 zadj(3,2) = (
z(1,2)*
z(1,3)-
z(1,1)*
z(2,3))
521 zadj(3,3) = (
z(1,1)*
z(2,2)-
z(1,2)*
z(1,2))
558 fmax = max(abs(
f(1)),abs(
f(2)),abs(
f(3)))
567 mx(1:3,1:3) =
z(1:3,1:3)
571 if (
detx.ne.0d0.and.
maxz.ne.0d0)
then
586 -
f(2)*
f(3)*
z(1,3) +
f(3)*
f(3)*
z(1,2)
588 +
f(2)*
f(2)*
z(3,1) -
f(3)*
f(2)*
z(2,1)
593 -
f(3)*
f(1)*
z(2,1) +
f(3)*
f(2)*
z(1,1)
708 dscale2 = max(abs(p10*p32),abs(p21*p30),abs(p20*p31),abs(m02*m02), \
709 abs(m12*m12),abs(m22*m22),abs(m32*m32))
711 d0est = max(abs(d0_coli(p10,p21,p32,p30,p20,p31,m02,m12,m22,m32)), \
716 if(dscale2.ne.0d0)
then
732 err_inf = acc_inf*d0est
736 if (
adetx.ne.0d0)
write(*,*)
'D0est',1d0/sqrt(
adetx)
737 if (
m2max.ne.0d0)
write(*,*)
'D0est',1d0/
m2max**2
738 if (
maxz.ne.0d0)
write(*,*)
'D0est',1d0/
maxz**2
739 write(*,*)
'D0est',d0est
749 dcount(0) = dcount(0)+1
752 if (
adetx.ne.0d0)
then
753 err_d0 = acc_def_d0*max( d0est, 1d0/sqrt(
adetx) )
755 err_d0 = acc_def_d0*d0est
760 err_req_d = acc_req_d * d0est
792 if (mod(rmax,2).eq.1)
then
793 err_pv(rmax) = max( w_pv**((rmax-1)/2) * v_pv * err_d0, &
794 w_pv**((rmax-1)/2) * z_pv * err_c0, z_pv * err_c(rmax-1) )
796 write(*,*)
'CalcDred w_pv: ',w_pv,v_pv,z_pv,err_d0,err_c0,rmax
798 write(*,*)
'CalcDred err_pv con: ',err_pv(rmax), w_pv**((rmax-1)/2) * v_pv * err_d0, &
799 w_pv**((rmax-1)/2) * z_pv * err_c0, z_pv * err_c(rmax-1)
802 err_pv(rmax) = max( w_pv**(rmax/2) * err_d0, &
803 w_pv**(rmax/2-1) * v_pv * z_pv * err_c0, z_pv * err_c(rmax-1) )
805 write(*,*)
'CalcDred w_pv: ',w_pv,v_pv,z_pv
806 write(*,*)
'CalcDred err_pv con: ',err_pv(rmax), w_pv**((rmax)/2) * err_d0, &
807 w_pv**(rmax/2-1) * v_pv * z_pv * err_c0, z_pv * err_c(rmax-1)
808 write(*,*)
'CalcDred err_pv con: ',err_pv(rmax), w_pv**((rmax)/2),err_d0, &
809 w_pv**(rmax/2-1) * v_pv,err_c0, z_pv,err_c(rmax-1)
847 write(*,*)
'CalcDred w_pv2: ',w_pv2,v_pv2,z_pv2
848 write(*,*)
'CalcDred w_pv2: ',mod(rmax,2).eq.1
851 if (mod(rmax,2).eq.1)
then
859 err_pv2(rmax) = max( err_d0 * max(hw_pv2**rmax, &
860 hw_pv2*v_pv2**((rmax-1)/2) ), &
861 err_c0 * z_pv2*max(w_pv2*hw_pv2**(rmax), &
862 max(1d0,w_pv2)*hw_pv2*v_pv2**((rmax-1)/2), &
863 v_pv2**((rmax+1)/2)), &
864 err_c(rmax-1) * z_pv2 * max(hw_pv2,hw_pv2*w_pv2,v_pv2) )
867 write(*,*)
'CalcDred err_pv2: ', &
868 err_pv2(rmax) , err_d0,err_d0*w_pv2**rmax,err_d0*v_pv2**((rmax-1)/2), &
869 err_d0*w_pv2*v_pv2**((rmax-1)/2)
870 write(*,*)
'CalcDred err_pv2: ' &
871 ,err_pv2(rmax) , err_d0 * max(1d0,w_pv2**rmax,v_pv2**((rmax-1)/2), &
872 w_pv2*v_pv2**((rmax-1)/2) ) &
873 , err_c0 * z_pv2*max(w_pv2**(rmax+1), &
874 w_pv2*v_pv2**((rmax-1)/2), &
875 v_pv2**((rmax+1)/2)) &
876 , err_c(rmax-1) * max(z_pv2*w_pv2, &
877 z_pv2*w_pv2**2,z_pv2*v_pv2)
878 write(*,*)
'CalcDred err_pv2: ', &
879 err_c0 * z_pv2*w_pv2**(rmax+1), &
880 err_c0 * z_pv2* w_pv2*v_pv2**((rmax-1)/2), &
881 err_c0 * z_pv2* v_pv2**((rmax+1)/2) , &
891 err_pv2(rmax) = max( err_d0 * max(hw_pv2**rmax,v_pv2**(rmax/2)), &
892 err_c0 * z_pv2 * max(w_pv2*hw_pv2**(rmax), &
893 v_pv2**(rmax/2),w_pv2*v_pv2**(rmax/2)), &
894 err_c(rmax-1) * z_pv2 * max(hw_pv2,hw_pv2*w_pv2, v_pv2) )
913 err_pv(rmax) = err_pv(rmax)/impest_d
914 err_pv2(rmax) = err_pv2(rmax)/impest_d
920 write(*,*)
'CalcDred: err_pv',err_pv(rmax),err_pv2(rmax),err_req_d
921 write(*,*)
'CalcDred: acc_pv',err_pv(rmax)/d0est,err_pv2(rmax)/d0est,acc_req_d
937 if(use_pv.or.use_pv2)
then
939 if (min(err_pv(rmax),err_pv2(rmax)).le.err_req_d)
then
941 if (err_pv(rmax).le.err_pv2(rmax))
then
944 write(*,*)
'CalcDred: call Dpv 1 ',rmax,id,err_pv(rmax)
948 call calcdpv1(d,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,derr1,derr2)
955 dcount(1) = dcount(1)+1
956 drcalc(0:rmax) = drcalc(0:rmax)+1
962 checkest=derr(rmax)/err_pv(rmax)
963 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
964 write(*,*)
'CalcDred: estimate err_pv imprecise',err_pv(rmax),derr(rmax)
973 write(*,*)
'CalcDred: call Dpv2 1',rmax,id,err_pv2(rmax)
977 call calcdpv2(d,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,derr1,derr2)
983 dcount(2) = dcount(2)+1
984 drcalc(0:rmax)=drcalc(0:rmax)+2
989 checkest=derr(rmax)/err_pv2(rmax)
990 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
991 write(*,*)
'CalcDred: estimate err_pv2 imprecise',err_pv2(rmax),derr(rmax)
1001 err_d0 = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
1002 err_req_d = acc_req_d * abs(d(0,0,0,0))
1007 dtyp = max(abs(d(0,0,0,0)), &
1008 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
1010 dtyp = abs(d(0,0,0,0))
1012 if(dtyp.eq.0d0) dtyp = d0est
1013 err_req_d = acc_req_d * dtyp
1017 write(*,*)
'CalcDred Derr1 after PV = ',derr1
1019 write(*,*)
'CalcDred Dacc1 after PV = ',derr1/dtyp
1020 write(*,*)
'CalcDred err1_D',derr1(rmax)
1021 write(*,*)
'CalcDred Derr2 after PV = ',derr2
1023 write(*,*)
'CalcDred Dacc2 after PV = ',derr2/dtyp
1024 write(*,*)
'CalcDred err2_D',derr2(rmax)
1028 if (derr1(rmax).lt.err_req_d)
then
1029 dcount(dcalc+dcountoffset0) = dcount(dcalc+dcountoffset0)+1
1063 if (
fac_g.ge.1)
then
1066 err_g_c = err_c(rmax)
1073 err_g_cr = max( err_c(rmax), err_c0 * u_g**rmax ) * z_g
1075 err_g_exp = u_g**(rmax-1) * dtyp
1082 err_g_c = err_c(rmax)
1089 write(*,*)
'CalcDred: after Gram pars',use_g,
fac_g,
x_g,u_g,z_g,err_g_cr,err_c(rmax),err_c0,err_g_exp
1091 write(*,*)
'CalcDred: after Gram pars',err_c(rmax), err_c0 * u_g**rmax
1109 err_gm_exp = err_inf
1110 err_gm_c = err_c(rmax)
1115 err_gm_cr = max( err_c(rmax), err_c0 * u_gm**rmax ) * z_gm
1116 err_gm_c = err_gm_cr
1117 err_gm_exp = u_gm**(rmax-1) * dtyp
1123 err_gm_exp = err_inf
1124 err_gm_c = err_c(rmax)
1128 write(*,*)
'CalcDred: after mod Gram pars',use_gm,
fac_gm,
x_gm,u_gm,z_gm,err_gm_cr,err_c(rmax),err_c0,err_gm_exp
1130 write(*,*)
'CalcDred: after mod Gram pars',err_c(rmax), err_c0 * u_gm**rmax
1138 err_gm_exp = err_inf
1139 err_gm_c = err_c(rmax)
1148 v1_gy = max(1d0,
v_gy)
1154 err_gy_exp = err_inf
1155 err_gy_c = err_c(rmax+1)
1162 err_gy_cr = max( err_c(rmax) * v1_gy, err_c(rmax+1) )
1163 err_gy_c = err_gy_cr * b_gy
1164 err_gy_exp = 1d0 * dtyp
1170 err_gy_exp = err_inf
1171 err_gy_c = err_c(rmax+1)
1179 write(*,*)
'CalcDred: after GramCay pars',use_gy,
fac_gy,
x_gy,
y_gy,
v_gy,b_gy,err_gy_cr,err_gy_exp
1190 v1_gp = max(1d0,v_gp)
1196 err_gp_exp = err_inf
1197 err_gp_c = err_c(rmax)
1204 err_gp_cr = max(err_c0 * v_gp**rmax , err_c(rmax)) * z_gp
1205 err_gp_c = err_gp_cr
1206 err_gp_exp = v1_gp**(rmax-1) * dtyp
1212 err_gp_exp = err_inf
1213 err_gp_c = err_c(rmax)
1221 write(*,*)
'CalcDred: after Mom pars',use_gp,
fac_gp,
w_gp,v_gp,z_gp,err_gp_cr,err_gp_exp
1229 y1_gr = max(1d0,y_gr)
1234 if (
fac_gr.ge.1.or.2*rmax.gt.rmax_c)
then
1236 err_gr_exp = err_inf
1237 err_gr_c = err_c(rmax)
1241 err_gr_cr = err_c(rmax)
1242 err_gr_c = err_gr_cr * a_gr
1243 err_gr_exp = y1_gr * dtyp
1249 err_gr_exp = err_inf
1250 err_gr_c = err_c(rmax)
1258 write(*,*)
'CalcDred: after revGram pars',use_gr,
fac_gr,x_gr,y_gr,y1_gr,a_gr,err_gr_cr,err_c(rmax),err_c0,err_gr_exp
1259 write(*,*)
'CalcDred: after revGram pars',err_gr_c,dtyp
1265 if (abs(m02).gt.
m2scale*dprec_cll)
then
1269 v1_gpf = max(1d0,
v_gpf)
1275 err_gpf_exp = err_inf
1276 err_gpf_c = err_c(rmax+1)
1281 b_gpf = 1d0/abs(m02)
1282 err_gpf_cr = max( err_c(rmax), err_c(rmax+1) )
1283 err_gpf_c = err_gpf_cr * b_gpf
1284 err_gpf_exp = 1d0 * dtyp
1290 err_gpf_exp = err_inf
1291 err_gpf_c = err_c(rmax+1)
1298 write(*,*)
'CalcDred: after pf pars',use_gpf,
fac_gpf,
x_gpf,
y_gpf,
v_gpf,b_gpf,err_gpf_cr,err_gpf_exp,err_gpf
1305 if(use_pv.or.use_pv2.or.use_g.or.use_gy.or.use_gp.or.use_gr.or.use_gm.or.use_gpf.eqv..false.)
then
1306 call seterrflag_coli(-6)
1307 call errout_coli(
'CalcDred',
' no reduction method works', &
1310 if (errorwriteflag)
then
1311 write(nerrout_coli,fmt10)
' CalcDred: p10 = ',p10
1312 write(nerrout_coli,fmt10)
' CalcDred: p21 = ',p21
1313 write(nerrout_coli,fmt10)
' CalcDred: p32 = ',p32
1314 write(nerrout_coli,fmt10)
' CalcDred: p30 = ',p30
1315 write(nerrout_coli,fmt10)
' CalcDred: p20 = ',p20
1316 write(nerrout_coli,fmt10)
' CalcDred: p31 = ',p31
1317 write(nerrout_coli,fmt10)
' CalcDred: m02 = ',m02
1318 write(nerrout_coli,fmt10)
' CalcDred: m12 = ',m12
1319 write(nerrout_coli,fmt10)
' CalcDred: m22 = ',m22
1320 write(nerrout_coli,fmt10)
' CalcDred: m32 = ',m32
1328 write(*,*)
'CalcDred: exit'
1351 if (err_g_exp.gt.err_g_c)
then
1353 err_g_exp = err_g_exp*
fac_g
1354 err_g_c = max(err_g_cr,err_c(rmax+g)*z_g*
x_g**g)
1355 err_g(rmax) = max(err_g_exp,err_g_c)
1356 if(err_g(rmax).lt.err_req_d)
then
1359 g = min(max(g+2,3*g/2),rmax_d-rmax)
1373 if (err_gm_exp.gt.err_gm_c)
then
1375 err_gm_exp = err_gm_exp*
fac_gm
1376 err_gm_c = max(err_gm_cr,err_c(rmax+gm)*z_gm*
x_gm**gm)
1377 err_gm(rmax) = max(err_gm_exp,err_gm_c)
1378 if(err_gm(rmax).lt.err_req_d)
then
1381 gm = min(max(gm+2,3*gm/2),rmax_d-rmax)
1395 if (mod(i,2).eq.1)
then
1398 write(*,*)
'CalcDred: it gy',use_gy,err_gy_exp,err_gy_c,err_gy(rmax),err_req_d
1402 if (err_gy_exp.gt.err_gy_c.and.err_gy(rmax).gt.err_req_d)
then
1404 err_gy_exp = err_gy_exp*
fac_gy
1405 err_gy_c = b_gy*max(err_gy_cr, &
1406 max(err_c(rmax+2*gy)*v1_gy,err_c(rmax+2*gy+1))*
y_gy**gy, &
1407 max(err_c(rmax+gy)*v1_gy,err_c(rmax+gy+1))*(max(
x_gy,
v_gy*
y_gy))**gy)
1408 err_gy(rmax) = max(err_gy_exp,err_gy_c)
1411 write(*,*)
'CalcDred i gy',i,gy,err_gy_exp,err_gy_c,err_gy(rmax)
1412 write(*,*)
'CalcDred i ',err_gy_cr, &
1413 max(err_c(rmax+2*gy)*v1_gy,err_c(rmax+2*gy+1))*
y_gy**gy, &
1414 max(err_c(rmax+gy)*v1_gy,err_c(rmax+gy+1))*(max(
x_gy,
v_gy*
y_gy))**gy
1415 write(*,*)
'CalcDred i ', b_gy*err_c(rmax+2*gy)*v1_gy*
y_gy**gy, &
1416 b_gy*err_c(rmax+2*gy+1)*
y_gy**gy
1417 write(*,*)
'CalcDred i ', &
1418 b_gy,err_c(rmax+2*gy+1),
y_gy**gy
1419 write(*,*)
'CalcDred i ', b_gy*
x_gy**gy*err_c(rmax+gy)*v1_gy, &
1420 b_gy*err_c(rmax+gy+1)*
x_gy**gy
1421 write(*,*)
'CalcDred i ', &
1422 b_gy,err_c(rmax+gy+1),
x_gy**gy,
x_gy,gy
1425 if(err_gy(rmax).lt.err_req_d)
then
1428 gy = min(max(gy+2,2*gy),(rmax_d-rmax)/2)
1443 if (err_gp_exp.gt.err_gp_c.and.err_gp(rmax).gt.err_req_d)
then
1445 err_gp_exp = err_gp_exp*
fac_gp
1446 err_gp_c = max(err_c(rmax+gp)*z_gp*
w_gp**gp,err_gp_cr)
1447 err_gp(rmax) = max(err_gp_exp,err_gp_c)
1448 if(err_gp(rmax).lt.err_req_d)
then
1451 gp = min(max(gp+2,3*gp/2),rmax_d-rmax)
1459 if (mod(i,2).eq.1.and.i.le.rmax_c-2*rmax)
then
1467 if (err_gr_exp.gt.err_gr_c.and.err_gr(rmax).gt.err_req_d)
then
1469 err_gr_exp = err_gr_exp*
fac_gr
1470 err_gr_c = a_gr*max(err_gr_cr, &
1471 max(err_c(rmax+gr),err_c(rmax+gr+1)*y_gr)*
fac_gr**gr)
1472 err_gr(rmax) = max(err_gr_exp,err_gr_c)
1476 if(err_gr(rmax).lt.err_req_d)
then
1481 gr = min(max(gr+2,2*gr),rmax_d-rmax,max(0,(rmax_c-2*rmax)/2))
1491 if (mod(i,2).eq.1)
then
1494 write(*,*)
'CalcDred: it gy',use_gy,err_gy_exp,err_gy_c,err_gy(rmax),err_req_d
1498 if (err_gpf_exp.gt.err_gpf_c.and.err_gpf(rmax).gt.err_req_d)
then
1500 err_gpf_exp = err_gpf_exp*
fac_gpf
1501 err_gpf_c = b_gpf*max(err_gpf_cr, &
1502 max(err_c(rmax+2*gpf)*v1_gpf,err_c(rmax+2*gpf+1))*
y_gpf**gpf, &
1503 max(err_c(rmax+gpf)*v1_gpf,err_c(rmax+gpf+1))*(max(
x_gpf,
v_gpf*
y_gpf))**gpf)
1504 err_gpf(rmax) = max(err_gpf_exp,err_gpf_c)
1507 write(*,*)
'CalcDred i gpf',i,gpf,err_gpf_exp,err_gpf_c,err_gpf(rmax)
1508 write(*,*)
'CalcDred i ',err_gpf_cr, &
1509 max(err_c(rmax+2*gpf)*v1_gpf,err_c(rmax+2*gpf+1))*
y_gpf**gpf, &
1510 max(err_c(rmax+gpf)*v1_gpf,err_c(rmax+gpf+1))*(max(
x_gpf,
v_gpf*
y_gpf))**gpf
1511 write(*,*)
'CalcDred i ', b_gpf*err_c(rmax+2*gpf)*v1_gpf*
y_gpf**gpf, &
1512 b_gpf*err_c(rmax+2*gpf+1)*
y_gpf**gpf
1513 write(*,*)
'CalcDred i ', &
1514 b_gpf,err_c(rmax+2*gpf+1),
y_gpf**gpf
1515 write(*,*)
'CalcDred i ', b_gpf*
x_gpf**gpf*err_c(rmax+gpf)*v1_gpf, &
1516 b_gpf*err_c(rmax+gpf+1)*
x_gpf**gpf
1517 write(*,*)
'CalcDred i ', &
1521 if(err_gpf(rmax).lt.err_req_d)
then
1524 gpf = min(max(gpf+2,2*gpf),(rmax_d-rmax)/2)
1541 err_g(rmax) = err_g(rmax)/impest_d
1542 err_gy(rmax) = err_gy(rmax)/impest_d
1543 err_gp(rmax) = err_gp(rmax)/impest_d
1544 err_gr(rmax) = err_gr(rmax)/impest_d
1545 err_gm(rmax) = err_gm(rmax)/impest_d
1546 err_gpf(rmax) = err_gpf(rmax)/impest_d
1549 write(*,*)
'iexp=',iexp
1551 write(*,*)
'errexp=',err_g_exp,err_gy_exp,err_gp_exp,err_gr_exp,err_gm_exp,err_gpf_exp,err_req_d
1552 write(*,*)
'errexptot=',i,g,err_g(rmax),gy,err_gy(rmax),gp,err_gp(rmax), &
1554 gr,err_gr(rmax),gm,err_gm(rmax),gpf,err_gpf(rmax)
1555 write(*,*)
'accexptot=',i,g,err_g(rmax)/dtyp,gy,err_gy(rmax)/dtyp, &
1556 gp,err_gp(rmax)/dtyp,gr,err_gr(rmax)/dtyp,gm,err_gm(rmax)/dtyp, &
1557 gpf,err_gpf(rmax)/dtyp
1577 call calcdg(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,g,g,id,derr1_alt,derr2_alt)
1579 derr_alt = derr2_alt
1581 derr_alt = derr1_alt
1583 dcount(3) = dcount(3)+1
1584 drcalc(0:rmax)=drcalc(0:rmax)+4
1586 drmethod_alt(0:rmax)=4
1589 checkest=derr_alt(rmax)/err_g(rmax)
1590 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1591 write(*,*)
'CalcDred: estimate err_g imprecise ',err_g(rmax),derr_alt(rmax)
1597 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1601 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1603 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1604 write(*,*)
'CalcDred after exp method=',drmethod
1608 call calcdgy(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,gy,gy,id,derr1_alt,derr2_alt)
1610 derr_alt = derr2_alt
1612 derr_alt = derr1_alt
1614 dcount(4) = dcount(4)+1
1615 drcalc(0:rmax)=drcalc(0:rmax)+8
1617 drmethod_alt(0:rmax)=8
1620 checkest=derr_alt(rmax)/err_gy(rmax)
1621 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1622 write(*,*)
'CalcDred: estimate err_gy imprecise',err_gy(rmax),derr_alt(rmax),checkest
1627 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1630 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1631 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1632 write(*,*)
'CalcDred after exp method=',drmethod
1640 call calcdgp(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,gp,gp,id,derr1_alt,derr2_alt)
1642 derr_alt = derr2_alt
1644 derr_alt = derr1_alt
1646 dcount(5) = dcount(5)+1
1647 drcalc(0:rmax)=drcalc(0:rmax)+16
1649 drmethod_alt(0:rmax)=16
1652 checkest=derr_alt(rmax)/err_gp(rmax)
1653 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1654 write(*,*)
'CalcDred: estimate err_gp imprecise',err_gp(rmax),derr_alt(rmax)
1659 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1662 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1663 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1664 write(*,*)
'CalcDred after exp method=',drmethod
1668 call calcdgr(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,gr,gr,id,derr1_alt,derr2_alt)
1670 derr_alt = derr2_alt
1672 derr_alt = derr1_alt
1674 dcount(6) = dcount(6)+1
1675 drcalc(0:rmax)=drcalc(0:rmax)+32
1677 drmethod_alt(0:rmax)=32
1680 checkest=derr_alt(rmax)/err_gr(rmax)
1681 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1682 write(*,*)
'CalcDred: estimate err_gr imprecise',err_gr(rmax),derr_alt(rmax)
1687 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1690 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1691 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1692 write(*,*)
'CalcDred after exp method=',drmethod
1697 call calcdgm(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,gm,gm,id,derr1_alt,derr2_alt)
1699 derr_alt = derr2_alt
1701 derr_alt = derr1_alt
1703 dcount(7) = dcount(7)+1
1704 drcalc(0:rmax)=drcalc(0:rmax)+64
1706 drmethod_alt(0:rmax)=64
1709 checkest=derr_alt(rmax)/err_gm(rmax)
1710 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1711 write(*,*)
'CalcDred: estimate err_gm imprecise',err_gm(rmax),derr_alt(rmax)
1716 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1719 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1720 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1721 write(*,*)
'CalcDred after exp method=',drmethod
1726 call calcdgpf(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,gpf,gpf,id,derr1_alt,derr2_alt)
1728 derr_alt = derr2_alt
1730 derr_alt = derr1_alt
1732 dcount(7) = dcount(7)+1
1733 drcalc(0:rmax)=drcalc(0:rmax)+64
1735 drmethod_alt(0:rmax)=64
1738 checkest=derr_alt(rmax)/err_gpf(rmax)
1739 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1740 write(*,*)
'CalcDred: estimate err_gpf imprecise',err_gpf(rmax),derr_alt(rmax),checkest
1745 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1748 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1749 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1750 write(*,*)
'CalcDred after exp method=',drmethod
1760 if(.not.lerr_d0.and.iexp.ne.0)
then
1762 err_d0 = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
1771 dtyp = max(abs(d(0,0,0,0)), &
1772 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
1774 dtyp = abs(d(0,0,0,0))
1776 err_req_d = acc_req_d * dtyp
1779 write(*,*)
'CalcDred ',rmax,derr1(rmax),err_req_d
1782 if (derr1(rmax).le.err_req_d)
then
1783 dcount(dcalc+dcountoffset0) = dcount(dcalc+dcountoffset0)+1
1793 write(*,*)
'CalcDred no optimal method'
1794 write(*,*)
'err_req_D=',err_req_d
1795 write(*,*)
'err_est=',err_pv(rmax),err_pv2(rmax),err_g(rmax),err_gy(rmax),err_gp(rmax),err_gr(rmax),err_gpf(rmax)
1802 if(use_pv.and.mod(drcalc(r),2).ne.1)
then
1806 if (mod(r,2).eq.1)
then
1807 err_pv(r) = max( w_pv**((r-1)/2) * v_pv * err_d0, &
1808 w_pv**((r-1)/2) * z_pv * err_c0, z_pv * err_c(r-1) )
1815 else if (r.ne.0)
then
1816 err_pv(r) = max( w_pv**(r/2) * err_d0, &
1817 w_pv**(r/2-1) * v_pv * z_pv * err_c0, z_pv * err_c(r-1) )
1827 err_pv(r) = err_pv(r)/impest_d
1830 if (use_pv2.and.mod(drcalc(r),4)-mod(drcalc(r),2).ne.2)
then
1836 if (mod(r,2).eq.1)
then
1844 err_pv2(r) = max( err_d0 * max(hw_pv2**r, &
1845 hw_pv2*v_pv2**((r-1)/2) ), &
1846 err_c0 * z_pv2* max(w_pv2*w_pv2**(r), &
1847 hw_pv2*v_pv2**((r-1)/2), &
1848 w_pv2*hw_pv2*v_pv2**((r-1)/2), &
1849 v_pv2**((r+1)/2)), &
1850 err_c(r-1) * z_pv2 * max(hw_pv2,w_pv2*hw_pv2**2,v_pv2) )
1863 else if (r.ne.0)
then
1869 err_pv2(r) = max( err_d0 * max(hw_pv2**r,v_pv2**(r/2)), &
1870 err_c0 * z_pv2 * max(w_pv2*hw_pv2**(r), &
1871 v_pv2**(r/2),w_pv2*v_pv2**(r/2)), &
1872 err_c(r-1) * z_pv2 * max(hw_pv2, w_pv2*hw_pv2**2, v_pv2) )
1887 err_pv2(r) = err_pv2(r)/impest_d
1890 if (mod(drcalc(r),8)-mod(drcalc(r),4).ne.4.and.use_g)
then
1892 err_g_cr = max( err_c(r), err_c0 * u_g**r ) * z_g
1894 err_g_exp = u_g**(r-1) * dtyp
1899 err_g_exp = err_g_exp*
fac_g
1900 err_g_c = max(err_g_cr,err_c(r+g)*z_g*
x_g**g)
1901 err_g(r) = max(err_g_exp,err_g_c)
1902 if (err_g_exp.lt.err_g_c.or.err_g(r).lt.err_req_d)
exit
1905 g = min(max(g+2,2*g),rmax_d-r)
1907 err_g(r) = err_g(r)/impest_d
1910 if (mod(drcalc(r),16)-mod(drcalc(r),8).ne.8.and.use_gy)
then
1912 err_gy_cr = max( err_c(r) * v1_gy, err_c(r+1) )
1913 err_gy_c = err_gy_cr * b_gy
1914 err_gy_exp = 1d0 * dtyp
1919 if (mod(i,2).eq.1)
then
1921 err_gy_exp = err_gy_exp*
fac_gy
1922 err_gy_c = b_gy*max(err_gy_cr, &
1923 max(err_c(r+2*gy)*v1_gy,err_c(r+2*gy+1))*
y_gy**gy, &
1924 max(err_c(r+gy)*v1_gy,err_c(r+gy+1))*(max(
x_gy,
v_gy*
y_gy))**gy)
1925 err_gy(r) = max(err_gy_exp,err_gy_c)
1926 if (err_gy_exp.lt.err_gy_c.or.err_gy(r).lt.err_req_d)
exit
1930 gy = min(max(gy+2,2*gy),(rmax_d-r)/2)
1932 err_gy(r) = err_gy(r)/impest_d
1935 if (mod(drcalc(r),32)-mod(drcalc(r),16).ne.16.and.use_gp)
then
1937 err_gp_cr = max(err_c0*v_gp**r,err_c(r))*z_gp
1938 err_gp_exp = v1_gp**(r-1) * dtyp
1943 err_gp_exp = err_gp_exp*
fac_gp
1944 err_gp_c = max(err_c(r+gp)*z_gp*
w_gp**gp,err_gp_cr)
1945 err_gp(r) = max(err_gp_exp,err_gp_c)
1946 if (err_gp_exp.lt.err_gp_c.or.err_gp(r).lt.err_req_d)
exit
1949 gp = min(max(gp+2,3*gp/2),rmax_d-r)
1951 err_gp(r) = err_gp(r)/impest_d
1954 if (mod(drcalc(r),64)-mod(drcalc(r),32).ne.32.and.use_gr)
then
1956 err_gr_cr = err_c(r)
1957 err_gr_c = err_gr_cr * a_gr
1958 err_gr_exp = y1_gr * dtyp
1962 do i=0,min(rmax_d-r,rmax_c-2*r)
1963 if (mod(i,2).eq.1)
then
1965 err_gr_exp = err_gr_exp*
fac_gr
1966 err_gr_c = a_gr*max(err_gr_cr, &
1967 max(err_c(r+gr),err_c(r+gr+1)*y_gr)*
fac_gr**gr)
1968 err_gr(r) = max(err_gr_exp,err_gr_c)
1971 write(*,*)
'CalcDgr err_gr',i,gr,err_gr_exp,err_gr_c,err_gr(r),err_req_d
1974 if (err_gr_exp.lt.err_gr_c.or.err_gr(r).lt.err_req_d)
exit
1980 gr = min(max(gr+2,2*gr),rmax_d-r,max(0,(rmax_c-2*r)/2))
1982 err_gr(r) = err_gr(r)/impest_d
1986 if (mod(drcalc(r),128)-mod(drcalc(r),64).ne.64.and.use_gpf)
then
1988 err_gpf_cr = max( err_c(r) * v1_gpf, err_c(r+1) )
1989 err_gpf_c = err_gpf_cr * b_gpf
1990 err_gpf_exp = 1d0 * dtyp
1995 if (mod(i,2).eq.1)
then
1997 err_gpf_exp = err_gpf_exp*
fac_gpf
1998 err_gpf_c = b_gpf*max(err_gpf_cr, &
1999 max(err_c(r+2*gpf)*v1_gpf,err_c(r+2*gpf+1))*
y_gpf**gpf, &
2000 max(err_c(r+gpf)*v1_gpf,err_c(r+gpf+1))*(max(
x_gpf,
v_gpf*
y_gpf))**gpf)
2001 err_gpf(r) = max(err_gpf_exp,err_gpf_c)
2002 if (err_gpf_exp.lt.err_gpf_c.or.err_gpf(r).lt.err_req_d)
exit
2006 gpf = min(max(gpf+2,2*gpf),(rmax_d-r)/2)
2008 err_gpf(r) = err_gpf(r)/impest_d
2014 if (mod(drcalc(r),128)-mod(drcalc(r),64).ne.64.and.use_gm)
then
2016 err_gm_cr = max( err_c(r), err_c0 * u_gm**r ) * z_gm
2017 err_gm_c = err_gm_cr
2018 err_gm_exp = u_gm**(r-1) * dtyp
2023 err_gm_exp = err_gm_exp*
fac_gm
2024 err_gm_c = max(err_gm_cr,err_c(r+gm)*z_gm*
x_gm**gm)
2025 err_gm(r) = max(err_gm_exp,err_gm_c)
2026 if (err_gm_exp.lt.err_gm_c.or.err_gm(r).lt.err_req_d)
exit
2029 gm = min(max(gm+2,2*gm),rmax_d-r)
2031 err_gm(r) = err_gm(r)/impest_d
2037 write(*,*)
'CalcDred: bef final loop ord methods',r,g,gy,gp,gr,gm,gpf
2038 write(*,*)
'CalcDred: bef final loop err methods',r,err_pv(r),err_pv2(r), &
2039 err_g(r),err_gy(r),err_gp(r),err_gr(r),err_gm(r),err_gpf(r)
2040 write(*,*)
'CalcDred: bef final loop acc methods',r,err_pv(r)/dtyp,err_pv2(r)/dtyp, &
2041 err_g(r)/dtyp,err_gy(r)/dtyp,err_gp(r)/dtyp, &
2042 err_gr(r)/dtyp,err_gm(r)/dtyp,err_gpf(r)/dtyp
2043 write(*,*)
'CalcDred: bef final loop method',r,drcalc(r),drmethod(r)
2047 if (min(err_pv(r),err_pv2(r)).le.min(err_g(r),err_gy(r),err_gp(r),err_gr(r),err_gpf(r)) &
2048 .and.min(err_pv(r),err_pv2(r)).lt.err_inf)
then
2050 if (use_pv.and.err_pv(r).le.err_pv2(r).and.mod(drcalc(r),2).ne.1)
then
2064 write(*,*)
'CalcDred: call Dpv 2',r,id
2073 call calcdpv1(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,id,derr1_alt,derr2_alt)
2075 call calcdpv1(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2076 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,id,derr1_alt(0:r),derr2_alt(0:r))
2079 derr_alt = derr2_alt
2081 derr_alt = derr1_alt
2083 dcount(11) = dcount(11)+1
2084 drcalc(0:r)=drcalc(0:r)+1
2087 checkest=derr_alt(r)/err_pv(r)
2090 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2091 write(*,*)
'CalcDred: estimate err_pv imprecise',err_pv(r),derr_alt(r)
2096 write(*,*)
'final loop r Dpv D(1,0,1,0)',r,d_alt(1,0,1,0),d(1,0,1,0)
2097 write(*,*)
'final loop r Dpv Derr',derr_alt(2),derr(2)
2099 err_pv(0:r)=derr_alt(0:r)
2101 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2102 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2105 dtyp = max(abs(d(0,0,0,0)), &
2106 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2108 dtyp = abs(d(0,0,0,0))
2110 err_req_d = acc_req_d * dtyp
2113 write(*,*)
'CalcDred: after pv 2nd try Dmethod_alt=',drmethod_alt
2114 write(*,*)
'CalcDred: after pv 2nd try Derr_alt(r)=',derr_alt
2115 write(*,*)
'CalcDred: after pv 2nd try Dacc_alt(r)=',derr_alt/dtyp
2116 write(*,*)
'CalcDred: after pv 2nd try Dmethod=',drmethod
2117 write(*,*)
'CalcDred: after pv 2nd try Derr(r)=',derr
2118 write(*,*)
'CalcDred: after pv 2nd try Dacc(r)=',derr/dtyp
2126 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2128 elseif (use_pv2.and.err_pv2(r).le.err_pv(r).and.mod(drcalc(r),4)-mod(drcalc(r),2).ne.2)
then
2142 write(*,*)
'CalcDred: call Dpv2 2',r,id
2145 call calcdpv2(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,id,derr1_alt,derr2_alt)
2147 call calcdpv2(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2148 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,id,derr1_alt(0:r),derr2_alt(0:r))
2152 derr_alt = derr2_alt
2154 derr_alt = derr1_alt
2156 dcount(12) = dcount(12)+1
2157 drcalc(0:r)=drcalc(0:r)+2
2160 checkest=derr_alt(r)/err_pv2(r)
2163 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2164 write(*,*)
'CalcDred: estimate err_pv2 imprecise',err_pv2(r),derr_alt(r)
2167 err_pv2(0:r)=derr_alt(0:r)
2171 write(*,*)
'final loop r Dpv2 D(1,0,1,0)',r,d_alt(1,0,1,0),d(1,0,1,0)
2172 write(*,*)
'final loop r Dpv2 Derr',derr_alt(2),derr(2)
2175 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2176 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2179 dtyp = max(abs(d(0,0,0,0)), &
2180 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2182 dtyp = abs(d(0,0,0,0))
2184 err_req_d = acc_req_d * dtyp
2186 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2189 write(*,*)
'CalcDred: after pv 2nd try Dmethod=',drmethod
2190 write(*,*)
'CalcDred: after pv 2nd try Derr(r)=',derr
2191 write(*,*)
'CalcDred: after pv 2nd try Dacc(r)=',derr/dtyp
2199 write(*,*)
'CalcDred: explore exps once more'
2202 if (use_g.and.err_g(r).le.min(err_gy(r),err_gp(r),err_gr(r),err_gpf(r)) &
2203 .and.mod(drcalc(r),8)-mod(drcalc(r),4).ne.4)
then
2217 call calcdg(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,g,rmax_d,id,derr1_alt,derr2_alt)
2219 call calcdg(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2220 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,g,rmax_d,id,derr1_alt(0:r),derr2_alt(0:r))
2223 derr_alt = derr2_alt
2225 derr_alt = derr1_alt
2227 dcount(13) = dcount(13)+1
2228 drcalc(0:r)=drcalc(0:r)+4
2231 checkest=derr_alt(r)/err_g(r)
2234 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2235 write(*,*)
'CalcDred: estimate err_g imprecise ',err_g(r),derr_alt(r)
2239 err_g(0:r)=derr_alt(0:r)
2241 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2242 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2245 dtyp = max(abs(d(0,0,0,0)), &
2246 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2248 dtyp = abs(d(0,0,0,0))
2250 err_req_d = acc_req_d * dtyp
2253 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2254 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2255 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2264 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2266 else if (use_gy.and.err_gy(r).le.min(err_g(r),err_gp(r),err_gr(r),err_gpf(r)) &
2267 .and.mod(drcalc(r),16)-mod(drcalc(r),8).ne.8)
then
2281 call calcdgy(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gy,(rmax_d)/2,id,derr1_alt,derr2_alt)
2283 call calcdgy(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2284 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gy,(rmax_d)/2,id,derr1_alt(0:r),derr2_alt(0:r))
2287 derr_alt = derr2_alt
2289 derr_alt = derr1_alt
2291 dcount(14) = dcount(14)+1
2292 drcalc(0:r)=drcalc(0:r)+8
2295 checkest=derr_alt(r)/err_gy(r)
2298 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2299 write(*,*)
'CalcDred: estimate err_gy imprecise',err_gy(r),derr_alt(r)
2303 err_gy(0:r)=derr_alt(0:r)
2305 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2306 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2309 dtyp = max(abs(d(0,0,0,0)), &
2310 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2312 dtyp = abs(d(0,0,0,0))
2314 err_req_d = acc_req_d * dtyp
2317 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2318 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2319 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2328 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2330 elseif (use_gp.and.err_gp(r).le.min(err_g(r),err_gy(r),err_gr(r),err_gpf(r)) &
2331 .and.mod(drcalc(r),32)-mod(drcalc(r),16).ne.16)
then
2345 call calcdgp(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gp,rmax_d,id,derr1_alt,derr2_alt)
2347 call calcdgp(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2348 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gp,rmax_d,id,derr1_alt(0:r),derr2_alt(0:r))
2351 derr_alt = derr2_alt
2353 derr_alt = derr1_alt
2355 dcount(15) = dcount(15)+1
2356 drcalc(0:r)=drcalc(0:r)+16
2358 drmethod_alt(0:r)=16
2359 checkest=derr_alt(r)/err_gp(r)
2362 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2363 write(*,*)
'CalcDred: estimate err_gp imprecise',err_gp(r),derr_alt(r)
2367 err_gp(0:r)=derr_alt(0:r)
2369 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2370 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2373 dtyp = max(abs(d(0,0,0,0)), &
2374 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2376 dtyp = abs(d(0,0,0,0))
2378 err_req_d = acc_req_d * dtyp
2381 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2382 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2383 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2386 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2388 elseif (use_gr.and.err_gr(r).le.min(err_g(r),err_gy(r),err_gp(r),err_gpf(r)) &
2389 .and.mod(drcalc(r),64)-mod(drcalc(r),32).ne.32)
then
2403 call calcdgr(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gr,rmax_d,id,derr1_alt,derr2_alt)
2405 call calcdgr(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2406 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gr,rmax_d,id,derr1_alt(0:r),derr2_alt(0:r))
2409 derr_alt = derr2_alt
2411 derr_alt = derr1_alt
2413 dcount(16) = dcount(16)+1
2414 drcalc(0:r)=drcalc(0:r)+32
2416 drmethod_alt(0:r)=32
2417 checkest=derr_alt(r)/err_gr(r)
2420 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2421 write(*,*)
'CalcDred: estimate err_gr imprecise',err_gr(r),derr_alt(r)
2425 err_gr(0:r)=derr_alt(0:r)
2427 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2428 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2431 dtyp = max(abs(d(0,0,0,0)), &
2432 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2434 dtyp = abs(d(0,0,0,0))
2436 err_req_d = acc_req_d * dtyp
2439 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2440 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2441 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2444 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2453 else if (use_gpf.and.err_gpf(r).le.min(err_g(r),err_gy(r),err_gp(r),err_gr(r)) &
2454 .and.mod(drcalc(r),128)-mod(drcalc(r),64).ne.64)
then
2457 call calcdgpf(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gpf,(rmax_d)/2,id,derr1_alt,derr2_alt)
2459 call calcdgpf(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2460 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gpf,(rmax_d)/2,id,derr1_alt(0:r),derr2_alt(0:r))
2463 derr_alt = derr2_alt
2465 derr_alt = derr1_alt
2467 dcount(17) = dcount(17)+1
2468 drcalc(0:r)=drcalc(0:r)+64
2470 drmethod_alt(0:r)=64
2471 checkest=derr_alt(r)/err_gpf(r)
2474 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2475 write(*,*)
'CalcDred: estimate err_gpf imprecise',err_gpf(r),derr_alt(r)
2479 err_gpf(0:r)=derr_alt(0:r)
2480 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2481 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2484 dtyp = max(abs(d(0,0,0,0)), &
2485 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2487 dtyp = abs(d(0,0,0,0))
2489 err_req_d = acc_req_d * dtyp
2492 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2493 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2494 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2503 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2506 else if (use_gm.and.err_gm(r).le.min(err_gy(r),err_gp(r),err_gr(r),err_g(r)) &
2507 .and.mod(drcalc(r),128)-mod(drcalc(r),64).ne.64)
then
2521 call calcdgm(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,g,rmax_d,id,derr1_alt,derr2_alt)
2523 call calcdgm(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2524 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,g,rmax_d,id,derr1_alt(0:r),derr2_alt(0:r))
2527 derr_alt = derr2_alt
2529 derr_alt = derr1_alt
2531 dcount(17) = dcount(17)+1
2532 drcalc(0:r)=drcalc(0:r)+64
2534 drmethod_alt(0:r)=64
2535 checkest=derr_alt(r)/err_gm(r)
2538 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2539 write(*,*)
'CalcDred: estimate err_g imprecise ',err_gm(r),derr_alt(r)
2543 err_gm(0:r)=derr_alt(0:r)
2545 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2546 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2549 dtyp = max(abs(d(0,0,0,0)), &
2550 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2552 dtyp = abs(d(0,0,0,0))
2554 err_req_d = acc_req_d * dtyp
2557 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2558 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2559 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2568 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2577 if(.not.lerr_d0)
then
2579 err_d0 = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
2595 write(*,*)
'CalcDred: final loop err methods',r,err_pv(r),err_pv2(r), &
2596 err_g(r),err_gy(r),err_gp(r),err_gr(r),err_gm(r),err_gpf(r)
2597 write(*,*)
'CalcDred: final loop acc methods',r,err_pv(r)/dtyp,err_pv2(r)/dtyp, &
2598 err_g(r)/dtyp,err_gy(r)/dtyp,err_gp(r)/dtyp, &
2599 err_gr(r)/dtyp,err_gm(r)/dtyp,err_gpf(r)/dtyp
2600 write(*,*)
'CalcDred: final loop method',r,drcalc(r),drmethod(r)
2605 norm = abs(d(0,0,0,0))
2611 norm = max(norm,abs(d(0,n1,n2,n3)))
2615 acc_d = derr(rmax)/norm
2617 dcount(dcalc+dcountoffset0) = dcount(dcalc+dcountoffset0)+1
2620 write(*,*)
'CalcDred final err_D=',derr
2621 write(*,*)
'CalcDred final acc_D=',derr/norm,critacc_coli
2622 write(*,*)
'CalcDred final method_D=',drmethod
2625 if (acc_d.gt.sqrt(reqacc_coli))
then
2626 dcount(dcalc+dcountoffset3) = dcount(dcalc+dcountoffset3)+1
2629 if (acc_d.gt.reqacc_coli)
then
2630 dcount(dcalc+dcountoffset1) = dcount(dcalc+dcountoffset1)+1
2633 if (acc_d.gt.critacc_coli)
then
2635 dcount(dcalc+dcountoffset2) = dcount(dcalc+dcountoffset2)+1
2638 write(*,*)
'CritPoint D',critacc_coli,acc_d
2639 write(*,*)
'CritPoint D',critpointcntd,maxcritpointd
2646 #ifdef CritPointsCOLI
2647 critpointcntd = critpointcntd + 1
2649 if (critpointcntd.le.maxcritpointd.and.monitoring)
then
2651 call critpointsout_coli(
'D_coli',acc_d)
2652 write(ncpout_coli,*)
'arguments of CalcDred_coli:'
2653 write(ncpout_coli,*)
'rank = ', rmax
2654 write(ncpout_coli,fmt1)
'p10 = ', p10
2655 write(ncpout_coli,fmt1)
'p21 = ', p21
2656 write(ncpout_coli,fmt1)
'p32 = ', p32
2657 write(ncpout_coli,fmt1)
'p30 = ', p30
2658 write(ncpout_coli,fmt1)
'p20 = ', p20
2659 write(ncpout_coli,fmt1)
'p31 = ', p31
2660 write(ncpout_coli,fmt1)
'm02 = ', m02
2661 write(ncpout_coli,fmt1)
'm12 = ', m12
2662 write(ncpout_coli,fmt1)
'm22 = ', m22
2663 write(ncpout_coli,fmt1)
'm32 = ', m32
2664 if (critpointcntd.eq.maxcritpointd)
then
2665 write(ncpout_coli,*)
2666 write(ncpout_coli,*)
2667 write(ncpout_coli,*)
2668 write(ncpout_coli,*)
'***********************************************************'
2669 write(ncpout_coli,*)
2670 write(ncpout_coli,*)
' Further output of bad D functions will be suppressed '
2677 write(*,*)
'CalcDred exit D(1,0,0,0)',r,d_alt(1,0,0,0),d(1,0,0,0)
2693 subroutine calcduv(Duv,Cuv_0,m02,f,rmax,id)
2695 integer,
intent(in) :: rmax,id
2696 double complex,
intent(in) :: m02,f(3)
2697 double complex,
intent(inout) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
2698 double complex,
intent(in) :: Cuv_0(0:rmax-1,0:rmax-1,0:rmax-1,0:rmax-1)
2699 integer :: r,n0,n1,n2,n3
2702 duv(0:min(rmax,1),:,:,:) = 0d0
2708 do n0=max(2,r-rmax),r/2
2713 duv(n0,n1,n2,n3) = (cuv_0(n0-1,n1,n2,n3) + 2*m02*duv(n0-1,n1,n2,n3) &
2714 + f(1)*duv(n0-1,n1+1,n2,n3) &
2715 + f(2)*duv(n0-1,n1,n2+1,n3) &
2716 + f(3)*duv(n0-1,n1,n2,n3+1)) / (2*(r-1))
2737 subroutine calcdpv1(D,Duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,Derr,Derr2)
2741 integer,
intent(in) :: rmax,id
2742 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
2743 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
2744 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
2745 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
2746 double complex :: C_0(0:rmax-1,0:rmax-1,0:rmax-1,0:rmax-1), Cuv_0(0:rmax-1,0:rmax-1,0:rmax-1,0:rmax-1)
2747 double complex :: C_i(0:rmax-1,0:rmax-1,0:rmax-1,3), Cuv_i(0:rmax-1,0:rmax-1,0:rmax-1,3)
2748 double complex :: D_alt(0:rmax,0:rmax,0:rmax,0:rmax)
2749 double precision :: Cerr_i(0:rmax-1,0:3),Cerr2_i(0:rmax-1,0:3)
2754 double complex :: Smod(3)
2755 double complex :: D0_coli, elimminf2_coli
2758 double precision :: D00_err(0:rmax),Dij_err(0:rmax),Cij_err(0:rmax-1)
2759 double precision :: D00_err2(0:rmax),Dij_err2(0:rmax),Cij_err2(0:rmax-1)
2760 integer :: rmaxC,r,n0,n1,n2,n3,nn0,nn1,nn2,nn3,i,j
2761 integer :: bin,k,nid(0:3)
2766 d(0,0,0,0) = d0_coli(p10,p21,p32,p30,p20,p31,m02,m12,m22,m32)
2770 derr(0) = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
2771 derr2(0) = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
2773 if (rmax.eq.0)
return
2798 if (mod(id/bin,2).eq.0)
then
2806 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0))
2807 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1))
2808 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2))
2809 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3))
2812 write(*,*)
'CalcDpv1 Cerr_i=',cerr_i(:,0)
2813 write(*,*)
'CalcDpv1 Cerr_i=',cerr_i(:,1)
2814 write(*,*)
'CalcDpv1 Cerr_i=',cerr_i(:,2)
2815 write(*,*)
'CalcDpv1 Cerr_i=',cerr_i(:,3)
2822 n0 = (rmaxc-n1-n2-n3)
2823 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
2824 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
2825 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
2826 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
2839 dij_err(0) = derr(0)
2840 cij_err = max(cerr_i(:,0),cerr_i(:,1),cerr_i(:,2),cerr_i(:,3))
2844 dij_err2(0) = derr2(0)
2845 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,1),cerr2_i(:,2),cerr2_i(:,3))
2848 write(*,*)
'CalcDpv1 Cij_err=',cij_err
2849 write(*,*)
'CalcDpv1 Dij_err(0)=',dij_err(0)
2850 write(*,*)
'CalcDpv1 test :', &
2868 d(n0,n1,n2,n3) = 4*duv(n0,n1,n2,n3) +
detx/
detz*d(n0-1,n1,n2,n3)
2870 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) &
2876 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) &
2879 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) &
2884 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) &
2887 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) &
2892 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) &
2895 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) &
2899 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) / (2*(r-1))
2917 else if (n2.ge.1)
then
2935 smod(1) = smod(1) - 2*nn1*d(n0+1,nn1-1,nn2,nn3)
2937 smod(1) = smod(1) + c_i(n0,nn2,nn3,1)
2941 smod(2) = smod(2) - 2*nn2*d(n0+1,nn1,nn2-1,nn3)
2943 smod(2) = smod(2) + c_i(n0,nn1,nn3,2)
2947 smod(3) = smod(3) - 2*nn3*d(n0+1,nn1,nn2,nn3-1)
2949 smod(3) = smod(3) + c_i(n0,nn1,nn2,3)
2952 d(n0,n1,n2,n3) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
2953 +
zadj(3,j)*smod(3) &
2954 -
zadjs(j)*c_0(n0,nn1,nn2,nn3) &
2971 if (n1.ge.1.and.n2+n3.ge.1)
then
2991 smod(1) = smod(1) - 2*nn1*d(n0+1,nn1-1,nn2,nn3)
2993 smod(1) = smod(1) + c_i(n0,nn2,nn3,1)
2997 smod(2) = smod(2) - 2*nn2*d(n0+1,nn1,nn2-1,nn3)
2999 smod(2) = smod(2) + c_i(n0,nn1,nn3,2)
3003 smod(3) = smod(3) - 2*nn3*d(n0+1,nn1,nn2,nn3-1)
3005 smod(3) = smod(3) + c_i(n0,nn1,nn2,3)
3008 d_alt(n0,n1,n2,n3) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
3009 +
zadj(3,j)*smod(3) &
3010 -
zadjs(j)*c_0(n0,nn1,nn2,nn3) &
3013 derr(r)=max(derr(r),abs(d(n0,n1,n2,n3)-d_alt(n0,n1,n2,n3)))
3014 derr2(r)=max(derr2(r),abs(d(n0,n1,n2,n3)-d_alt(n0,n1,n2,n3)))
3029 d00_err(r) = max(2*abs(m02)*dij_err(r-2), cerr_i(r-2,0), &
3035 dij_err(r) = max(
maxzadjf*dij_err(r-1), &
3040 d00_err2(r) = max(2*abs(m02)*dij_err2(r-2), cerr2_i(r-2,0), &
3042 maxzadjf/
adetz*max(2*d00_err2(r-1),cij_err2(r-2)))/(2*(r-1))
3046 dij_err2(r) = max(
maxzadjf*dij_err2(r-1), &
3050 write(*,*)
'Dij_err(r)', r,dij_err(r),d00_err(r)
3055 write(*,*)
'Dij_err2(r)', r,dij_err2(r),d00_err2(r)
3072 d(n0,n1,n2,n3) = (c_0(n0-1,n1,n2,n3) + 2*
mm02*d(n0-1,n1,n2,n3) &
3073 + 4*duv(n0,n1,n2,n3) &
3074 +
f(1)*d(n0-1,n1+1,n2,n3) +
f(2)*d(n0-1,n1,n2+1,n3) &
3075 +
f(3)*d(n0-1,n1,n2,n3+1)) / (2*(r-1))
3082 write(*,*)
'CalcDpv1 Derrsym',derr
3083 write(*,*)
'CalcDpv1 Daccsym',derr/abs(d(0,0,0,0))
3088 write(*,*)
'CalcDpv1 Dijerr',dij_err(1:rmax)
3089 write(*,*)
'CalcDpv1 Dijacc',dij_err(1:rmax)/abs(d(0,0,0,0))
3091 write(*,*)
'CalcDpv1 Derr2sym',derr2
3092 write(*,*)
'CalcDpv1 Dacc2sym',derr2/abs(d(0,0,0,0))
3097 write(*,*)
'CalcDpv1 Dijerr2',dij_err2(1:rmax)
3098 write(*,*)
'CalcDpv1 Dijacc2',dij_err2(1:rmax)/abs(d(0,0,0,0))
3102 derr2 = max(derr2,dij_err2(0:rmax))
3103 derr = max(derr,dij_err(0:rmax))
3106 write(*,*)
'CalcDpv1 D(0,0,0,0) = ',d(0,0,0,0)
3108 write(*,*)
'CalcDpv1 D(0,1,1,1) = ',d(0,1,1,1)
3111 write(*,*)
'CalcDpv1 Derr',derr
3112 write(*,*)
'CalcDpv1 Dacc',derr/abs(d(0,0,0,0))
3113 write(*,*)
'CalcDpv1 Derr2',derr2
3114 write(*,*)
'CalcDpv1 Dacc2',derr2/abs(d(0,0,0,0))
3130 subroutine calcdpv1o(D,Duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,Derr,Derr2)
3134 integer,
intent(in) :: rmax,id
3135 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
3136 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
3137 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
3138 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
3139 double complex,
allocatable :: C_0(:,:,:,:), Cuv_0(:,:,:,:)
3140 double complex,
allocatable :: C_i(:,:,:,:), Cuv_i(:,:,:,:)
3141 double complex,
allocatable :: D_alt(:,:,:,:)
3142 double precision,
allocatable :: Cerr_i(:,:),Cerr2_i(:,:)
3143 double complex :: Smod(3)
3144 double complex :: D0_coli, elimminf2_coli
3145 double precision,
allocatable :: D00_err(:),Dij_err(:),Cij_err(:)
3146 double precision,
allocatable :: D00_err2(:),Dij_err2(:),Cij_err2(:)
3147 integer :: rmaxC,r,n0,n1,n2,n3,nn0,nn1,nn2,nn3,i,j
3148 integer :: bin,k,nid(0:3)
3153 d(0,0,0,0) = d0_coli(p10,p21,p32,p30,p20,p31,m02,m12,m22,m32)
3157 derr(0) = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
3158 derr2(0) = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
3160 if (rmax.eq.0)
return
3165 allocate(c_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
3166 allocate(cuv_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
3167 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
3168 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
3169 allocate(cerr_i(0:rmaxc,0:3))
3170 allocate(cerr2_i(0:rmaxc,0:3))
3173 allocate(d00_err(0:rmax))
3174 allocate(dij_err(0:rmax))
3175 allocate(cij_err(0:rmaxc))
3177 allocate(d00_err2(0:rmax))
3178 allocate(dij_err2(0:rmax))
3179 allocate(cij_err2(0:rmaxc))
3185 if (mod(id/bin,2).eq.0)
then
3193 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0))
3194 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1))
3195 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2))
3196 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3))
3199 write(*,*)
'CalcDpv1o Cerr_i=',cerr_i(:,0)
3200 write(*,*)
'CalcDpv1o Cerr_i=',cerr_i(:,1)
3201 write(*,*)
'CalcDpv1o Cerr_i=',cerr_i(:,2)
3202 write(*,*)
'CalcDpv1o Cerr_i=',cerr_i(:,3)
3209 n0 = (rmaxc-n1-n2-n3)
3210 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
3211 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
3212 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
3213 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
3276 dij_err(0) = derr(0)
3277 cij_err = max(cerr_i(:,0),cerr_i(:,1),cerr_i(:,2),cerr_i(:,3))
3281 dij_err2(0) = derr2(0)
3282 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,1),cerr2_i(:,2),cerr2_i(:,3))
3285 write(*,*)
'CalcDpv1o Cij_err=',cij_err
3286 write(*,*)
'CalcDpv1o Dij_err(0)=',dij_err(0)
3289 allocate(d_alt(0:rmax,0:rmax,0:rmax,0:rmax))
3300 d(n0,n1,n2,n3) = (c_0(n0-1,n1,n2,n3) + 2*
mm02*d(n0-1,n1,n2,n3) + 4*duv(n0,n1,n2,n3) &
3301 +
f(1)*d(n0-1,n1+1,n2,n3) +
f(2)*d(n0-1,n1,n2+1,n3) &
3302 +
f(3)*d(n0-1,n1,n2,n3+1)) / (2*(r-1))
3320 else if (n2.ge.1)
then
3333 smod(i) = -c_0(n0,nn1,nn2,nn3)-
f(i)*d(n0,nn1,nn2,nn3)
3337 smod(1) = smod(1) - 2*nn1*d(n0+1,nn1-1,nn2,nn3)
3339 smod(1) = smod(1) + c_i(n0,nn2,nn3,1)
3343 smod(2) = smod(2) - 2*nn2*d(n0+1,nn1,nn2-1,nn3)
3345 smod(2) = smod(2) + c_i(n0,nn1,nn3,2)
3349 smod(3) = smod(3) - 2*nn3*d(n0+1,nn1,nn2,nn3-1)
3351 smod(3) = smod(3) + c_i(n0,nn1,nn2,3)
3354 d(n0,n1,n2,n3) =
zinv(1,j)*smod(1) +
zinv(2,j)*smod(2) &
3368 if (n1.ge.1.and.n2+n3.ge.1)
then
3383 smod(i) = -c_0(n0,nn1,nn2,nn3)-
f(i)*d(n0,nn1,nn2,nn3)
3387 smod(1) = smod(1) - 2*nn1*d(n0+1,nn1-1,nn2,nn3)
3389 smod(1) = smod(1) + c_i(n0,nn2,nn3,1)
3393 smod(2) = smod(2) - 2*nn2*d(n0+1,nn1,nn2-1,nn3)
3395 smod(2) = smod(2) + c_i(n0,nn1,nn3,2)
3399 smod(3) = smod(3) - 2*nn3*d(n0+1,nn1,nn2,nn3-1)
3401 smod(3) = smod(3) + c_i(n0,nn1,nn2,3)
3404 d_alt(n0,n1,n2,n3) =
zinv(1,j)*smod(1) +
zinv(2,j)*smod(2) &
3407 derr(r)=max(derr(r),abs(d(n0,n1,n2,n3)-d_alt(n0,n1,n2,n3)))
3408 derr2(r)=max(derr2(r),abs(d(n0,n1,n2,n3)-d_alt(n0,n1,n2,n3)))
3423 d00_err(r) = max(2*abs(m02)*dij_err(r-2), cerr_i(r-2,0), &
3432 dij_err(r) = max(
maxzadjf*dij_err(r-1), &
3437 d00_err2(r) = max(2*abs(m02)*dij_err2(r-2), cerr2_i(r-2,0), &
3439 maxzadjf/
adetz*max(2*d00_err2(r-1),cij_err2(r-2)))/(2*(r-1))
3446 dij_err2(r) = max(
maxzadjf*dij_err2(r-1), &
3450 write(*,*)
'Dij_err(r)', r,dij_err(r),d00_err(r)
3455 write(*,*)
'Dij_err2(r)', r,dij_err2(r),d00_err2(r)
3472 d(n0,n1,n2,n3) = (c_0(n0-1,n1,n2,n3) + 2*
mm02*d(n0-1,n1,n2,n3) &
3473 + 4*duv(n0,n1,n2,n3) &
3474 +
f(1)*d(n0-1,n1+1,n2,n3) +
f(2)*d(n0-1,n1,n2+1,n3) &
3475 +
f(3)*d(n0-1,n1,n2,n3+1)) / (2*(r-1))
3482 write(*,*)
'CalcDpv1o Derrsym',derr
3483 write(*,*)
'CalcDpv1o Daccsym',derr/abs(d(0,0,0,0))
3488 write(*,*)
'CalcDpv1o Dijerr',dij_err(1:rmax)
3489 write(*,*)
'CalcDpv1o Dijacc',dij_err(1:rmax)/abs(d(0,0,0,0))
3491 write(*,*)
'CalcDpv1o Derr2sym',derr2
3492 write(*,*)
'CalcDpv1o Dacc2sym',derr2/abs(d(0,0,0,0))
3497 write(*,*)
'CalcDpv1o Dijerr2',dij_err2(1:rmax)
3498 write(*,*)
'CalcDpv1o Dijacc2',dij_err2(1:rmax)/abs(d(0,0,0,0))
3502 derr2 = max(derr2,dij_err2(0:rmax))
3503 derr = max(derr,dij_err(0:rmax))
3506 write(*,*)
'CalcDpv1o D(0,0,0,0) = ',d(0,0,0,0)
3508 write(*,*)
'CalcDpv1o D(0,1,1,1) = ',d(0,1,1,1)
3511 write(*,*)
'CalcDpv1o Derr',derr
3512 write(*,*)
'CalcDpv1o Dacc',derr/abs(d(0,0,0,0))
3513 write(*,*)
'CalcDpv1o Derr2',derr2
3514 write(*,*)
'CalcDpv1o Dacc2',derr2/abs(d(0,0,0,0))
3531 subroutine calcdpv(D,Duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,Derr,Derr2)
3535 integer,
intent(in) :: rmax,id
3536 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
3537 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
3538 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
3539 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
3540 double complex,
allocatable :: C_0(:,:,:,:), Cuv_0(:,:,:,:)
3541 double complex,
allocatable :: C_i(:,:,:,:), Cuv_i(:,:,:,:)
3542 double complex,
allocatable :: D_alt(:,:,:,:)
3543 double precision,
allocatable :: Cerr_i(:,:),Cerr2_i(:,:)
3544 double complex :: Smod(3)
3545 double complex :: D0_coli, elimminf2_coli
3546 double precision,
allocatable :: D00_err(:),Dij_err(:),Cij_err(:)
3547 double precision,
allocatable :: D00_err2(:),Dij_err2(:),Cij_err2(:)
3548 integer :: rmaxC,r,n0,n1,n2,n3,nn0,nn1,nn2,nn3,i,j
3549 integer :: bin,k,nid(0:3)
3554 d(0,0,0,0) = d0_coli(p10,p21,p32,p30,p20,p31,m02,m12,m22,m32)
3558 derr(0) = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
3559 derr2(0) = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
3561 if (rmax.eq.0)
return
3566 allocate(c_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
3567 allocate(cuv_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
3568 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
3569 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
3570 allocate(cerr_i(0:rmaxc,0:3))
3571 allocate(cerr2_i(0:rmaxc,0:3))
3574 allocate(d00_err(0:rmax))
3575 allocate(dij_err(0:rmax))
3576 allocate(cij_err(0:rmaxc))
3578 allocate(d00_err2(0:rmax))
3579 allocate(dij_err2(0:rmax))
3580 allocate(cij_err2(0:rmaxc))
3586 if (mod(id/bin,2).eq.0)
then
3594 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0))
3595 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1))
3596 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2))
3597 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3))
3600 write(*,*)
'CalcDpv Cerr_i=',cerr_i(:,0)
3601 write(*,*)
'CalcDpv Cerr_i=',cerr_i(:,1)
3602 write(*,*)
'CalcDpv Cerr_i=',cerr_i(:,2)
3603 write(*,*)
'CalcDpv Cerr_i=',cerr_i(:,3)
3610 n0 = (rmaxc-n1-n2-n3)
3611 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
3612 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
3613 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
3614 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
3677 dij_err(0) = derr(0)
3678 cij_err = max(cerr_i(:,0),cerr_i(:,1),cerr_i(:,2),cerr_i(:,3))
3682 dij_err2(0) = derr2(0)
3683 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,1),cerr2_i(:,2),cerr2_i(:,3))
3686 write(*,*)
'CalcDpv Cij_err=',cij_err
3687 write(*,*)
'CalcDpv Dij_err(0)=',dij_err(0)
3690 allocate(d_alt(0:rmax,0:rmax,0:rmax,0:rmax))
3695 if (mod(r,2).eq.0)
then
3698 d(n0,0,0,0) = (c_0(n0-1,0,0,0) + 2*
mm02*d(n0-1,0,0,0) + 4*duv(n0,0,0,0) &
3699 +
f(1)*d(n0-1,1,0,0) +
f(2)*d(n0-1,0,1,0) &
3700 +
f(3)*d(n0-1,0,0,1)) / (2*(r-1))
3714 else if (n2.ge.1)
then
3727 smod(i) = -c_0(n0,nn1,nn2,nn3)-
f(i)*d(n0,nn1,nn2,nn3)
3731 smod(1) = smod(1) - 2*nn1*d(n0+1,nn1-1,nn2,nn3)
3733 smod(1) = smod(1) + c_i(n0,nn2,nn3,1)
3737 smod(2) = smod(2) - 2*nn2*d(n0+1,nn1,nn2-1,nn3)
3739 smod(2) = smod(2) + c_i(n0,nn1,nn3,2)
3743 smod(3) = smod(3) - 2*nn3*d(n0+1,nn1,nn2,nn3-1)
3745 smod(3) = smod(3) + c_i(n0,nn1,nn2,3)
3748 d(n0,n1,n2,n3) =
zinv(1,j)*smod(1) +
zinv(2,j)*smod(2) &
3762 if (n1.ge.1.and.n2+n3.ge.1)
then
3777 smod(i) = -c_0(n0,nn1,nn2,nn3)-
f(i)*d(n0,nn1,nn2,nn3)
3781 smod(1) = smod(1) - 2*nn1*d(n0+1,nn1-1,nn2,nn3)
3783 smod(1) = smod(1) + c_i(n0,nn2,nn3,1)
3787 smod(2) = smod(2) - 2*nn2*d(n0+1,nn1,nn2-1,nn3)
3789 smod(2) = smod(2) + c_i(n0,nn1,nn3,2)
3793 smod(3) = smod(3) - 2*nn3*d(n0+1,nn1,nn2,nn3-1)
3795 smod(3) = smod(3) + c_i(n0,nn1,nn2,3)
3798 d_alt(n0,n1,n2,n3) =
zinv(1,j)*smod(1) +
zinv(2,j)*smod(2) &
3801 derr(r)=max(derr(r),abs(d(n0,n1,n2,n3)-d_alt(n0,n1,n2,n3)))
3802 derr2(r)=max(derr2(r),abs(d(n0,n1,n2,n3)-d_alt(n0,n1,n2,n3)))
3813 d00_err(r) = max(abs(m02)*dij_err(r-2), cerr_i(r-2,0), &
3819 dij_err(r) = max(
maxzadjf*dij_err(r-1), &
3823 d00_err2(r) = max(abs(m02)*dij_err2(r-2), cerr2_i(r-2,0), &
3829 dij_err2(r) = max(
maxzadjf*dij_err2(r-1), &
3841 d(n0,n1,n2,n3) = (c_0(n0-1,n1,n2,n3) + 2*
mm02*d(n0-1,n1,n2,n3) &
3842 + 4*duv(n0,n1,n2,n3) &
3843 +
f(1)*d(n0-1,n1+1,n2,n3) +
f(2)*d(n0-1,n1,n2+1,n3) &
3844 +
f(3)*d(n0-1,n1,n2,n3+1)) / (2*(r-1))
3851 write(*,*)
'CalcDpv Derrsym',derr
3852 write(*,*)
'CalcDpv Daccsym',derr/abs(d(0,0,0,0))
3858 write(*,*)
'CalcDpv Dijerr',dij_err(1:rmax)
3859 write(*,*)
'CalcDpv Dijacc',dij_err(1:rmax)/abs(d(0,0,0,0))
3862 derr2 = max(derr2,dij_err2(0:rmax))
3863 derr = max(derr,dij_err(0:rmax))
3866 write(*,*)
'CalcDpv Derr',derr
3867 write(*,*)
'CalcDpv Dacc',derr/abs(d(0,0,0,0))
3871 write(*,*)
'CalcDpv Derr ',derr
3872 write(*,*)
'CalcDpv Derr2',derr2
3886 subroutine calcdpv2(D,Duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,Derr,Derr2)
3890 integer,
intent(in) :: rmax,id
3891 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
3892 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
3893 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
3894 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
3895 double complex,
allocatable :: C_0(:,:,:,:), Cuv_0(:,:,:,:)
3896 double complex,
allocatable :: C_i(:,:,:,:), Cuv_i(:,:,:,:)
3897 double complex,
allocatable :: D_alt(:,:,:,:)
3898 double precision,
allocatable :: Cerr_i(:,:),Cerr2_i(:,:)
3899 double complex :: D0_coli, elimminf2_coli
3900 double complex :: Daux(1:rmax/2+1,0:rmax-1,0:rmax-1,0:rmax-1), Smod(3)
3901 double precision,
allocatable :: D00_err(:),Dij_err(:),Cij_err(:)
3902 double precision,
allocatable :: D00_err2(:),Dij_err2(:),Cij_err2(:)
3903 integer :: rmaxC,r,n0,n1,n2,n3,k
3904 integer :: bin,nid(0:3)
3907 write(*,*)
'CalcDpv2 in'
3912 d(0,0,0,0) = d0_coli(p10,p21,p32,p30,p20,p31,m02,m12,m22,m32)
3916 derr(0) = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
3917 derr2(0) = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
3919 if (rmax.eq.0)
return
3924 allocate(c_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
3925 allocate(cuv_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
3926 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
3927 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
3928 allocate(cerr_i(0:rmaxc,0:3))
3929 allocate(cerr2_i(0:rmaxc,0:3))
3932 allocate(d00_err(0:rmax+1))
3933 allocate(dij_err(0:rmax))
3934 allocate(cij_err(0:rmaxc))
3936 allocate(d00_err2(0:rmax+1))
3937 allocate(dij_err2(0:rmax))
3938 allocate(cij_err2(0:rmaxc))
3945 if (mod(id/bin,2).eq.0)
then
3952 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0))
3953 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1))
3954 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2))
3955 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3))
3961 n0 = (rmaxc-n1-n2-n3)
3962 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
3963 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
3964 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
3965 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
3997 dij_err(0) = derr(0)
3998 cij_err = max(cerr_i(:,0),cerr_i(:,1),cerr_i(:,2),cerr_i(:,3))
4002 dij_err2(0) = derr2(0)
4003 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,1),cerr2_i(:,2),cerr2_i(:,3))
4016 write(*,*)
'CalcDpv2 Dij_err(0)=',dij_err(0)
4017 write(*,*)
'CalcDpv2 Dij_acc(0)=',dij_err(0)/d(0,0,0,0)
4018 write(*,*)
'CalcDpv2 Cij_err=',cij_err
4021 allocate(d_alt(0:rmax,0:rmax,0:rmax,0:rmax))
4032 smod(k) = -c_0(n0-1,n1,n2,n3)
4036 smod(1) = smod(1) - 2*n1*d(n0,n1-1,n2,n3)
4038 smod(1) = smod(1) + c_i(n0-1,n2,n3,1)
4042 smod(2) = smod(2) - 2*n2*d(n0,n1,n2-1,n3)
4044 smod(2) = smod(2) + c_i(n0-1,n1,n3,2)
4048 smod(3) = smod(3) - 2*n3*d(n0,n1,n2,n3-1)
4050 smod(3) = smod(3) + c_i(n0-1,n1,n2,3)
4053 daux(n0,n1,n2,n3) = (d(n0-1,n1,n2,n3) -
mxinv(1,0)*smod(1) &
4066 d(n0,n1,n2,n3) = (daux(n0,n1,n2,n3) + 4d0*duv(n0,n1,n2,n3) &
4067 + c_0(n0-1,n1,n2,n3))/(r-1)/2d0
4121 smod(k) = -c_0(0,n1,n2,n3)
4125 smod(1) = smod(1) - 2*n1*d(1,n1-1,n2,n3)
4127 smod(1) = smod(1) + c_i(0,n2,n3,1)
4131 smod(2) = smod(2) - 2*n2*d(1,n1,n2-1,n3)
4133 smod(2) = smod(2) + c_i(0,n1,n3,2)
4137 smod(3) = smod(3) - 2*n3*d(1,n1,n2,n3-1)
4139 smod(3) = smod(3) + c_i(0,n1,n2,3)
4142 daux(1,n1,n2,n3) = (d(0,n1,n2,n3) -
mxinv(1,0)*smod(1) &
4145 d(0,n1+1,n2,n3) =
mxinv(0,1)*daux(1,n1,n2,n3) &
4147 d(0,n1,n2+1,n3) =
mxinv(0,2)*daux(1,n1,n2,n3) &
4149 d_alt(0,n1,n2,n3+1) =
mxinv(0,3)*daux(1,n1,n2,n3) &
4153 d(0,0,0,r) = d_alt(0,0,0,r)
4158 derr(r)=max(derr(r),abs(d(0,n1,n2,n3+1)-d_alt(0,n1,n2,n3+1)))
4159 derr2(r)=max(derr2(r),abs(d(0,n1,n2,n3+1)-d_alt(0,n1,n2,n3+1)))
4169 d00_err(r+1) = max(cerr_i(r-1,0),
adetx/
adetz*dij_err(r-1), &
4171 dij_err(r) = max(
maxzadjf*max(2*r*d00_err(r+1),cerr_i(r-1,0)), &
4174 d00_err2(r+1) = max(cerr2_i(r-1,0),
adetx/
adetz*dij_err2(r-1), &
4176 dij_err2(r) = max(
maxzadjf*max(2*r*d00_err2(r+1),cerr2_i(r-1,0)), &
4180 write(*,*)
'CalcDpv2 Cerr_i ',r-1, cerr_i(r-1,0)
4181 write(*,*)
'CalcDpv2 Cij_err ',r-1, cij_err(r-1)
4182 write(*,*)
'CalcDpv2 D00_err ',r+1, d00_err(r+1)
4183 write(*,*)
'CalcDpv2 Dij_err ',r, dij_err(r)
4184 write(*,*)
'CalcDpv2 Cerr2_i ',r-1, cerr2_i(r-1,0)
4185 write(*,*)
'CalcDpv2 Cij_err2',r-1, cij_err2(r-1)
4186 write(*,*)
'CalcDpv2 D00_err2',r+1, d00_err2(r+1)
4187 write(*,*)
'CalcDpv2 Dij_err2',r, dij_err2(r)
4203 d(n0,n1,n2,n3) = (c_0(n0-1,n1,n2,n3) + 2*
mm02*d(n0-1,n1,n2,n3) &
4204 + 4*duv(n0,n1,n2,n3) &
4205 +
f(1)*d(n0-1,n1+1,n2,n3) +
f(2)*d(n0-1,n1,n2+1,n3) &
4206 +
f(3)*d(n0-1,n1,n2,n3+1)) / (2*(r-1))
4220 do n0=max(2,r-rmax),r/2
4226 smod(k) = -c_0(n0-1,n1,n2,n3)
4230 smod(1) = smod(1) - 2*n1*d(n0,n1-1,n2,n3)
4232 smod(1) = smod(1) + c_i(n0-1,n2,n3,1)
4236 smod(2) = smod(2) - 2*n2*d(n0,n1,n2-1,n3)
4238 smod(2) = smod(2) + c_i(n0-1,n1,n3,2)
4242 smod(3) = smod(3) - 2*n3*d(n0,n1,n2,n3-1)
4244 smod(3) = smod(3) + c_i(n0-1,n1,n2,3)
4247 daux(n0,n1,n2,n3) = (d(n0-1,n1,n2,n3) -
mxinv(1,0)*smod(1) &
4260 d(n0,n1,n2,n3) = (daux(n0,n1,n2,n3) + 4d0*duv(n0,n1,n2,n3) &
4261 + c_0(n0-1,n1,n2,n3))/(r-1)/2d0
4275 write(*,*)
'CalcDpv2 Derrsym',derr
4276 write(*,*)
'CalcDpv2 Daccsym',derr/abs(d(0,0,0,0))
4277 write(*,*)
'CalcDpv2 Derr2sym',derr2
4278 write(*,*)
'CalcDpv2 Dacc2sym',derr2/abs(d(0,0,0,0))
4280 write(*,*)
'CalcDpv2 Dijerr',dij_err
4281 write(*,*)
'CalcDpv2 Dijacc',dij_err/abs(d(0,0,0,0))
4284 derr2 = max(derr2,dij_err2(0:rmax))
4285 derr = max(derr,dij_err(0:rmax))
4296 write(*,*)
'CalcDpv2 Derr ',derr
4297 write(*,*)
'CalcDpv2 Dacc ',derr/abs(d(0,0,0,0))
4298 write(*,*)
'CalcDpv2 Derr2',derr2
4299 write(*,*)
'CalcDpv2 Dacc2',derr2/abs(d(0,0,0,0))
4317 subroutine calcdg(D,Duv,p10,p21,p32,p30,p20,p31, &
4318 m02,m12,m22,m32,rmax,ordg_min,ordg_max,id,Derr,Derr2)
4322 integer,
intent(in) :: rmax,ordg_min,ordg_max,id
4323 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
4324 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
4325 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
4326 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
4327 double complex :: Zadjfj,Zadj2(4), Zadjkl, Xtilde
4328 double complex,
allocatable :: Dexpg(:,:,:,:,:), DuvExpg(:,:,:,:)
4329 double complex,
allocatable :: C_0(:,:,:,:), Cuv_0(:,:,:,:), Shat(:,:,:,:,:)
4330 double complex,
allocatable :: C_i(:,:,:,:), Cuv_i(:,:,:,:)
4331 double complex,
allocatable :: D_alt(:,:,:,:)
4332 double precision,
allocatable :: Cerr_i(:,:),Cerr2_i(:,:)
4333 double complex :: Smod(3), Skl, DexpgAux
4334 double complex :: cC0f, elimminf2_coli
4335 double precision,
allocatable :: D00_err(:),Dij_err(:),Cij_err(:),acc_req_Cextra(:)
4336 double precision,
allocatable :: D00_err2(:),Dij_err2(:),Cij_err2(:)
4337 double precision :: maxDexpg(0:1,0:rmax+ordg_min+1,0:ordg_max),truncfacexp
4338 integer :: rmaxC,rmaxExp,gtrunc,r,n0,n1,n2,n3,k,l,i,j,m,n,g,rg
4339 integer :: inds0(3), inds(3), inds2(2,4)
4340 integer :: bin,nid(0:3)
4341 logical :: errorwriteflag
4344 write(*,*)
'CalcDg in, ord',rmax,ordg_min,ordg_max
4348 rmaxc = rmax + ordg_min
4349 allocate(c_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
4350 allocate(cuv_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
4351 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
4352 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
4353 allocate(cerr_i(0:rmaxc,0:3))
4354 allocate(cerr2_i(0:rmaxc,0:3))
4355 allocate(acc_req_cextra(0:rmaxc))
4361 if (mod(id/bin,2).eq.0)
then
4370 acc_req_cextra(0:rmax) = acc_req_cind
4371 if (
x_g.ne.0d0)
then
4373 acc_req_cextra(r)= acc_req_cextra(r-1)/
x_g
4376 acc_req_cextra(rmax+1:rmaxc) = acc_inf
4379 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0),rmax,acc_req_cextra)
4380 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1),rmax,acc_req_cextra)
4381 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2),rmax,acc_req_cextra)
4382 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3),rmax,acc_req_cextra)
4388 n0 = (rmaxc-n1-n2-n3)
4389 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
4390 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
4391 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
4392 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
4446 allocate(shat(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc,3))
4454 shat(n0,n1,n2,n3,:) = -c_0(n0,n1,n2,n3)
4457 shat(n0,n1,n2,n3,1) = shat(n0,n1,n2,n3,1) + c_i(n0,n2,n3,1)
4461 shat(n0,n1,n2,n3,2) = shat(n0,n1,n2,n3,2) + c_i(n0,n1,n3,2)
4465 shat(n0,n1,n2,n3,3) = shat(n0,n1,n2,n3,3) + c_i(n0,n1,n2,3)
4488 inds2 = reshape((/2,2,2,3,3,2,3,3/),shape(inds2))
4498 inds2 = reshape((/1,1,1,3,3,1,3,3/),shape(inds2))
4508 inds2 = reshape((/1,1,1,2,2,1,2,2/),shape(inds2))
4518 inds2 = reshape((/2,1,2,3,3,1,3,3/),shape(inds2))
4528 inds2 = reshape((/2,1,2,2,3,1,3,2/),shape(inds2))
4537 inds2 = reshape((/1,1,1,2,3,1,3,2/),shape(inds2))
4553 allocate(dexpg(0:rmaxexp/2,0:rmaxexp,0:rmaxexp,0:rmaxexp,0:ordg_max))
4557 allocate(duvexpg(0:rmaxexp,0:rmaxexp,0:rmaxexp,0:rmaxexp))
4559 duv(0:rmax,0:rmax,0:rmax,0:rmax) = duvexpg(0:rmax,0:rmax,0:rmax,0:rmax)
4562 allocate(d00_err(0:rmaxexp))
4563 allocate(dij_err(0:rmaxexp))
4564 allocate(cij_err(0:rmaxc))
4566 allocate(d00_err2(0:rmaxexp))
4567 allocate(dij_err2(0:rmaxexp))
4568 allocate(cij_err2(0:rmaxc))
4574 cij_err = max(cerr_i(:,0),cerr_i(:,1),cerr_i(:,2),cerr_i(:,3))
4579 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,1),cerr2_i(:,2),cerr2_i(:,3))
4582 write(*,*)
'CalcDg Cij_err = ',cij_err
4583 write(*,*)
'CalcDg C0_err = ', cerr_i(0,0),cerr_i(0,1),cerr_i(0,2),cerr_i(0,3)
4584 write(*,*)
'CalcDg C0 = ', c_0(0,0,0,0),c_i(0,0,0,1),c_i(0,0,0,2),c_i(0,0,0,3)
4596 rloop:
do r=1,rmaxexp
4602 if (r.gt.rmax+gtrunc+1)
exit rloop
4605 write(*,*)
'CalcDg rloop',r
4625 dexpgaux = 2d0*
zadj(k,l)*c_0(n0-1,n1,n2,n3) &
4626 + xtilde*dexpg(n0-1,n1,n2,n3,0) &
4627 + 4d0*
zadj(k,l)*duvexpg(n0,n1,n2,n3)
4632 dexpgaux = dexpgaux +
zadj(i,l)*shat(n0-1,inds(1),inds(2),inds(3),i)
4638 dexpgaux = dexpgaux -
zadj(k,l)*shat(n0-1,inds(1),inds(2),inds(3),i)
4645 skl =
f(n)*shat(n0-1,inds0(1),inds0(2),inds0(3),m)
4648 if (inds(m).ge.1)
then
4650 skl = skl - 2d0*
f(n)*inds0(m)*dexpg(n0,inds(1),inds(2),inds(3),0)
4651 if (inds(n).ge.1)
then
4653 skl = skl - 4d0*inds0(m)*(inds(n)+1)*dexpg(n0+1,inds(1),inds(2),inds(3),0)
4657 if (inds(n).ge.1)
then
4659 skl = skl + 2d0*inds0(n)*shat(n0,inds(1),inds(2),inds(3),m) &
4660 - 2d0*
f(m)*inds0(n)*dexpg(n0,inds(1),inds(2),inds(3),0)
4663 dexpgaux = dexpgaux - zadj2(i)*skl
4667 dexpg(n0,n1,n2,n3,0) = dexpgaux/(2d0*zadjkl)/(2d0*(r-n0))
4670 maxdexpg(1,r,0) = maxdexpg(1,r,0) + abs(dexpg(n0,n1,n2,n3,0) )
4673 if (r-n0.le.rmax)
then
4674 d(n0,n1,n2,n3) = dexpg(n0,n1,n2,n3,0)
4684 maxdexpg(0,r-1,0)=0d0
4689 smod = shat(0,n1,n2,n3,:)
4691 smod(1) = smod(1) - 2d0*n1*dexpg(1,n1-1,n2,n3,0)
4694 smod(2) = smod(2) - 2d0*n2*dexpg(1,n1,n2-1,n3,0)
4697 smod(3) = smod(3) - 2d0*n3*dexpg(1,n1,n2,n3-1,0)
4700 dexpg(0,n1,n2,n3,0) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
4701 +
zadj(3,j)*smod(3))/zadjfj
4702 maxdexpg(0,r-1,0) = maxdexpg(0,r-1,0) + abs(dexpg(0,n1,n2,n3,0))
4703 if (r.le.rmax+1)
then
4704 d(0,n1,n2,n3) = dexpg(0,n1,n2,n3,0)
4709 if(n0.eq.0.and.n1.eq.0.and.n2.eq.0.and.n3.eq.0)
then
4710 write(*,*)
'D2(0,0,0,0)= ',0,d(n0,n1,n2,n3),
detz/zadjfj
4743 write(*,*)
'CalcDg maxDexpg 0',r-1, maxdexpg(0,r-1,0)
4746 if(r.le.rmax+1)
then
4748 derr(r-1) =
fac_g*maxdexpg(0,r-1,0)
4754 d00_err(r) = max(cij_err(r-1),cij_err(r-2), &
4758 dij_err(r-1)=
maxzadj*max(cij_err(r-1),2*d00_err(r))/abs(zadjfj)
4761 d00_err2(r) = max(cij_err2(r-1),cij_err2(r-2), &
4766 dij_err2(r-1)=
maxzadj*max(cij_err2(r-1),2*d00_err2(r))/abs(zadjfj)
4774 gloop:
do g=1,min(gtrunc,r-1)
4778 write(*,*)
'gloop ',g,rg
4784 maxdexpg(1,rg,g) = 0d0
4797 dexpgaux = xtilde*dexpg(n0-1,n1,n2,n3,g) &
4798 -
detz*dexpg(n0-1,inds(1),inds(2),inds(3),g-1)
4808 if (inds(m).ge.1)
then
4810 skl = skl - 2d0*
f(n)*inds0(m)*dexpg(n0,inds(1),inds(2),inds(3),g)
4811 if (inds(n).ge.1)
then
4813 skl = skl - 4d0*inds0(m)*(inds(n)+1)*dexpg(n0+1,inds(1),inds(2),inds(3),g)
4817 if (inds(n).ge.1)
then
4819 skl = skl - 2d0*
f(m)*inds0(n)*dexpg(n0,inds(1),inds(2),inds(3),g)
4822 dexpgaux = dexpgaux - zadj2(i)*skl
4826 dexpg(n0,n1,n2,n3,g) = dexpgaux/(2d0*zadjkl)/(2d0*(rg-n0))
4830 maxdexpg(1,rg,g) = maxdexpg(1,rg,g) + abs(dexpg(n0,n1,n2,n3,g))
4832 if (g.eq.1.and.abs(dexpg(1,n1,n2,n3,g)).gt. &
4833 truncfacexp*max(1/
m2scale,maxdexpg(1,rg,g-1)) .or. &
4834 g.ge.2.and.abs(dexpg(1,n1,n2,n3,g)).gt. &
4835 truncfacexp*maxdexpg(1,rg,g-1))
then
4838 write(*,*)
'CalcDg exit gloop',n0,n1,n2,n3,g,abs(dexpg(n0,n1,n2,n3,g)),maxdexpg(1,rg,g-1),truncfacexp
4852 if (rg-n0.le.rmax)
then
4856 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) + dexpg(n0,n1,n2,n3,g)
4870 maxdexpg(0,rg-1,g) = 0d0
4877 smod(1) = smod(1) - 2d0*n1*dexpg(1,n1-1,n2,n3,g)
4880 smod(2) = smod(2) - 2d0*n2*dexpg(1,n1,n2-1,n3,g)
4883 smod(3) = smod(3) - 2d0*n3*dexpg(1,n1,n2,n3-1,g)
4890 dexpg(0,n1,n2,n3,g) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
4891 +
zadj(3,j)*smod(3) &
4892 -
detz*dexpg(0,inds(1),inds(2),inds(3),g-1))/zadjfj
4894 maxdexpg(0,rg-1,g) = maxdexpg(0,rg-1,g) + abs(dexpg(0,n1,n2,n3,g))
4905 if (g.eq.1.and.abs(dexpg(0,n1,n2,n3,g)).gt. &
4906 truncfacexp*max(1/
m2scale**2,maxdexpg(0,rg-1,g-1)) .or. &
4907 g.ge.2.and.abs(dexpg(0,n1,n2,n3,g)).gt. &
4908 truncfacexp*maxdexpg(0,rg-1,g-1))
then
4911 write(*,*)
'CalcDg exit gloop',0,n1,n2,n3,g,abs(dexpg(0,n1,n2,n3,g)),maxdexpg(0,rg-1,g-1),truncfacexp
4928 d00_err(rg) = max( d00_err(rg), &
4929 max( abs(
detz)*dij_err(rg),abs(xtilde)*dij_err(rg-2), &
4930 maxzadj2f*d00_err(rg-1) ) / abs(zadjkl) &
4933 dij_err(rg-1)=max(dij_err(rg-1), &
4934 max(2*
maxzadj*d00_err(rg),abs(
detz)*dij_err(rg))/abs(zadjfj) )
4937 d00_err2(rg) = max( d00_err2(rg), &
4938 max( abs(
detz)*dij_err2(rg),abs(xtilde)*dij_err2(rg-2), &
4939 maxzadj2f*d00_err2(rg-1) ) / abs(zadjkl) &
4942 dij_err2(rg-1)=max(dij_err2(rg-1), &
4943 max(2*
maxzadj*d00_err2(rg),abs(
detz)*dij_err2(rg))/abs(zadjfj) )
4947 if (rg-n0.le.rmax)
then
4951 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) + dexpg(n0,n1,n2,n3,g)
4958 if ((rg.le.rmax+1))
then
4963 d(0,n1,n2,n3) = d(0,n1,n2,n3) + dexpg(0,n1,n2,n3,g)
4965 if(abs(dexpg(0,n1,n2,n3,g-1)).ne.0d0)
then
4966 derr(rg-1)=max(derr(rg-1),abs(dexpg(0,n1,n2,n3,g))*min(1d0,abs(dexpg(0,n1,n2,n3,g))/abs(dexpg(0,n1,n2,n3,g-1))))
4968 derr(rg-1)=max(derr(rg-1),abs(dexpg(0,n1,n2,n3,g)))
4980 if(dij_err2(rg-1).gt.3d0*derr(rg-1))
then
4982 if(dij_err(rg-1).gt.3d0*derr(rg-1))
then
4984 gtrunc = min(g,gtrunc)
4987 write(*,*)
'CalcDg exit err',r,rg-1,g,gtrunc,dij_err(rg-1),derr(rg-1)
4997 write(*,*)
'CalcDg D(0,0,0,0) = ',r,d(0,0,0,0)
5000 write(*,*)
'CalcDg D(0,0,1,0) = ',r,d(0,0,1,0)
5002 if(r.gt.2.and.rmax.ge.2)
then
5003 write(*,*)
'CalcDg D(1,0,0,0) = ',r,d(1,0,0,0)
5004 write(*,*)
'CalcDg D(0,2,0,0) = ',r,d(0,2,0,0)
5005 write(*,*)
'CalcDg D(0,0,1,1) = ',r,d(0,0,1,1)
5007 write(*,*)
'CalcDg D(0,0,2,0) = ',r,d(0,0,2,0)
5009 if(r.gt.3.and.rmax.ge.3)
then
5010 write(*,*)
'CalcDg D(1,0,1,0) = ',r,d(1,0,1,0)
5011 write(*,*)
'CalcDg D(1,1,0,0) = ',r,d(1,1,0,0)
5012 write(*,*)
'CalcDg D(1,0,1,0) = ',r,d(1,0,1,0)
5013 write(*,*)
'CalcDg D(1,0,0,1) = ',r,d(1,0,0,1)
5015 write(*,*)
'CalcDg D(0,3,0,0) = ',r,d(0,3,0,0)
5016 write(*,*)
'CalcDg D(0,2,1,0) = ',r,d(0,2,1,0)
5017 write(*,*)
'CalcDg D(0,2,0,1) = ',r,d(0,2,0,1)
5018 write(*,*)
'CalcDg D(0,0,3,0) = ',r,d(0,0,3,0)
5019 write(*,*)
'CalcDg D(0,1,1,1) = ',r,d(0,1,1,1)
5020 write(*,*)
'CalcDg D(0,0,2,1) = ',r,d(0,0,2,1)
5022 write(*,*)
'CalcDg Dij_err',r,dij_err
5023 write(*,*)
'CalcDg Dij_acc',r,dij_err/abs(d(0,0,0,0))
5025 write(*,*)
'CalcDg err',r,derr
5026 write(*,*)
'CalcDg acc',r,derr/abs(d(0,0,0,0))
5029 derr2 = max(derr,dij_err2(0:rmax))
5030 derr = max(derr,dij_err(0:rmax))
5033 write(*,*)
'CalcDg exit r',r,maxval(derr),maxval(derr2),acc_req_d*abs(d(0,0,0,0))
5039 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0)
then
5047 d(n0,n1,n2,rg-2*n0-n1-n2)=0d0
5055 d(0,n1,n2,r-n1-n2)=0d0
5061 111
format(a22,2(
'(',g24.17,
',',g24.17,
') ':))
5062 call seterrflag_coli(-5)
5063 call errout_coli(
'CalcDg',
' exit rloop for D', &
5065 if (errorwriteflag)
then
5066 write(nerrout_coli,100)
' CalcDg: exit rloop for D ', &
5067 ' should not appear'
5068 write(nerrout_coli,111)
' CalcDg: p10 = ',p10
5069 write(nerrout_coli,111)
' CalcDg: p21 = ',p21
5070 write(nerrout_coli,111)
' CalcDg: p32 = ',p32
5071 write(nerrout_coli,111)
' CalcDg: p30 = ',p30
5072 write(nerrout_coli,111)
' CalcDg: p20 = ',p20
5073 write(nerrout_coli,111)
' CalcDg: p31 = ',p31
5074 write(nerrout_coli,111)
' CalcDg: m02 = ',m02
5075 write(nerrout_coli,111)
' CalcDg: m12 = ',m12
5076 write(nerrout_coli,111)
' CalcDg: m22 = ',m22
5077 write(nerrout_coli,111)
' CalcDg: m32 = ',m32
5082 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0.and.r.ge.rmax)
then
5100 write(*,*)
'CalcDg exp rmax+1',r,n0,n1,n2,n3, d(n0,n1,n2,n3)
5102 d(n0,n1,n2,n3) = (c_0(n0-1,n1,n2,n3) + 2*
mm02*d(n0-1,n1,n2,n3) &
5103 + 4*duv(n0,n1,n2,n3) &
5104 +
f(1)*d(n0-1,n1+1,n2,n3) +
f(2)*d(n0-1,n1,n2+1,n3) &
5105 +
f(3)*d(n0-1,n1,n2,n3+1)) / (2*(r-1))
5107 write(*,*)
'CalcDg dir rmax+1',r,n0,n1,n2,n3, d(n0,n1,n2,n3)
5116 write(*,*)
'CalcDg D(0,0,0,0) = ',d(0,0,0,0)
5118 write(*,*)
'CalcDg D(1,0,1,0) = ',d(1,0,1,0)
5121 write(*,*)
'CalcDg final err ',derr
5122 write(*,*)
'CalcDg final acc ',derr/abs(d(0,0,0,0))
5123 write(*,*)
'CalcDg final err2',derr2
5124 write(*,*)
'CalcDg final acc2',derr2/abs(d(0,0,0,0))
5140 subroutine calcdgm(D,Duv,p10,p21,p32,p30,p20,p31, &
5141 m02,m12,m22,m32,rmax,ordgm_min,ordgm_max,id,Derr,Derr2)
5145 integer,
intent(in) :: rmax,ordgm_min,ordgm_max,id
5146 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
5147 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
5148 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
5149 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
5150 double complex :: Zadjfj,Zadj2(4), Zadjkl, Xtilde
5151 double complex,
allocatable :: Dexpgm(:,:,:,:,:), DuvExpgm(:,:,:,:)
5152 double complex,
allocatable :: C_0(:,:,:,:), Cuv_0(:,:,:,:), Shat(:,:,:,:,:)
5153 double complex,
allocatable :: C_i(:,:,:,:), Cuv_i(:,:,:,:)
5154 double complex,
allocatable :: D_alt(:,:,:,:)
5155 double precision,
allocatable :: Cerr_i(:,:), Cerr2_i(:,:)
5156 double complex :: Smod(3), Skl, DexpgmAux
5157 double complex :: cC0f, elimminf2_coli
5158 double precision,
allocatable :: D00_err(:),Dij_err(:),Cij_err(:),acc_req_Cextra(:)
5159 double precision,
allocatable :: D00_err2(:),Dij_err2(:),Cij_err2(:)
5160 double precision :: maxDexpgm(0:1,0:rmax+ordgm_min+1,0:ordgm_max),truncfacexp
5161 integer :: rmaxC,rmaxExp,gtrunc,r,n0,n1,n2,n3,k,l,i,j,m,n,gm,rgm
5162 integer :: inds0(3), inds(3), inds2(2,4)
5163 integer :: bin,nid(0:3)
5164 logical :: errorwriteflag
5167 write(*,*)
'CalcDgm in, ord',rmax,ordgm_min,ordgm_max
5171 rmaxc = rmax + ordgm_min
5172 allocate(c_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
5173 allocate(cuv_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
5174 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
5175 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
5176 allocate(cerr_i(0:rmaxc,0:3))
5177 allocate(cerr2_i(0:rmaxc,0:3))
5178 allocate(acc_req_cextra(0:rmaxc))
5184 if (mod(id/bin,2).eq.0)
then
5193 acc_req_cextra(0:rmax) = acc_req_cind
5194 if (
x_gm.ne.0d0)
then
5196 acc_req_cextra(r)= acc_req_cextra(r-1)/
x_gm
5199 acc_req_cextra(rmax+1,rmaxc) = acc_inf
5202 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0),rmax,acc_req_cextra)
5203 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1),rmax,acc_req_cextra)
5204 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2),rmax,acc_req_cextra)
5205 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3),rmax,acc_req_cextra)
5211 n0 = (rmaxc-n1-n2-n3)
5212 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
5213 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
5214 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
5215 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
5221 allocate(shat(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc,3))
5229 shat(n0,n1,n2,n3,:) = -c_0(n0,n1,n2,n3)
5232 shat(n0,n1,n2,n3,1) = shat(n0,n1,n2,n3,1) + c_i(n0,n2,n3,1)
5236 shat(n0,n1,n2,n3,2) = shat(n0,n1,n2,n3,2) + c_i(n0,n1,n3,2)
5240 shat(n0,n1,n2,n3,3) = shat(n0,n1,n2,n3,3) + c_i(n0,n1,n2,3)
5265 inds2 = reshape((/2,2,2,3,3,2,3,3/),shape(inds2))
5275 inds2 = reshape((/1,1,1,3,3,1,3,3/),shape(inds2))
5285 inds2 = reshape((/1,1,1,2,2,1,2,2/),shape(inds2))
5295 inds2 = reshape((/2,1,2,3,3,1,3,3/),shape(inds2))
5305 inds2 = reshape((/2,1,2,2,3,1,3,2/),shape(inds2))
5314 inds2 = reshape((/1,1,1,2,3,1,3,2/),shape(inds2))
5330 allocate(dexpgm(0:rmaxexp/2,0:rmaxexp,0:rmaxexp,0:rmaxexp,0:ordgm_max))
5334 allocate(duvexpgm(0:rmaxexp,0:rmaxexp,0:rmaxexp,0:rmaxexp))
5336 duv(0:rmax,0:rmax,0:rmax,0:rmax) = duvexpgm(0:rmax,0:rmax,0:rmax,0:rmax)
5339 allocate(d00_err(0:rmaxexp))
5340 allocate(dij_err(0:rmaxexp))
5341 allocate(cij_err(0:rmaxc))
5343 allocate(d00_err2(0:rmaxexp))
5344 allocate(dij_err2(0:rmaxexp))
5345 allocate(cij_err2(0:rmaxc))
5351 cij_err = max(cerr_i(:,0),cerr_i(:,1),cerr_i(:,2),cerr_i(:,3))
5356 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,1),cerr2_i(:,2),cerr2_i(:,3))
5359 write(*,*)
'CalcDgm Cij_err = ',cij_err
5360 write(*,*)
'CalcDgm C0_err = ', cerr_i(0,0),cerr_i(0,1),cerr_i(0,2),cerr_i(0,3)
5361 write(*,*)
'CalcDgm C0 = ', c_0(0,0,0,0),c_i(0,0,0,1),c_i(0,0,0,2),c_i(0,0,0,3)
5373 rloop:
do r=1,rmaxexp
5379 if (r.gt.rmax+gtrunc+1)
exit rloop
5382 write(*,*)
'CalcDgm rloop',r
5392 maxdexpgm(1,r,0)=0d0
5402 dexpgmaux = 2d0*
zadj(k,l)*c_0(n0-1,n1,n2,n3) &
5403 + xtilde*dexpgm(n0-1,n1,n2,n3,0) &
5404 + 4d0*
zadj(k,l)*duvexpgm(n0,n1,n2,n3)
5409 dexpgmaux = dexpgmaux +
zadj(i,l)*shat(n0-1,inds(1),inds(2),inds(3),i)
5415 dexpgmaux = dexpgmaux -
zadj(k,l)*shat(n0-1,inds(1),inds(2),inds(3),i)
5422 skl =
f(n)*shat(n0-1,inds0(1),inds0(2),inds0(3),m)
5425 if (inds(m).ge.1)
then
5427 skl = skl - 2d0*
f(n)*inds0(m)*dexpgm(n0,inds(1),inds(2),inds(3),0)
5428 if (inds(n).ge.1)
then
5430 skl = skl - 4d0*inds0(m)*(inds(n)+1)*dexpgm(n0+1,inds(1),inds(2),inds(3),0)
5434 if (inds(n).ge.1)
then
5436 skl = skl + 2d0*inds0(n)*shat(n0,inds(1),inds(2),inds(3),m) &
5437 - 2d0*
f(m)*inds0(n)*dexpgm(n0,inds(1),inds(2),inds(3),0)
5440 dexpgmaux = dexpgmaux - zadj2(i)*skl
5444 dexpgm(n0,n1,n2,n3,0) = dexpgmaux/(2d0*zadjkl)/(2d0*(r-n0))
5447 maxdexpgm(1,r,0) = maxdexpgm(1,r,0) + abs(dexpgm(n0,n1,n2,n3,0) )
5450 if (r-n0.le.rmax)
then
5451 d(n0,n1,n2,n3) = dexpgm(n0,n1,n2,n3,0)
5461 maxdexpgm(0,r-1,0)=0d0
5466 smod = shat(0,n1,n2,n3,:)
5468 smod(1) = smod(1) - 2d0*n1*dexpgm(1,n1-1,n2,n3,0)
5471 smod(2) = smod(2) - 2d0*n2*dexpgm(1,n1,n2-1,n3,0)
5474 smod(3) = smod(3) - 2d0*n3*dexpgm(1,n1,n2,n3-1,0)
5477 dexpgm(0,n1,n2,n3,0) = (
zadjf(1)*smod(1) +
zadjf(2)*smod(2) &
5479 maxdexpgm(0,r-1,0) = maxdexpgm(0,r-1,0) + abs(dexpgm(0,n1,n2,n3,0))
5480 if (r.le.rmax+1)
then
5481 d(0,n1,n2,n3) = dexpgm(0,n1,n2,n3,0)
5486 if(n0.eq.0.and.n1.eq.0.and.n2.eq.0.and.n3.eq.0)
then
5487 write(*,*)
'D2(0,0,0,0)= ',0,d(n0,n1,n2,n3),
detz*
fmax/
zadjff
5489 write(*,*)
'D2(0,0,0,0)= ',smod
5492 write(*,*)
'D2(0,0,0,0)= ',smod(1)*
zadjf(1)/
zadjff, &
5504 if(r.le.rmax+1)
then
5506 derr(r-1) =
fac_gm*maxdexpgm(0,r-1,0)
5511 d00_err(r) = max(cij_err(r-1),cij_err(r-2), &
5518 d00_err2(r) = max(cij_err2(r-1),cij_err2(r-2), &
5522 dij_err2(r-1)=
maxzadjf*max(cij_err2(r-1),2*d00_err2(r))/abs(
zadjff)
5529 gloop:
do gm=1,min(gtrunc,r-1)
5537 maxdexpgm(1,rgm,gm) = 0d0
5550 dexpgmaux = xtilde*dexpgm(n0-1,n1,n2,n3,gm) &
5551 -
detz*dexpgm(n0-1,inds(1),inds(2),inds(3),gm-1)
5561 if (inds(m).ge.1)
then
5563 skl = skl - 2d0*
f(n)*inds0(m)*dexpgm(n0,inds(1),inds(2),inds(3),gm)
5564 if (inds(n).ge.1)
then
5566 skl = skl - 4d0*inds0(m)*(inds(n)+1)*dexpgm(n0+1,inds(1),inds(2),inds(3),gm)
5570 if (inds(n).ge.1)
then
5572 skl = skl - 2d0*
f(m)*inds0(n)*dexpgm(n0,inds(1),inds(2),inds(3),gm)
5575 dexpgmaux = dexpgmaux - zadj2(i)*skl
5579 dexpgm(n0,n1,n2,n3,gm) = dexpgmaux/(2d0*zadjkl)/(2d0*(rgm-n0))
5583 maxdexpgm(1,rgm,gm) = maxdexpgm(1,rgm,gm) + abs(dexpgm(n0,n1,n2,n3,gm))
5585 if (gm.eq.1.and.abs(dexpgm(1,n1,n2,n3,gm)).gt. &
5586 truncfacexp*max(1/
m2scale,maxdexpgm(1,rgm,gm-1)) .or. &
5587 gm.ge.2.and.abs(dexpgm(1,n1,n2,n3,gm)).gt. &
5588 truncfacexp*maxdexpgm(1,rgm,gm-1))
then
5591 write(*,*)
'CalcDgm exit gloop',n0,n1,n2,n3,gm,abs(dexpgm(n0,n1,n2,n3,gm)),maxdexpgm(1,rgm,gm-1),truncfacexp
5605 if (rgm-n0.le.rmax)
then
5609 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) + dexpgm(n0,n1,n2,n3,gm)
5623 maxdexpgm(0,rgm-1,gm) = 0d0
5630 smod(1) = smod(1) - 2d0*n1*dexpgm(1,n1-1,n2,n3,gm)
5633 smod(2) = smod(2) - 2d0*n2*dexpgm(1,n1,n2-1,n3,gm)
5636 smod(3) = smod(3) - 2d0*n3*dexpgm(1,n1,n2,n3-1,gm)
5643 dexpgm(0,n1,n2,n3,gm) = (
zadjf(1)*smod(1) +
zadjf(2)*smod(2) &
5644 +
zadjf(3)*smod(3) &
5646 f(1)*dexpgm(0,inds(1)+1,inds(2),inds(3),gm-1) &
5647 +
f(2)*dexpgm(0,inds(1),inds(2)+1,inds(3),gm-1) &
5648 +
f(3)*dexpgm(0,inds(1),inds(2),inds(3)+1,gm-1)) &
5651 maxdexpgm(0,rgm-1,gm) = maxdexpgm(0,rgm-1,gm) + abs(dexpgm(0,n1,n2,n3,gm))
5662 if (gm.eq.1.and.abs(dexpgm(0,n1,n2,n3,gm)).gt. &
5663 truncfacexp*max(1/
m2scale**2,maxdexpgm(0,rgm-1,gm-1)) .or. &
5664 gm.ge.2.and.abs(dexpgm(0,n1,n2,n3,gm)).gt. &
5665 truncfacexp*maxdexpgm(0,rgm-1,gm-1))
then
5668 write(*,*)
'CalcDgm exit gloop',0,n1,n2,n3,gm,abs(dexpgm(0,n1,n2,n3,gm)),maxdexpgm(0,rgm-1,gm-1),truncfacexp
5679 d00_err(rgm) = max( d00_err(rgm), &
5680 max( abs(m02)*dij_err(rgm-2), &
5681 max( abs(
detz)*dij_err(rgm),abs(xtilde)*dij_err(rgm-2), &
5682 maxzadj2f*d00_err(rgm-1) ) / abs(zadjkl) ) &
5685 dij_err(rgm-1)=max(dij_err(rgm-1), &
5688 d00_err2(rgm) = max( d00_err2(rgm), &
5689 max( abs(m02)*dij_err2(rgm-2), &
5690 max( abs(
detz)*dij_err2(rgm),abs(xtilde)*dij_err2(rgm-2), &
5691 maxzadj2f*d00_err2(rgm-1) ) / abs(zadjkl) ) &
5694 dij_err2(rgm-1)=max(dij_err2(rgm-1), &
5699 if (rgm-n0.le.rmax)
then
5703 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) + dexpgm(n0,n1,n2,n3,gm)
5710 if ((rgm.le.rmax+1))
then
5715 d(0,n1,n2,n3) = d(0,n1,n2,n3) + dexpgm(0,n1,n2,n3,gm)
5718 if(abs(dexpgm(0,n1,n2,n3,gm-1)).ne.0d0)
then
5719 derr(rgm-1)=max(derr(rgm-1),abs(dexpgm(0,n1,n2,n3,gm)) &
5720 *min(1d0,abs(dexpgm(0,n1,n2,n3,gm))/abs(dexpgm(0,n1,n2,n3,gm-1))))
5722 derr(rgm-1)=max(derr(rgm-1),abs(dexpgm(0,n1,n2,n3,gm)))
5734 if(dij_err2(rgm-1).gt.3d0*derr(rgm-1))
then
5736 if(dij_err(rgm-1).gt.3d0*derr(rgm-1))
then
5738 gtrunc = min(gm,gtrunc)
5741 write(*,*)
'CalcDgm exit err',r,rgm-1,gm,gtrunc,dij_err(rgm-1),derr(rgm-1)
5751 write(*,*)
'CalcDgm D(0,0,0,0) = ',r,d(0,0,0,0)
5754 write(*,*)
'CalcDgm D(0,0,1,0) = ',r,d(0,0,1,0)
5756 if(r.gt.2.and.rmax.ge.2)
then
5757 write(*,*)
'CalcDgm D(1,0,0,0) = ',r,d(1,0,0,0)
5758 write(*,*)
'CalcDgm D(0,2,0,0) = ',r,d(0,2,0,0)
5759 write(*,*)
'CalcDgm D(0,0,1,1) = ',r,d(0,0,1,1)
5761 write(*,*)
'CalcDgm D(0,0,2,0) = ',r,d(0,0,2,0)
5763 if(r.gt.3.and.rmax.ge.3)
then
5764 write(*,*)
'CalcDgm D(1,0,1,0) = ',r,d(1,0,1,0)
5765 write(*,*)
'CalcDgm D(1,1,0,0) = ',r,d(1,1,0,0)
5766 write(*,*)
'CalcDgm D(1,0,1,0) = ',r,d(1,0,1,0)
5767 write(*,*)
'CalcDgm D(1,0,0,1) = ',r,d(1,0,0,1)
5769 write(*,*)
'CalcDgm D(0,3,0,0) = ',r,d(0,3,0,0)
5770 write(*,*)
'CalcDgm D(0,2,1,0) = ',r,d(0,2,1,0)
5771 write(*,*)
'CalcDgm D(0,2,0,1) = ',r,d(0,2,0,1)
5772 write(*,*)
'CalcDgm D(0,0,3,0) = ',r,d(0,0,3,0)
5773 write(*,*)
'CalcDgm D(0,1,1,1) = ',r,d(0,1,1,1)
5774 write(*,*)
'CalcDgm D(0,0,2,1) = ',r,d(0,0,2,1)
5776 write(*,*)
'CalcDgm Dij_err',r,dij_err
5777 write(*,*)
'CalcDgm Dij_acc',r,dij_err/abs(d(0,0,0,0))
5779 write(*,*)
'CalcDgm err',r,derr
5780 write(*,*)
'CalcDgm acc',r,derr/abs(d(0,0,0,0))
5783 derr2 = max(derr,dij_err2(0:rmax))
5784 derr = max(derr,dij_err(0:rmax))
5793 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0)
then
5801 d(n0,n1,n2,rgm-2*n0-n1-n2)=0d0
5809 d(0,n1,n2,r-n1-n2)=0d0
5815 111
format(a22,2(
'(',g24.17,
',',g24.17,
') ':))
5816 call seterrflag_coli(-5)
5817 call errout_coli(
'CalcDgm',
' exit rloop for D', &
5819 if (errorwriteflag)
then
5820 write(nerrout_coli,100)
' CalcDgm: exit rloop for D ', &
5821 ' should not appear'
5822 write(nerrout_coli,111)
' CalcDgm: p10 = ',p10
5823 write(nerrout_coli,111)
' CalcDgm: p21 = ',p21
5824 write(nerrout_coli,111)
' CalcDgm: p32 = ',p32
5825 write(nerrout_coli,111)
' CalcDgm: p30 = ',p30
5826 write(nerrout_coli,111)
' CalcDgm: p20 = ',p20
5827 write(nerrout_coli,111)
' CalcDgm: p31 = ',p31
5828 write(nerrout_coli,111)
' CalcDgm: m02 = ',m02
5829 write(nerrout_coli,111)
' CalcDgm: m12 = ',m12
5830 write(nerrout_coli,111)
' CalcDgm: m22 = ',m22
5831 write(nerrout_coli,111)
' CalcDgm: m32 = ',m32
5836 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0.and.r.ge.rmax)
then
5852 d(n0,n1,n2,n3) = (c_0(n0-1,n1,n2,n3) + 2*
mm02*d(n0-1,n1,n2,n3) &
5853 + 4*duv(n0,n1,n2,n3) &
5854 +
f(1)*d(n0-1,n1+1,n2,n3) +
f(2)*d(n0-1,n1,n2+1,n3) &
5855 +
f(3)*d(n0-1,n1,n2,n3+1)) / (2*(r-1))
5871 end subroutine calcdgm
5881 subroutine calcdgr(D,Duv,p10,p21,p32,p30,p20,p31, &
5882 m02,m12,m22,m32,rmax,ordgr_min,ordgr_max,id,Derr,Derr2)
5886 integer,
intent(in) :: rmax,ordgr_min,ordgr_max,id
5887 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
5888 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
5889 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
5890 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
5891 double complex :: Zadjfj,Zadj2(3,3), Zadjkl, Xtilde
5892 double complex,
allocatable :: Dexpgr(:,:,:,:,:), DuvExpgr(:,:,:,:)
5893 double complex,
allocatable :: C_0(:,:,:,:), Cuv_0(:,:,:,:), Shat(:,:,:,:,:)
5894 double complex,
allocatable :: C_i(:,:,:,:), Cuv_i(:,:,:,:)
5895 double complex,
allocatable :: D_alt(:,:,:,:)
5896 double precision,
allocatable :: Cerr_i(:,:),Cerr2_i(:,:)
5897 double complex :: Smod(3), Skl, Daux
5898 double complex :: cC0f, elimminf2_coli
5899 double precision,
allocatable :: D00_err(:),Dij_err(:),Cij_err(:),acc_req_Cextra(:)
5900 double precision,
allocatable :: D00_err2(:),Dij_err2(:),Cij_err2(:)
5902 double precision :: maxDexpgr(0:1,0:2*(rmax+ordgr_min),0:ordgr_max),truncfacexp
5903 integer :: rmaxC,rmaxExp,gtrunc,r,n0,n1,n2,n3,k,l,i,j,m,n,g,rg,nt,mt,nn,nnt,nntt
5904 integer :: inds0(3), inds1(3), inds(3)
5905 integer :: bin,nid(0:3)
5906 logical :: errorwriteflag
5909 write(*,*)
'CalcDgr in, ord',rmax,ordgr_min,ordgr_max
5913 rmaxc = 2*rmax + 2*ordgr_min
5914 allocate(c_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
5915 allocate(cuv_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
5916 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
5917 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
5918 allocate(cerr_i(0:rmaxc,0:3))
5919 allocate(cerr2_i(0:rmaxc,0:3))
5920 allocate(acc_req_cextra(0:rmaxc))
5926 if (mod(id/bin,2).eq.0)
then
5935 acc_req_cextra(0:rmax) = acc_req_cind
5940 acc_req_cextra(r)= acc_req_cextra(r-1)/
fac_gr
5943 acc_req_cextra(rmax+1:rmaxc) = acc_inf
5946 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0),rmax,acc_req_cextra)
5947 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1),rmax,acc_req_cextra)
5948 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2),rmax,acc_req_cextra)
5949 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3),rmax,acc_req_cextra)
5952 write(*,*)
'CalcDgr Cerr 0 =',cerr_i(0:rmaxc,0)
5953 write(*,*)
'CalcDgr Cerr 1 =',cerr_i(0:rmaxc,1)
5954 write(*,*)
'CalcDgr Cerr 2 =',cerr_i(0:rmaxc,2)
5955 write(*,*)
'CalcDgr Cerr 3 =',cerr_i(0:rmaxc,3)
5956 write(*,*)
'CalcDgr Cacc 0 =',cerr_i(0:rmaxc,0)/abs(c_0(0,0,0,0))
5957 write(*,*)
'CalcDgr Cacc 1 =',cerr_i(0:rmaxc,1)/abs(c_i(0,0,0,1))
5958 write(*,*)
'CalcDgr Cacc 2 =',cerr_i(0:rmaxc,2)/abs(c_i(0,0,0,2))
5959 write(*,*)
'CalcDgr Cacc 3 =',cerr_i(0:rmaxc,3)/abs(c_i(0,0,0,3))
5966 n0 = (rmaxc-n1-n2-n3)
5967 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
5968 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
5969 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
5970 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
6024 allocate(shat(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc,3))
6032 shat(n0,n1,n2,n3,:) = -c_0(n0,n1,n2,n3)
6035 shat(n0,n1,n2,n3,1) = shat(n0,n1,n2,n3,1) + c_i(n0,n2,n3,1)
6037 if(n0.eq.3.and.r.eq.6)
then
6044 shat(n0,n1,n2,n3,2) = shat(n0,n1,n2,n3,2) + c_i(n0,n1,n3,2)
6046 if(n0.eq.3.and.r.eq.6)
then
6054 shat(n0,n1,n2,n3,3) = shat(n0,n1,n2,n3,3) + c_i(n0,n1,n2,3)
6056 if(n0.eq.3.and.r.eq.6)
then
6087 zadj2(2,2) = -
z(3,3)
6090 zadj2(3,3) = -
z(2,2)
6101 zadj2(2,3) = -
z(3,1)
6102 zadj2(3,1) = -
z(2,3)
6122 zadj2(2,1) = -
z(3,2)
6125 zadj2(3,2) = -
z(2,1)
6137 zadj2(3,3) = -
z(2,2)
6138 zadj2(2,2) = -
z(3,3)
6149 zadj2(3,1) = -
z(2,3)
6152 zadj2(2,3) = -
z(3,1)
6163 zadj2(3,2) = -
z(2,1)
6164 zadj2(2,1) = -
z(3,2)
6176 zadj2(3,2) = -
z(1,3)
6179 zadj2(1,3) = -
z(3,2)
6190 zadj2(3,3) = -
z(1,1)
6191 zadj2(1,1) = -
z(3,3)
6202 zadj2(3,1) = -
z(1,2)
6205 zadj2(1,2) = -
z(3,1)
6233 allocate(dexpgr(0:rmaxexp/2,0:rmaxexp,0:rmaxexp,0:rmaxexp,0:ordgr_max))
6236 allocate(duvexpgr(0:(rmaxexp+1),0:rmaxexp+1,0:rmaxexp+1,0:rmaxexp+1))
6243 duv(0:rmax,0:rmax,0:rmax,0:rmax) = duvexpgr(0:rmax,0:rmax,0:rmax,0:rmax)
6246 allocate(d00_err(0:rmaxexp))
6247 allocate(dij_err(0:rmaxexp))
6248 allocate(cij_err(0:rmaxc))
6250 allocate(d00_err2(0:rmaxexp))
6251 allocate(dij_err2(0:rmaxexp))
6252 allocate(cij_err2(0:rmaxc))
6258 cij_err = max(cerr_i(:,0),cerr_i(:,1),cerr_i(:,2),cerr_i(:,3))
6263 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,1),cerr2_i(:,2),cerr2_i(:,3))
6266 write(*,*)
'CalcDgr Cij_err = ',cij_err
6267 write(*,*)
'CalcDgr C0_err = ', cerr_i(0,0),cerr_i(0,1),cerr_i(0,2),cerr_i(0,3)
6268 write(*,*)
'CalcDgr C0 = ', c_0(0,0,0,0),c_i(0,0,0,1),c_i(0,0,0,2),c_i(0,0,0,3)
6280 rloop:
do r=0,rmaxexp/2
6286 if (r.gt.rmax+gtrunc)
exit rloop
6289 write(*,*)
'CalcDgr rloop',r,rmaxexp,rmaxc
6299 maxdexpgr(1,r,0)=0d0
6316 write(*,*)
'CalcDgr inds0',n0,inds0
6324 write(*,*)
'CalcDgr inds1',n0,inds1
6328 daux = -
zadj(k,l)*c_0(n0-1,inds1(1),inds1(2),inds1(3))
6342 write(*,*)
'CalcDgr C_0 1c',n0-1,inds1(1),inds1(2),inds1(3),c_0(n0-1,inds1(1),inds1(2),inds1(3))
6343 write(*,*)
'CalcDgr Daux 1c',-
zadj(k,l)*c_0(n0-1,inds1(1),inds1(2),inds1(3))
6344 write(*,*)
'CalcDgr Daux 1s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6350 daux = daux -
zadj(i,l)*shat(n0-1,inds(1),inds(2),inds(3),i)
6352 write(*,*)
'CalcDgr Daux 2ci', -
zadj(i,l)*shat(n0-1,inds(1),inds(2),inds(3),i)
6357 write(*,*)
'CalcDgr Daux 2s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6363 daux = daux +
zadj(k,l)*shat(n0-1,inds(1),inds(2),inds(3),i)
6365 write(*,*)
'CalcDgr Daux 3ci',
zadj(k,l)*shat(n0-1,inds(1),inds(2),inds(3),i)
6371 write(*,*)
'CalcDgr Daux 3s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6374 daux = daux + 2*(nn+1) *zadj2(n ,m )*shat(n0,inds0(1),inds0(2),inds0(3),m) &
6375 + 2*(nn+1) *zadj2(n ,mt)*shat(n0,inds0(1),inds0(2),inds0(3),mt)
6379 write(*,*)
'CalcDgr Daux 4ca', 2*(nn+1) *zadj2(n ,m )*shat(n0,inds0(1),inds0(2),inds0(3),m)
6380 write(*,*)
'CalcDgr Daux 4cb', 2*(nn+1) *zadj2(n ,mt)*shat(n0,inds0(1),inds0(2),inds0(3),mt)
6381 write(*,*)
'CalcDgr Daux 4s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6388 inds(nt) = inds(nt)-1
6389 daux = daux + 2*nnt*zadj2(nt,m )*shat(n0,inds(1),inds(2),inds(3),m) &
6390 + 2*nnt*zadj2(nt,mt)*shat(n0,inds(1),inds(2),inds(3),mt)
6391 daux = daux - 2*nnt*
zadj2f(k,nt,l)*dexpgr(n0,inds(1),inds(2),inds(3),0)
6394 write(*,*)
'CalcDgr Daux 5ci', 2*nnt*zadj2(nt,m )*shat(n0,inds(1),inds(2),inds(3),m)
6395 write(*,*)
'CalcDgr Daux 5ci', 2*nnt*zadj2(nt,mt)*shat(n0,inds(1),inds(2),inds(3),mt)
6396 write(*,*)
'CalcDgr Daux 5ci', 2*nnt*
zadj2f(k,nt,l)*dexpgr(n0,inds(1),inds(2),inds(3),0)
6397 write(*,*)
'CalcDgr Daux 5cii', 2*nnt*zadj2(nt,m ),shat(n0,inds(1),inds(2),inds(3),m)
6398 write(*,*)
'CalcDgr Daux 5cii', 2*nnt*zadj2(nt,mt),shat(n0,inds(1),inds(2),inds(3),mt)
6399 write(*,*)
'CalcDgr Daux 5cii', 2*nnt*
zadj2f(k,nt,l),dexpgr(n0,inds(1),inds(2),inds(3),0)
6400 write(*,*)
'CalcDgr Daux 5s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6407 if (inds(n).gt.1)
then
6409 daux = daux - 4*(nn+1)*nn * zadj2(n,m ) * dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6411 write(*,*)
'CalcDgr Daux 6c',4*(nn+1)*nn* zadj2(n,m ) *dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6412 write(*,*)
'CalcDgr Daux 6s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6416 if (inds(n).gt.0.and.inds(m).gt.0)
then
6419 daux = daux - 4*(nn+1)*(inds(m)+1)* zadj2(n,m ) * dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6421 write(*,*)
'CalcDgr Daux 6c',-4*(nn+1)*(inds(m)+1)* zadj2(n,m ) *dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6422 write(*,*)
'CalcDgr Daux 6s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6430 if (inds(nt).gt.1)
then
6431 inds(nt) = inds(nt)-2
6432 daux = daux - 4*nnt*(nnt-1) * zadj2(nt,m ) * dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6434 write(*,*)
'CalcDgr Daux 7c',4*nnt*(nnt-1) * zadj2(nt,m ) *dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6435 write(*,*)
'CalcDgr Daux 7s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6439 if (inds(nt).gt.0.and.inds(m).gt.0)
then
6440 inds(nt) = inds(nt)-1
6442 daux = daux - 4*nnt*(inds(m)+1)* zadj2(nt,m )* dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6444 write(*,*)
'CalcDgr Daux 7c',4*nnt*(inds(m)+1)* zadj2(nt,m )* dexpgr(n0,inds(1),inds(2),inds(3),0)
6445 write(*,*)
'CalcDgr Daux 7s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6453 if (inds(n).gt.1)
then
6455 daux = daux - 4*(nn+1)*nn * zadj2(n ,mt)* dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6457 write(*,*)
'CalcDgr Daux 8c',- 4*(nn+1)*nn * zadj2(n ,mt)* dexpgr(n0+1,inds(1),inds(2),inds(3),0) &
6458 , n0+1,inds(1),inds(2),inds(3)
6459 write(*,*)
'CalcDgr Daux 8s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6463 if (inds(n).gt.0.and.inds(mt).gt.0)
then
6465 inds(mt) = inds(mt)-1
6466 daux = daux - 4*(nn+1)*(inds(mt)+1)* zadj2(n ,mt)* dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6468 write(*,*)
'CalcDgr Daux 8c',- 4*(nn+1)*(inds(mt)+1)* zadj2(n ,mt)* dexpgr(n0+1,inds(1),inds(2),inds(3),0) &
6469 ,n0+1,inds(1),inds(2),inds(3)
6470 write(*,*)
'CalcDgr Daux 8s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6478 if (inds(nt).gt.1)
then
6479 inds(nt) = inds(nt)-2
6480 daux = daux - 4*nnt*(nnt-1) * zadj2(nt,mt)* dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6482 write(*,*)
'CalcDgr Daux 9c', - 4*nnt*(nnt-1) * zadj2(nt,mt) * dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6483 write(*,*)
'CalcDgr Daux 9s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6487 if (inds(nt).gt.0.and.inds(mt).gt.0)
then
6488 inds(nt) = inds(nt)-1
6489 inds(mt) = inds(mt)-1
6490 daux = daux - 4*nnt*(inds(mt)+1) * zadj2(nt,mt)* dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6492 write(*,*)
'CalcDgr Daux 9c',- 4*nnt*(inds(mt)+1) * zadj2(nt,mt) * dexpgr(n0+1,inds(1),inds(2),inds(3),0)
6493 write(*,*)
'CalcDgr Daux 9s',daux,daux/(2*(nn+1)*
zadj2f(k,n,l))
6498 dexpgr(n0,inds0(1),inds0(2),inds0(3),0) = daux/(2*(nn+1)*
zadj2f(k,n,l))
6501 write(*,*)
'CalcDgr Dexpgr',n0,inds0(1),inds0(2),inds0(3),dexpgr(n0,inds0(1),inds0(2),inds0(3),0)
6505 maxdexpgr(1,r,0) = maxdexpgr(1,r,0) + abs(dexpgr(n0,inds0(1),inds0(2),inds0(3),0) )
6510 d(n0,inds0(1),inds0(2),inds0(3)) = dexpgr(n0,inds0(1),inds0(2),inds0(3),0)
6520 maxdexpgr(0,r,0)=0d0
6525 smod = shat(0,n1,n2,n3,:)
6527 smod(1) = smod(1) - 2d0*n1*dexpgr(1,n1-1,n2,n3,0)
6530 smod(2) = smod(2) - 2d0*n2*dexpgr(1,n1,n2-1,n3,0)
6533 smod(3) = smod(3) - 2d0*n3*dexpgr(1,n1,n2,n3-1,0)
6536 dexpgr(0,n1,n2,n3,0) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
6537 +
zadj(3,j)*smod(3))/zadjfj
6538 maxdexpgr(0,r,0) = maxdexpgr(0,r,0) + abs(dexpgr(0,n1,n2,n3,0))
6540 d(0,n1,n2,n3) = dexpgr(0,n1,n2,n3,0)
6550 if(n0.eq.0.and.n1.eq.3.and.n2.eq.0.and.n3.eq.0)
then
6551 write(*,*)
'Smod(0,3,0,0,1)= ',shat(0,n1,n2,n3,1)
6552 write(*,*)
'Smod(0,3,0,0,2)= ',shat(0,n1,n2,n3,2)
6553 write(*,*)
'Smod(0,3,0,0,3)= ',shat(0,n1,n2,n3,3)
6554 write(*,*)
'D(1,2,0,0)= ',0,dexpgr(1,2,n2,n3,0)
6555 write(*,*)
'D(0,3,0,0)= ',0,d(n0,n1,n2,n3)
6568 derr(r) =
fac_gr*maxdexpgr(0,r,0)
6575 dij_err(r)=
maxzadj*max(cij_err(r),2*d00_err(r+1))/abs(zadjfj)
6580 dij_err2(r)=
maxzadj*max(cij_err2(r),2*d00_err2(r+1))/abs(zadjfj)
6587 gloop:
do g=1,min(gtrunc,r)
6591 write(*,*)
'CalcDgr: gloop ',r,rg,g
6597 maxdexpgr(1,rg,g) = 0d0
6600 do nnt=rg-n0-nn,0,-1
6614 daux = 2*
zadj(k,l) * (2+rg-n0) * dexpgr(n0,inds1(1),inds1(2),inds1(3),g-1)
6625 inds(k) = inds(k) + 1
6626 inds(l) = inds(l) + 1
6627 daux = daux +
detz * dexpgr(n0-1,inds(1),inds(2),inds(3),g-2)
6636 inds(k) = inds(k) + 1
6637 daux = daux +
zadjf(l) * dexpgr(n0-1,inds(1),inds(2),inds(3),g-1)
6649 inds(nt) = inds(nt)-1
6650 daux = daux - 2*nnt*
zadj2f(k,nt,l)*dexpgr(n0,inds(1),inds(2),inds(3),g)
6661 if (inds(n).gt.1)
then
6663 daux = daux - 4*(nn+1)*nn * zadj2(n,m ) * dexpgr(n0+1,inds(1),inds(2),inds(3),g)
6670 if (inds(n).gt.0.and.inds(m).gt.0)
then
6673 daux = daux - 4*(nn+1)*(inds(m)+1)* zadj2(n,m ) * dexpgr(n0+1,inds(1),inds(2),inds(3),g)
6684 if (inds(nt).gt.1)
then
6685 inds(nt) = inds(nt)-2
6686 daux = daux - 4*nnt*(nnt-1) * zadj2(nt,m ) * dexpgr(n0+1,inds(1),inds(2),inds(3),g)
6693 if (inds(nt).gt.0.and.inds(m).gt.0)
then
6694 inds(nt) = inds(nt)-1
6696 daux = daux - 4*nnt*(inds(m)+1)* zadj2(nt,m )* dexpgr(n0+1,inds(1),inds(2),inds(3),g)
6706 if (inds(n).gt.1)
then
6708 daux = daux - 4*(nn+1)*nn * zadj2(n ,mt)* dexpgr(n0+1,inds(1),inds(2),inds(3),g)
6716 if (inds(n).gt.0.and.inds(mt).gt.0)
then
6718 inds(mt) = inds(mt)-1
6719 daux = daux - 4*(nn+1)*(inds(mt)+1)* zadj2(n ,mt)* dexpgr(n0+1,inds(1),inds(2),inds(3),g)
6729 if (inds(nt).gt.1)
then
6730 inds(nt) = inds(nt)-2
6731 daux = daux - 4*nnt*(nnt-1) * zadj2(nt,mt)* dexpgr(n0+1,inds(1),inds(2),inds(3),g)
6738 if (inds(nt).gt.0.and.inds(mt).gt.0)
then
6739 inds(nt) = inds(nt)-1
6740 inds(mt) = inds(mt)-1
6741 daux = daux - 4*nnt*(inds(mt)+1) * zadj2(nt,mt)* dexpgr(n0+1,inds(1),inds(2),inds(3),g)
6751 dexpgr(n0,inds0(1),inds0(2),inds0(3),g) = daux/(2*(nn+1)*
zadj2f(k,n,l))
6754 maxdexpgr(1,rg,g) = maxdexpgr(1,rg,g) + abs(dexpgr(n0,inds0(1),inds0(2),inds0(3),g))
6757 if (g.eq.1.and.abs(dexpgr(1,inds0(1),inds0(2),inds0(3),g)).gt. &
6758 truncfacexp*max(1/
m2scale,maxdexpgr(1,rg,g-1)) .or. &
6759 g.ge.2.and.abs(dexpgr(1,inds0(1),inds0(2),inds0(3),g)).gt. &
6760 truncfacexp*maxdexpgr(1,rg,g-1))
then
6778 if (rg.le.rmax)
then
6781 if (rg.le.rmax)
then
6785 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) + dexpgr(n0,n1,n2,n3,g)
6800 maxdexpgr(0,rg,g) = 0d0
6807 smod(1) = smod(1) - 2d0*n1*dexpgr(1,n1-1,n2,n3,g)
6810 smod(2) = smod(2) - 2d0*n2*dexpgr(1,n1,n2-1,n3,g)
6813 smod(3) = smod(3) - 2d0*n3*dexpgr(1,n1,n2,n3-1,g)
6820 dexpgr(0,n1,n2,n3,g) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
6821 +
zadj(3,j)*smod(3) &
6822 -
detz*dexpgr(0,inds(1),inds(2),inds(3),g-1))/zadjfj
6824 maxdexpgr(0,rg,g) = maxdexpgr(0,rg,g) + abs(dexpgr(0,n1,n2,n3,g))
6835 if (g.eq.1.and.abs(dexpgr(0,n1,n2,n3,g)).gt. &
6836 truncfacexp*max(1/
m2scale,maxdexpgr(0,rg,g-1)) .or. &
6837 g.ge.2.and.abs(dexpgr(0,n1,n2,n3,g)).gt. &
6838 truncfacexp*maxdexpgr(0,rg,g-1))
then
6841 write(*,*)
'CalcDgr exit gloop',0,n1,n2,n3,g,abs(dexpgr(0,n1,n2,n3,g)),maxdexpgr(0,rg,g-1),truncfacexp
6852 d00_err(rg+1) = max( d00_err(rg+1), &
6853 max(
maxzadj*(2+rg-2*n0)*d00_err(rg+2), &
6854 abs(
detz)*dij_err(rg+2), &
6858 dij_err(rg)=max(dij_err(rg), &
6859 max(2*
maxzadj*d00_err(rg+1),abs(
detz)*dij_err(rg))/abs(zadjfj) )
6862 d00_err2(rg+1) = max( d00_err2(rg+1), &
6863 max(
maxzadj*(2+rg-2*n0)*d00_err2(rg+2), &
6864 abs(
detz)*dij_err2(rg+2), &
6868 dij_err2(rg)=max(dij_err2(rg), &
6869 max(2*
maxzadj*d00_err2(rg+1),abs(
detz)*dij_err2(rg))/abs(zadjfj) )
6872 if (rg.le.rmax)
then
6875 if (rg.le.rmax)
then
6879 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) + dexpgr(n0,n1,n2,n3,g)
6887 if (rg.le.rmax)
then
6892 d(0,n1,n2,n3) = d(0,n1,n2,n3) + dexpgr(0,n1,n2,n3,g)
6893 if(abs(dexpgr(0,n1,n2,n3,g-1)).ne.0d0)
then
6895 derr(rg)=max(derr(rg),abs(dexpgr(0,n1,n2,n3,g))*min(1d0,abs(dexpgr(0,n1,n2,n3,g))/abs(dexpgr(0,n1,n2,n3,g-1))))
6897 derr(rg)=max(derr(rg),abs(dexpgr(0,n1,n2,n3,g)))
6909 if(dij_err2(rg).gt.3d0*derr(rg))
then
6911 if(dij_err(rg).gt.3d0*derr(rg))
then
6913 gtrunc = min(g,gtrunc)
6916 write(*,*)
'CalcDgr exit err',r,rg,g,gtrunc,dij_err(rg),derr(rg)
6926 write(*,*)
'CalcDgr D(0,0,0,0) = ',r,d(0,0,0,0)
6928 write(*,*)
'CalcDgr D(1,0,0,0) = ',r,d(1,0,0,0)
6929 write(*,*)
'CalcDgr D(0,1,0,0) = ',r,d(0,1,0,0)
6930 write(*,*)
'CalcDgr D(0,0,1,0) = ',r,d(0,0,1,0)
6931 write(*,*)
'CalcDgr D(0,0,0,1) = ',r,d(0,0,0,1)
6933 if(r.ge.2.and.rmax.ge.2)
then
6934 write(*,*)
'CalcDgr D(1,1,0,0) = ',r,d(1,1,0,0)
6935 write(*,*)
'CalcDgr D(1,0,1,0) = ',r,d(1,0,1,0)
6936 write(*,*)
'CalcDgr D(1,0,0,1) = ',r,d(1,0,0,1)
6937 write(*,*)
'CalcDgr D(0,2,0,0) = ',r,d(0,2,0,0)
6939 write(*,*)
'CalcDgr D(0,0,2,0) = ',r,d(0,0,2,0)
6941 if(r.ge.3.and.rmax.ge.2)
then
6944 write(*,*)
'CalcDgr D(1,2,0,0) = ',r,d(1,2,0,0)
6945 write(*,*)
'CalcDgr D(1,0,2,0) = ',r,d(1,0,2,0)
6946 write(*,*)
'CalcDgr D(0,3,0,0) = ',r,d(0,3,0,0)
6947 write(*,*)
'CalcDgr D(0,2,1,0) = ',r,d(0,2,1,0)
6948 write(*,*)
'CalcDgr D(0,0,3,0) = ',r,d(0,0,3,0)
6949 write(*,*)
'CalcDgr D(0,1,1,1) = ',r,d(0,1,1,1)
6950 write(*,*)
'CalcDgr D(0,0,2,1) = ',r,d(0,0,2,1)
6952 write(*,*)
'CalcDgr Dij_err',r,dij_err
6953 write(*,*)
'CalcDgr Dij_acc',r,dij_err/abs(d(0,0,0,0))
6955 write(*,*)
'CalcDgr err',r,derr
6956 write(*,*)
'CalcDgr acc',r,derr/abs(d(0,0,0,0))
6959 derr2 = max(derr,dij_err2(0:rmax))
6960 derr = max(derr,dij_err(0:rmax))
6969 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0)
then
6975 d(n0,n1,n2,rg-n0-n1-n2)=0d0
6982 111
format(a22,2(
'(',g24.17,
',',g24.17,
') ':))
6983 call seterrflag_coli(-5)
6984 call errout_coli(
'CalcDgr',
' exit rloop for D', &
6986 if (errorwriteflag)
then
6987 write(nerrout_coli,100)
' CalcDgr: exit rloop for D ', &
6988 ' should not appear'
6989 write(nerrout_coli,111)
' CalcDgr: p10 = ',p10
6990 write(nerrout_coli,111)
' CalcDgr: p21 = ',p21
6991 write(nerrout_coli,111)
' CalcDgr: p32 = ',p32
6992 write(nerrout_coli,111)
' CalcDgr: p30 = ',p30
6993 write(nerrout_coli,111)
' CalcDgr: p20 = ',p20
6994 write(nerrout_coli,111)
' CalcDgr: p31 = ',p31
6995 write(nerrout_coli,111)
' CalcDgr: m02 = ',m02
6996 write(nerrout_coli,111)
' CalcDgr: m12 = ',m12
6997 write(nerrout_coli,111)
' CalcDgr: m22 = ',m22
6998 write(nerrout_coli,111)
' CalcDgr: m32 = ',m32
7003 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0.and.r.ge.rmax)
then
7017 write(*,*)
'CalcDgr final err',derr
7018 write(*,*)
'CalcDgr final acc',derr/abs(d(0,0,0,0))
7035 subroutine calcdgx(D,Duv,p10,p21,p32,p30,p20,p31, &
7036 m02,m12,m22,m32,rmax,ordgx_min,ordgx_max,id,Derr,Derr2)
7040 integer,
intent(in) :: rmax,ordgx_min,ordgx_max,id
7041 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
7042 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
7043 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
7044 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
7045 double complex :: Zadjfj,Zadj2(4), Zadjkl, Xtilde
7046 double complex,
allocatable :: Dexpgx(:,:,:,:,:), DuvExpgx(:,:,:,:)
7047 double complex,
allocatable :: C_0(:,:,:,:), Cuv_0(:,:,:,:), Shat(:,:,:,:,:)
7048 double complex,
allocatable :: C_i(:,:,:,:), Cuv_i(:,:,:,:)
7049 double complex,
allocatable :: D_alt(:,:,:,:)
7050 double precision,
allocatable :: Cerr_i(:,:),Cerr2_i(:,:)
7051 double complex :: Smod(3), Skl, Daux, DexpgAux
7052 double complex :: cC0f, elimminf2_coli
7053 double precision,
allocatable :: D00_err(:),Dij_err(:),Cij_err(:),acc_req_Cextra(:)
7054 double precision,
allocatable :: D00_err2(:),Dij_err2(:),Cij_err2(:)
7055 double precision :: maxDexpgx(0:1,0:rmax+ordgx_min+1,0:ordgx_max),truncfacexp
7056 integer :: rmaxC,rmaxExp,gtrunc,r,n0,n1,n2,n3,k,l,i,j,m,n,g,rg,lt,ltt,nl,nlt,nltt
7057 integer :: inds0(3), inds(3), inds2(2,4)
7058 integer :: bin,nid(0:3)
7059 logical :: errorwriteflag
7062 write(*,*)
'CalcDgx in, ord',rmax,ordgx_min,ordgx_max
7066 rmaxc = rmax + ordgx_min
7067 allocate(c_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
7068 allocate(cuv_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
7069 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
7070 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
7071 allocate(cerr_i(0:rmaxc,0:3))
7072 allocate(cerr2_i(0:rmaxc,0:3))
7073 allocate(acc_req_cextra(0:rmaxc))
7079 if (mod(id/bin,2).eq.0)
then
7088 acc_req_cextra(0:rmax) = acc_req_cind
7089 if (
x_g.ne.0d0)
then
7091 acc_req_cextra(r)= acc_req_cextra(r-1)/
x_g
7094 acc_req_cextra(rmax+1:rmaxc) = acc_inf
7097 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0),rmax,acc_req_cextra)
7098 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1),rmax,acc_req_cextra)
7099 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2),rmax,acc_req_cextra)
7100 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3),rmax,acc_req_cextra)
7106 n0 = (rmaxc-n1-n2-n3)
7107 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
7108 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
7109 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
7110 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
7164 allocate(shat(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc,3))
7172 shat(n0,n1,n2,n3,:) = -c_0(n0,n1,n2,n3)
7175 shat(n0,n1,n2,n3,1) = shat(n0,n1,n2,n3,1) + c_i(n0,n2,n3,1)
7179 shat(n0,n1,n2,n3,2) = shat(n0,n1,n2,n3,2) + c_i(n0,n1,n3,2)
7183 shat(n0,n1,n2,n3,3) = shat(n0,n1,n2,n3,3) + c_i(n0,n1,n2,3)
7276 write(*,*)
'CalcDgx i,j,l',i,j,l,lt,ltt
7278 write(*,*)
'CalcDgx pars', abs(
zadjf(j)),abs(
xadj(i,j))
7291 allocate(dexpgx(0:rmaxexp/2,0:rmaxexp,0:rmaxexp,0:rmaxexp,0:ordgx_max))
7295 allocate(duvexpgx(0:rmaxexp,0:rmaxexp,0:rmaxexp,0:rmaxexp))
7297 duv(0:rmax,0:rmax,0:rmax,0:rmax) = duvexpgx(0:rmax,0:rmax,0:rmax,0:rmax)
7300 allocate(d00_err(0:rmaxexp))
7301 allocate(dij_err(0:rmaxexp))
7302 allocate(cij_err(0:rmaxc))
7304 allocate(d00_err2(0:rmaxexp))
7305 allocate(dij_err2(0:rmaxexp))
7306 allocate(cij_err2(0:rmaxc))
7312 cij_err = max(cerr_i(:,0),cerr_i(:,1),cerr_i(:,2),cerr_i(:,3))
7317 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,1),cerr2_i(:,2),cerr2_i(:,3))
7320 write(*,*)
'CalcDgx Cij_err = ',cij_err
7321 write(*,*)
'CalcDgx C0_err = ', cerr_i(0,0),cerr_i(0,1),cerr_i(0,2),cerr_i(0,3)
7322 write(*,*)
'CalcDgx C0 = ', c_0(0,0,0,0),c_i(0,0,0,1),c_i(0,0,0,2),c_i(0,0,0,3)
7335 rloop:
do r=1,rmaxexp
7341 if (r.gt.rmax+gtrunc+1)
exit rloop
7344 write(*,*)
'CalcDgx rloop',r
7352 maxdexpgx(1,r,0)=0d0
7363 daux =
zadj2f(i,j,1)*shat(0,inds(1),inds(2),inds(3),1) &
7364 +
zadj2f(i,j,2)*shat(0,inds(1),inds(2),inds(3),2) &
7365 +
zadj2f(i,j,3)*shat(0,inds(1),inds(2),inds(3),3)
7369 daux = daux -
zadj(i,j)*(c_0(0,inds(1),inds(2),inds(3)) &
7370 +4*duvexpgx(1,inds(1),inds(2),inds(3)))
7374 daux = daux - 2*nlt*
zadj2f(i,j,lt)*dexpgx(1,inds(1),inds(2),inds(3),0)
7379 daux = daux - 2*nltt*
zadj2f(i,j,ltt)*dexpgx(1,inds(1),inds(2),inds(3),0)
7382 dexpgx(1,inds0(1),inds0(2),inds0(3),0) = daux/(2*(nl+1)*
zadj2f(i,j,l))
7384 maxdexpgx(1,r,0) = maxdexpgx(1,r,0) + abs(dexpgx(1,inds0(1),inds0(2),inds0(3),0) )
7387 d(1,inds0(1),inds0(2),inds0(3)) = dexpgx(1,inds0(1),inds0(2),inds0(3),0)
7396 maxdexpgx(0,r-1,0)=0d0
7401 smod = shat(0,n1,n2,n3,:)
7403 smod(1) = smod(1) - 2d0*n1*dexpgx(1,n1-1,n2,n3,0)
7406 smod(2) = smod(2) - 2d0*n2*dexpgx(1,n1,n2-1,n3,0)
7409 smod(3) = smod(3) - 2d0*n3*dexpgx(1,n1,n2,n3-1,0)
7412 dexpgx(0,n1,n2,n3,0) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
7413 +
zadj(3,j)*smod(3))/zadjfj
7414 maxdexpgx(0,r-1,0) = maxdexpgx(0,r-1,0) + abs(dexpgx(0,n1,n2,n3,0))
7415 if (r.le.rmax+1)
then
7416 d(0,n1,n2,n3) = dexpgx(0,n1,n2,n3,0)
7433 if(r-1.le.rmax)
then
7435 derr(r-1) =
fac_g*maxdexpgx(0,r-1,0)
7442 dij_err(r-1)=
maxzadj*max(cij_err(r-1),2*d00_err(r))/abs(zadjfj)
7447 dij_err2(r-1)=
maxzadj*max(cij_err2(r-1),2*d00_err2(r))/abs(zadjfj)
7454 gloop:
do g=1,min(gtrunc,r-1)
7457 write(*,*)
'gloop ',g,rg
7460 maxdexpgx(1,rg,g) = 0d0
7470 daux = -
xadj(i,j)*dexpgx(0,inds(1),inds(2),inds(3),g-1) &
7471 +
zadj(i,j)*2*rg*dexpgx(1,inds(1),inds(2),inds(3),g-1)
7473 write(*,*)
'CalcDgx con Xij',-
xadj(i,j)*dexpgx(0,inds(1),inds(2),inds(3),g-1)/(2*(nl+1)*
zadj2f(i,j,l))
7474 write(*,*)
'CalcDgx con Zij',+
zadj(i,j)*2*(1+rg)*dexpgx(1,inds(1),inds(2),inds(3),g-1)/(2*(nl+1)*
zadj2f(i,j,l))
7477 daux = daux - zadjfj*dexpgx(0,inds(1),inds(2),inds(3),g-1)
7478 write(*,*)
'CalcDgx con Zadj2f', - zadjfj*dexpgx(0,inds(1),inds(2),inds(3),g-1)/(2*(nl+1)*
zadj2f(i,j,l))
7484 daux = daux - 2*nlt*
zadj2f(i,j,lt)*dexpgx(1,inds(1),inds(2),inds(3),g)
7490 daux = daux - 2*nltt*
zadj2f(i,j,ltt)*dexpgx(1,inds(1),inds(2),inds(3),g)
7493 dexpgx(1,inds0(1),inds0(2),inds0(3),g) = daux/(2*(nl+1)*
zadj2f(i,j,l))
7495 maxdexpgx(1,rg,g) = maxdexpgx(1,rg,g) + abs(dexpgx(1,inds0(1),inds0(2),inds0(3),g) )
7497 write(*,*)
'CalcDgx gloop 00',g,rg,nl,nlt,nltt,dexpgx(1,inds0(1),inds0(2),inds0(3),g)
7500 if (g.eq.1.and.abs(dexpgx(1,inds0(1),inds0(2),inds0(3),g)).gt. &
7501 truncfacexp*max(1/
m2scale,maxdexpgx(1,rg,g-1)) .or. &
7502 g.ge.2.and.abs(dexpgx(1,inds0(1),inds0(2),inds0(3),g)).gt. &
7503 truncfacexp*maxdexpgx(1,rg,g-1))
then
7506 write(*,*)
'CalcDgx cycle loop',1,inds0(1),inds0(2),inds0(3),g, &
7507 abs(dexpgx(1,inds0(1),inds0(2),inds0(3),g)),abs(dexpgx(1,inds0(1),inds0(2),inds0(3),g-1)),maxdexpgx(1,rg,g-1)
7518 if (rg.le.rmax)
then
7522 d(1,n1,n2,n3) = d(1,n1,n2,n3) + dexpgx(1,n1,n2,n3,g)
7535 maxdexpgx(0,rg-1,g) = 0d0
7542 smod(1) = smod(1) - 2d0*n1*dexpgx(1,n1-1,n2,n3,g)
7545 smod(2) = smod(2) - 2d0*n2*dexpgx(1,n1,n2-1,n3,g)
7548 smod(3) = smod(3) - 2d0*n3*dexpgx(1,n1,n2,n3-1,g)
7555 dexpgx(0,n1,n2,n3,g) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
7556 +
zadj(3,j)*smod(3) &
7557 -
detz*dexpgx(0,inds(1),inds(2),inds(3),g-1))/zadjfj
7559 maxdexpgx(0,rg-1,g) = maxdexpgx(0,rg-1,g) + abs(dexpgx(0,n1,n2,n3,g))
7570 if (g.eq.1.and.abs(dexpgx(0,n1,n2,n3,g)).gt. &
7571 truncfacexp*max(1/
m2scale**2,maxdexpgx(0,rg,g-1)) .or. &
7572 g.ge.2.and.abs(dexpgx(0,n1,n2,n3,g)).gt. &
7573 truncfacexp*maxdexpgx(0,rg,g-1))
then
7576 write(*,*)
'CalcDgx exit gloop',0,n1,n2,n3,g,abs(dexpgx(0,n1,n2,n3,g)),maxdexpgx(0,rg-1,g-1),truncfacexp
7587 d00_err(rg) = max( d00_err(rg), &
7588 max( abs(m02)*dij_err(rg-2), &
7589 max(
maxzadjf*dij_err(rg),abs(xtilde)*dij_err(rg-1), &
7593 dij_err(rg-1)=max(dij_err(rg-1), &
7594 max(2*
maxzadj*d00_err(rg),abs(
detz)*dij_err(rg))/abs(zadjfj) )
7597 d00_err2(rg) = max( d00_err2(rg), &
7598 max( abs(m02)*dij_err2(rg-2), &
7599 max(
maxzadjf*dij_err2(rg),abs(xtilde)*dij_err2(rg-1), &
7603 dij_err2(rg-1)=max(dij_err2(rg-1), &
7604 max(2*
maxzadj*d00_err2(rg),abs(
detz)*dij_err2(rg))/abs(zadjfj) )
7607 if (rg.le.rmax)
then
7611 d(0,n1,n2,n3) = d(0,n1,n2,n3) + dexpgx(0,n1,n2,n3,g)
7617 if ((rg.le.rmax+1))
then
7622 d(0,n1,n2,n3) = d(0,n1,n2,n3) + dexpgx(0,n1,n2,n3,g)
7623 if(abs(dexpgx(0,n1,n2,n3,g-1)).ne.0d0)
then
7625 derr(rg-1)=max(derr(rg-1),abs(dexpgx(0,n1,n2,n3,g))*min(1d0,abs(dexpgx(0,n1,n2,n3,g))/abs(dexpgx(0,n1,n2,n3,g-1))))
7627 derr(rg-1)=max(derr(rg-1),abs(dexpgx(0,n1,n2,n3,g)))
7639 if(dij_err2(rg-1).gt.3d0*derr(rg-1))
then
7641 if(dij_err(rg-1).gt.3d0*derr(rg-1))
then
7643 gtrunc = min(g,gtrunc)
7646 write(*,*)
'CalcDgx exit err',r,rg-1,g,gtrunc,dij_err(rg-1),derr(rg-1)
7656 write(*,*)
'CalcDgx D(0,0,0,0) = ',r,d(0,0,0,0)
7658 write(*,*)
'CalcDgx D(1,0,0,0) = ',r,d(1,0,0,0)
7659 write(*,*)
'CalcDgx D(0,1,0,0) = ',r,d(0,1,0,0)
7660 write(*,*)
'CalcDgx D(0,0,1,0) = ',r,d(0,0,1,0)
7662 if(r.gt.2.and.rmax.ge.2)
then
7663 write(*,*)
'CalcDgx D(1,1,0,0) = ',r,d(1,1,0,0)
7666 write(*,*)
'CalcDgx D(0,0,2,0) = ',r,d(0,0,2,0)
7668 if(r.gt.3.and.rmax.ge.2)
then
7669 write(*,*)
'CalcDgx D(1,0,1,0) = ',r,d(1,0,1,0)
7670 write(*,*)
'CalcDgx D(1,1,0,0) = ',r,d(1,1,0,0)
7672 write(*,*)
'CalcDgx D(0,3,0,0) = ',r,d(0,3,0,0)
7673 write(*,*)
'CalcDgx D(0,2,1,0) = ',r,d(0,2,1,0)
7674 write(*,*)
'CalcDgx D(0,0,3,0) = ',r,d(0,0,3,0)
7675 write(*,*)
'CalcDgx D(0,1,1,1) = ',r,d(0,1,1,1)
7676 write(*,*)
'CalcDgx D(0,0,2,1) = ',r,d(0,0,2,1)
7678 write(*,*)
'CalcDgx Dij_err',r,dij_err
7679 write(*,*)
'CalcDgx Dij_acc',r,dij_err/abs(d(0,0,0,0))
7681 write(*,*)
'CalcDgx err',r,derr
7682 write(*,*)
'CalcDgx acc',r,derr/abs(d(0,0,0,0))
7685 derr2 = max(derr,dij_err2(0:rmax))
7686 derr = max(derr,dij_err(0:rmax))
7696 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0)
then
7698 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0.and.r.ge.rmax)
then
7704 d(n0,n1,n2,rg-2*n0-n1-n2)=0d0
7723 write(*,*)
'CalcDgx final err',derr
7724 write(*,*)
'CalcDgx final acc',derr/abs(d(0,0,0,0))
7747 subroutine calcdgy(D,Duv,p10,p21,p32,p30,p20,p31, &
7748 m02,m12,m22,m32,rmax,ordgy_min,ordgy_max,id,Derr,Derr2)
7752 integer,
intent(in) :: rmax,ordgy_min,ordgy_max,id
7753 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
7754 double complex ::Zadj2(4)
7755 double complex,
allocatable :: Dexpgy(:,:,:,:,:), DuvExpgy(:,:,:,:)
7756 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
7757 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
7758 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
7759 double complex,
allocatable :: C_0(:,:,:,:), C_i(:,:,:,:), Shat(:,:,:,:,:)
7760 double complex,
allocatable :: Cuv_0(:,:,:,:), Cuv_i(:,:,:,:)
7761 double complex,
allocatable :: D_alt(:,:,:,:)
7762 double precision,
allocatable :: Cerr_i(:,:),Cerr2_i(:,:)
7763 double complex :: Smod(3), Daux, elimminf2_coli
7764 double precision,
allocatable :: D00_err(:),Dij_err(:),Cij_err(:),acc_req_Cextra(:)
7765 double precision,
allocatable :: D00_err2(:),Dij_err2(:),Cij_err2(:)
7766 double precision :: maxDexpgy(0:1,0:rmax+2*ordgy_min,0:ordgy_max),truncfacexp,acc_aux
7767 integer :: rmaxC,rmaxExp,gtrunc,r,n0,n1,n2,n3,a,b,i,g,rg,m,n
7768 integer :: inds0(3),inds(3),inds2(2,4),at,bt,k,l,lt,ltt,nl,nlt,nltt
7769 integer :: bin,nid(0:3)
7770 logical :: errorwriteflag
7773 write(*,*)
'CalcDgy in, ord',rmax,ordgy_min,ordgy_max
7777 rmaxc = rmax + 2*ordgy_min + 1
7778 allocate(c_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
7779 allocate(cuv_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
7780 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
7781 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
7782 allocate(cerr_i(0:rmaxc,0:3))
7783 allocate(cerr2_i(0:rmaxc,0:3))
7784 allocate(acc_req_cextra(0:rmaxc))
7790 if (mod(id/bin,2).eq.0)
then
7799 acc_req_cextra(0:rmax+1) = acc_req_cind
7801 if (
y_gy.ne.0d0)
then
7803 acc_req_cextra(rmax+2*g) = acc_req_cextra(rmax+2*g-2)/
y_gy
7804 acc_req_cextra(rmax+2*g+1) = acc_req_cextra(rmax+2*g-1)/
y_gy
7806 acc_req_cextra(rmax+g+1) = min(acc_req_cextra(rmax+g+1),acc_aux)
7808 else if(
x_gy.ne.0d0)
then
7810 acc_aux = acc_aux/
x_gy
7811 acc_req_cextra(rmax+g+1) = acc_aux
7814 acc_req_cextra(rmax+2:rmax+2*ordgy_min+1) = acc_inf
7820 write(*,*)
'CalcDgy: accreq_Cextra',acc_req_cextra
7823 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0),rmax,acc_req_cextra)
7824 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1),rmax,acc_req_cextra)
7825 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2),rmax,acc_req_cextra)
7826 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3),rmax,acc_req_cextra)
7829 write(*,*)
'CalcDgy Cerr 0',cerr_i(:,0)
7830 write(*,*)
'CalcDgy Cerr 1',cerr_i(:,1)
7831 write(*,*)
'CalcDgy Cerr 2',cerr_i(:,2)
7832 write(*,*)
'CalcDgy Cerr 3',cerr_i(:,3)
7840 n0 = (rmaxc-n1-n2-n3)
7841 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
7842 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
7843 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
7844 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
7914 allocate(shat(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc,3))
7922 shat(n0,n1,n2,n3,:) = -c_0(n0,n1,n2,n3)
7925 shat(n0,n1,n2,n3,1) = shat(n0,n1,n2,n3,1) + c_i(n0,n2,n3,1)
7929 shat(n0,n1,n2,n3,2) = shat(n0,n1,n2,n3,2) + c_i(n0,n1,n3,2)
7933 shat(n0,n1,n2,n3,3) = shat(n0,n1,n2,n3,3) + c_i(n0,n1,n2,3)
7937 if(n0.eq.0.and.n1.eq.0.and.n2.eq.0.and.n3.eq.1)
then
7938 write(*,*)
'CalcDgy 0 C_0',c_0(n0,n1,n2,n3)
7939 write(*,*)
'CalcDgy 0 C_1',c_i(n0,n2,n3,1)
7940 write(*,*)
'CalcDgy 0 C_2',c_i(n0,n1,n3,2)
7941 write(*,*)
'CalcDgy 0 C_3',c_i(n0,n1,n2,3)
7942 write(*,*)
'CalcDgy 0 Sh1',shat(n0,n1,n2,n3,1)
7943 write(*,*)
'CalcDgy 0 Sh2',shat(n0,n1,n2,n3,2)
7944 write(*,*)
'CalcDgy 0 Sh3',shat(n0,n1,n2,n3,3)
7959 inds2 = reshape((/2,2,2,3,3,2,3,3/),shape(inds2))
7969 inds2 = reshape((/1,1,1,3,3,1,3,3/),shape(inds2))
7979 inds2 = reshape((/1,1,1,2,2,1,2,2/),shape(inds2))
7989 inds2 = reshape((/2,1,2,3,3,1,3,3/),shape(inds2))
7999 inds2 = reshape((/2,1,2,2,3,1,3,2/),shape(inds2))
8008 inds2 = reshape((/1,1,1,2,3,1,3,2/),shape(inds2))
8059 write(*,*)
'CalcDgy: Zadj',k,l,
zadj(k,l)
8060 write(*,*)
'CalcDgy: Xadj',a,b,
xadj(a,b)
8066 allocate(dexpgy(0:max(rmax/2,1),0:rmaxexp-2,0:rmaxexp-2,0:rmaxexp-2,0:ordgy_max))
8070 allocate(duvexpgy(0:rmaxexp,0:rmaxexp,0:rmaxexp,0:rmaxexp))
8072 duv(0:rmax,0:rmax,0:rmax,0:rmax) = duvexpgy(0:rmax,0:rmax,0:rmax,0:rmax)
8075 allocate(d00_err(0:rmaxexp))
8076 allocate(dij_err(0:rmaxexp))
8077 allocate(cij_err(0:rmaxc))
8079 allocate(d00_err2(0:rmaxexp))
8080 allocate(dij_err2(0:rmaxexp))
8081 allocate(cij_err2(0:rmaxc))
8087 cij_err = max(cerr_i(:,0),cerr_i(:,1),cerr_i(:,2),cerr_i(:,3))
8092 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,1),cerr2_i(:,2),cerr2_i(:,3))
8107 rloop:
do r=0,rmaxexp-2
8109 if (r.gt.rmax+2*gtrunc+2)
exit rloop
8116 maxdexpgy(1,r,0)=0d0
8128 daux =
zadj(k,1)*shat(0,inds(1),inds(2),inds(3),1) &
8129 +
zadj(k,2)*shat(0,inds(1),inds(2),inds(3),2) &
8130 +
zadj(k,3)*shat(0,inds(1),inds(2),inds(3),3)
8134 daux = daux - 2*nlt*
zadj(k,lt)*dexpgy(1,inds(1),inds(2),inds(3),0)
8140 daux = daux - 2*nltt*
zadj(k,ltt)*dexpgy(1,inds(1),inds(2),inds(3),0)
8143 dexpgy(1,inds0(1),inds0(2),inds0(3),0) = daux/(2*(nl+1)*
zadj(k,l))
8145 maxdexpgy(1,r,0) = maxdexpgy(1,r,0) + abs(dexpgy(1,inds0(1),inds0(2),inds0(3),0) )
8148 if (r+1.le.rmax)
then
8149 d(1,inds0(1),inds0(2),inds0(3)) = dexpgy(1,inds0(1),inds0(2),inds0(3),0)
8158 maxdexpgy(0,r,0)=0d0
8165 daux = (2d0*(1+r)*dexpgy(1,n1,n2,n3,0) - 4*duvexpgy(1,n1,n2,n3) &
8166 - c_0(0,n1,n2,n3))*
zadj(a,b)
8168 smod = shat(0,n1,n2,n3,:)
8171 if(n1.eq.0.and.n2.eq.2.and.n3.eq.0)
then
8172 write(*,*)
'CalcDgy 0 Smod',smod
8173 write(*,*)
'CalcDgy 0 Daux',daux
8178 smod(1) = smod(1) - 2d0*n1*dexpgy(1,n1-1,n2,n3,0)
8181 smod(2) = smod(2) - 2d0*n2*dexpgy(1,n1,n2-1,n3,0)
8184 smod(3) = smod(3) - 2d0*n3*dexpgy(1,n1,n2,n3-1,0)
8188 if(n1.eq.0.and.n2.eq.2.and.n3.eq.0)
then
8189 write(*,*)
'CalcDgy 0',r,a,b,
zadjf(b)/
xadj(a,b)
8191 write(*,*)
'CalcDgy 0 line1',r,daux/
xadj(a,b)
8198 daux = daux + zadj2(i)*
f(n)*smod(m)
8201 if(n1.eq.0.and.n2.eq.2.and.n3.eq.0)
then
8202 write(*,*)
'CalcDgy 0 2f',r,i,zadj2(i)*
f(n)*smod(m)/
xadj(a,b)
8208 dexpgy(0,n1,n2,n3,0) = daux/
xadj(a,b)
8211 if(n1.eq.1.and.n2.eq.1.and.n3.eq.1)
then
8212 write(*,*)
'CalcDgy D_0',r,dexpgy(0,n1,n2,n3,0)
8216 maxdexpgy(0,r,0) = maxdexpgy(0,r,0) + abs(dexpgy(0,n1,n2,n3,0))
8218 d(0,n1,n2,n3) = dexpgy(0,n1,n2,n3,0)
8227 derr(r) =
fac_gy*maxdexpgy(0,r,0)
8231 d00_err(r+2) = cij_err(r+1)/2d0
8232 dij_err(r)=max(
maxzadj/
maxxadj*max(2*(r+1)*d00_err(r+2),cerr_i(r,0)), &
8235 d00_err2(r+2) = cij_err2(r+1)/2d0
8236 dij_err2(r)=max(
maxzadj/
maxxadj*max(2*(r+1)*d00_err2(r+2),cerr2_i(r,0)), &
8244 gloop:
do g=1,min(gtrunc,r/2)
8248 maxdexpgy(1,rg,g) = 0d0
8258 daux = -
zadjf(k)*dexpgy(0,inds(1),inds(2),inds(3),g-1)
8261 daux = daux -
detz*dexpgy(0,inds(1),inds(2),inds(3),g-1)
8267 daux = daux - 2*nlt*
zadj(k,lt)*dexpgy(1,inds(1),inds(2),inds(3),g)
8273 daux = daux - 2*nltt*
zadj(k,ltt)*dexpgy(1,inds(1),inds(2),inds(3),g)
8276 dexpgy(1,inds0(1),inds0(2),inds0(3),g) = daux/(2*(nl+1)*
zadj(k,l))
8278 maxdexpgy(1,rg,g) = maxdexpgy(1,rg,g) + abs(dexpgy(1,inds0(1),inds0(2),inds0(3),g) )
8286 if (g.eq.1.and.abs(dexpgy(1,inds0(1),inds0(2),inds0(3),g)).gt. &
8287 1d1*truncfacexp*max(1/
m2scale,maxdexpgy(1,rg,g-1)) .or. &
8288 g.ge.2.and.abs(dexpgy(1,inds0(1),inds0(2),inds0(3),g)).gt. &
8289 truncfacexp*maxdexpgy(1,rg,g-1))
then
8292 write(*,*)
'CalcDgy exit gloop',1,inds0(1),inds0(2),inds0(3),g, &
8293 abs(dexpgy(1,inds0(1),inds0(2),inds0(3),g)),abs(dexpgy(1,inds0(1),inds0(2),inds0(3),g-1)),maxdexpgy(1,rg,g-1)
8307 if (rg+1.le.rmax)
then
8311 d(1,n1,n2,n3) = d(1,n1,n2,n3) + dexpgy(1,n1,n2,n3,g)
8318 maxdexpgy(0,rg,g) = 0d0
8327 daux = 2*(1+rg)*dexpgy(1,n1,n2,n3,g)*
zadj(a,b) &
8328 -
zadjf(b)*dexpgy(0,inds(1),inds(2),inds(3),g-1)
8332 smod(1) = smod(1) - 2d0*n1*dexpgy(1,n1-1,n2,n3,g)
8335 smod(2) = smod(2) - 2d0*n2*dexpgy(1,n1,n2-1,n3,g)
8338 smod(3) = smod(3) - 2d0*n3*dexpgy(1,n1,n2,n3-1,g)
8344 daux = daux + zadj2(i)*
f(n)*smod(m)
8347 dexpgy(0,n1,n2,n3,g) = daux/
xadj(a,b)
8349 maxdexpgy(0,rg,g) = maxdexpgy(0,rg,g) + abs(dexpgy(0,n1,n2,n3,g))
8355 if (g.eq.1.and.abs(dexpgy(0,n1,n2,n3,g)).gt. &
8356 truncfacexp*max(1/
m2scale**2,maxdexpgy(0,rg,g-1)) .or. &
8357 g.ge.2.and.abs(dexpgy(0,n1,n2,n3,g)).gt. &
8358 truncfacexp*maxdexpgy(0,rg,g-1))
then
8361 write(*,*)
'CalcDgy cycle loop',n1,n2,n3,g,abs(dexpgy(0,n1,n2,n3,g)),abs(dexpgy(0,n1,n2,n3,g-1)),maxdexpgy(0,rg,g-1)
8375 d00_err(rg+2) = max(d00_err(rg+2), &
8384 d00_err2(rg+2) = max(d00_err2(rg+2), &
8394 if (rg+2.le.rmax)
then
8398 d(1,n1,n2,n3) = d(1,n1,n2,n3) + dexpgy(1,n1,n2,n3,g)
8404 if ((rg.le.rmax))
then
8409 d(0,n1,n2,n3) = d(0,n1,n2,n3) + dexpgy(0,n1,n2,n3,g)
8410 if(abs(dexpgy(0,n1,n2,n3,g-1)).ne.0d0)
then
8412 derr(rg)=max(derr(rg),abs(dexpgy(0,n1,n2,n3,g))*min(1d0,abs(dexpgy(0,n1,n2,n3,g))/abs(dexpgy(0,n1,n2,n3,g-1))))
8414 derr(rg)=max(derr(rg),abs(dexpgy(0,n1,n2,n3,g)))
8427 if(dij_err2(rg).gt.3d0*derr(rg))
then
8429 if(dij_err(rg).gt.3d0*derr(rg))
then
8431 gtrunc = min(g,gtrunc)
8435 write(*,*)
'CalcDgy exit err',rg,g,gtrunc
8436 write(*,*)
'CalcDgy exit err',dij_err(rg),derr(rg)
8446 write(*,*)
'CalcDgy D(1,0,0,0)',r,d(1,0,0,0)
8447 write(*,*)
'CalcDgy D(0,0,0,0)',r,d(0,0,0,0)
8448 write(*,*)
'CalcDgy D(0,0,0,1)',r,d(0,0,0,1)
8449 if (r.ge.2.and.rmax.ge.2)
then
8450 write(*,*)
'CalcDgy D(0,0,0,2)',r,d(0,0,0,2)
8452 if (r.ge.3.and.rmax.ge.3)
then
8453 write(*,*)
'CalcDgy D(1,0,0,1)',r,d(1,0,0,1)
8454 write(*,*)
'CalcDgy D(0,1,0,2)',r,d(0,1,0,2)
8455 write(*,*)
'CalcDgy D(0,0,0,3)',r,d(0,0,0,3)
8456 write(*,*)
'CalcDgy D(0,1,1,1)',r,d(0,1,1,1)
8457 write(*,*)
'CalcDgy D(0,2,1,0)',r,d(0,2,1,0)
8460 write(*,*)
'CalcDgy Dij_err',r,dij_err
8461 write(*,*)
'CalcDgy Dij_acc',r,dij_err/abs(d(0,0,0,0))
8463 write(*,*)
'CalcDgy err',r,g,derr
8464 write(*,*)
'CalcDgy acc',r,g,derr/abs(d(0,0,0,0))
8467 derr2 = max(derr,dij_err2(0:rmax))
8468 derr = max(derr,dij_err(0:rmax))
8473 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0)
then
8478 d(0,n1,n2,rg-n1-n2)=0d0
8485 d(1,n1,n2,rg-2-n1-n2)=0d0
8491 111
format(a22,2(
'(',g24.17,
',',g24.17,
') ':))
8492 call seterrflag_coli(-5)
8493 call errout_coli(
'CalcDgy',
' exit rloop for D', &
8495 if (errorwriteflag)
then
8496 write(nerrout_coli,100)
' CalcDgy: exit rloop for D ', &
8497 ' should not appear'
8498 write(nerrout_coli,111)
' CalcDgy: p10 = ',p10
8499 write(nerrout_coli,111)
' CalcDgy: p21 = ',p21
8500 write(nerrout_coli,111)
' CalcDgy: p32 = ',p32
8501 write(nerrout_coli,111)
' CalcDgy: p30 = ',p30
8502 write(nerrout_coli,111)
' CalcDgy: p20 = ',p20
8503 write(nerrout_coli,111)
' CalcDgy: p31 = ',p31
8504 write(nerrout_coli,111)
' CalcDgy: m02 = ',m02
8505 write(nerrout_coli,111)
' CalcDgy: m12 = ',m12
8506 write(nerrout_coli,111)
' CalcDgy: m22 = ',m22
8507 write(nerrout_coli,111)
' CalcDgy: m32 = ',m32
8512 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0.and.r.ge.rmax)
then
8527 do n0=2,max(rmax,r/2)
8529 do nlt=r-2*n0-nl,0,-1
8530 nltt = r-2*n0-nl-nlt
8538 daux =
zadj(k,1)*shat(n0-1,inds(1),inds(2),inds(3),1) &
8539 +
zadj(k,2)*shat(n0-1,inds(1),inds(2),inds(3),2) &
8540 +
zadj(k,3)*shat(n0-1,inds(1),inds(2),inds(3),3) &
8541 -
zadjf(k)*d(n0-1,inds(1),inds(2),inds(3))
8544 daux = daux -
detz*d(n0-1,inds(1),inds(2),inds(3))
8549 daux = daux - 2*nlt*
zadj(k,lt)*d(n0,inds(1),inds(2),inds(3))
8554 daux = daux - 2*nltt*
zadj(k,ltt)*d(n0,inds(1),inds(2),inds(3))
8557 d(n0,inds0(1),inds0(2),inds0(3)) = daux/(2*(nl+1)*
zadj(k,l))
8574 write(*,*)
'CalcDgy exp rmax+1',r,n0,n1,n2,n3, d(n0,n1,n2,n3)
8576 d(n0,n1,n2,n3) = (c_0(n0-1,n1,n2,n3) + 2*
mm02*d(n0-1,n1,n2,n3) &
8577 + 4*duv(n0,n1,n2,n3) &
8578 +
f(1)*d(n0-1,n1+1,n2,n3) +
f(2)*d(n0-1,n1,n2+1,n3) &
8579 +
f(3)*d(n0-1,n1,n2,n3+1)) / (2*(r-1))
8581 write(*,*)
'CalcDgy dir rmax+1',r,n0,n1,n2,n3, d(n0,n1,n2,n3)
8591 write(*,*)
'CalcDgy D(1,0,0,0) fin',d(1,0,0,0)
8592 write(*,*)
'CalcDgy D(0,0,2,0) fin',d(0,0,2,0)
8593 write(*,*)
'CalcDgy D(0,0,0,2) fin',d(0,0,0,2)
8595 write(*,*)
'CalcDgy D(1,0,1,0) fin',d(1,0,1,0)
8596 write(*,*)
'CalcDgy D(0,1,1,1) fin',d(0,1,1,1)
8597 write(*,*)
'CalcDgy D(0,0,3,0) fin',d(0,0,3,0)
8601 write(*,*)
'CalcDgy final err',derr
8602 write(*,*)
'CalcDgy final acc',derr/abs(d(0,0,0,0))
8620 subroutine calcdgp(D,Duv,p10,p21,p32,p30,p20,p31, &
8621 m02,m12,m22,m32,rmax,ordgp_min,ordgp_max,id,Derr,Derr2)
8625 integer,
intent(in) :: rmax,ordgp_min,ordgp_max,id
8626 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
8627 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
8628 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
8629 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
8630 double complex,
allocatable :: Dexpgp(:,:,:,:,:), DuvExpgp(:,:,:,:)
8631 double complex,
allocatable :: C_0(:,:,:,:), Cuv_0(:,:,:,:), Shat(:,:,:)
8632 double complex,
allocatable :: C_k(:,:,:), Cuv_k(:,:,:)
8633 double complex,
allocatable :: D_alt(:,:,:,:)
8634 double precision,
allocatable :: Cerr_i(:,:),Cerr2_i(:,:)
8635 double complex :: Smod, fk, elimminf2_coli
8636 double precision,
allocatable :: D00_err(:),Dij_err(:),Cij_err(:),acc_req_Cextra(:)
8637 double precision,
allocatable :: D00_err2(:),Dij_err2(:),Cij_err2(:)
8638 double precision :: maxDexpgp(0:1,0:rmax+ordgp_min+1,0:ordgp_max),truncfacexp
8639 integer :: rmaxC,rmaxExp,gtrunc,r,n0,n1,n2,n3,k,l,g,rg
8640 integer :: bin,nid(0:3),i
8641 logical :: errorwriteflag
8644 write(*,*)
'CalcDgp in, ord',rmax,ordgp_min,ordgp_max
8676 if (abs(
f(1)).ge.max(abs(
f(2)),abs(
f(3))))
then
8678 else if (abs(
f(2)).ge.max(abs(
f(1)),abs(
f(3))))
then
8687 rmaxc = rmax + ordgp_min
8688 allocate(c_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
8689 allocate(cuv_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
8690 allocate(c_k(0:rmaxc,0:rmaxc,0:rmaxc))
8691 allocate(cuv_k(0:rmaxc,0:rmaxc,0:rmaxc))
8692 allocate(cerr_i(0:rmaxc,0:3))
8693 allocate(cerr2_i(0:rmaxc,0:3))
8694 allocate(acc_req_cextra(0:rmaxc))
8700 if (mod(id/bin,2).eq.0)
then
8709 acc_req_cextra(0:rmax) = acc_req_cind
8710 if(
w_gp.ne.0d0)
then
8712 acc_req_cextra(r)= acc_req_cextra(r-1)/
w_gp
8715 acc_req_cextra(rmax+1:rmaxc)=acc_inf
8718 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0),rmax,acc_req_cextra)
8720 call calcc(c_k(:,:,:),cuv_k(:,:,:),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1),rmax,acc_req_cextra)
8721 else if (k.eq.2)
then
8722 call calcc(c_k(:,:,:),cuv_k(:,:,:),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2),rmax,acc_req_cextra)
8723 else if (k.eq.3)
then
8724 call calcc(c_k(:,:,:),cuv_k(:,:,:),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3),rmax,acc_req_cextra)
8731 n0 = (rmaxc-n1-n2-n3)
8732 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
8733 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
8734 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
8735 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
8742 allocate(shat(0:rmaxc,0:rmaxc,0:rmaxc))
8749 shat(n1,n2,n3) = -c_0(0,n1,n2,n3)
8751 if ((k.eq.1).and.(n1.eq.0))
then
8752 shat(n1,n2,n3) = shat(n1,n2,n3) + c_k(0,n2,n3)
8753 else if ((k.eq.2).and.(n2.eq.0))
then
8754 shat(n1,n2,n3) = shat(n1,n2,n3) + c_k(0,n1,n3)
8755 else if ((k.eq.3).and.(n3.eq.0))
then
8756 shat(n1,n2,n3) = shat(n1,n2,n3) + c_k(0,n1,n2)
8767 allocate(dexpgp(0:rmaxexp/2,0:rmaxexp,0:rmaxexp,0:rmaxexp,0:ordgp_max))
8770 allocate(duvexpgp(0:rmaxexp,0:rmaxexp,0:rmaxexp,0:rmaxexp))
8772 duv(0:rmax,0:rmax,0:rmax,0:rmax) = duvexpgp(0:rmax,0:rmax,0:rmax,0:rmax)
8775 allocate(d00_err(0:rmaxexp))
8776 allocate(dij_err(0:rmaxexp))
8777 allocate(cij_err(0:rmaxc))
8779 allocate(d00_err2(0:rmaxexp))
8780 allocate(dij_err2(0:rmaxexp))
8781 allocate(cij_err2(0:rmaxc))
8795 cij_err = max(cerr_i(:,0),cerr_i(:,k))
8796 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,k))
8808 rloop:
do r=1,rmaxexp
8810 if (r.gt.rmax+gtrunc+1)
exit rloop
8819 maxdexpgp(1,r,0)=0d0
8825 dexpgp(n0,n1,n2,n3,0) = (2d0*duvexpgp(n0,n1,n2,n3) + c_0(n0-1,n1,n2,n3) &
8826 +
mm02*dexpgp(n0-1,n1,n2,n3,0))/((r-n0)+1)/2d0
8829 maxdexpgp(1,r,0) = maxdexpgp(1,r,0) + abs(dexpgp(n0,n1,n2,n3,0) )
8832 if (r-n0.le.rmax)
then
8833 d(n0,n1,n2,n3) = dexpgp(n0,n1,n2,n3,0)
8841 write(*,*)
'CalcDgp 0 D(1,0,1,0)= ',r,d(1,0,1,0)
8848 maxdexpgp(0,r-1,0)=0d0
8853 smod = shat(n1,n2,n3)
8854 if ((k.eq.1).and.(n1.ge.1))
then
8855 smod = smod - 2d0*n1*dexpgp(1,n1-1,n2,n3,0)
8856 else if ((k.eq.2).and.(n2.ge.1))
then
8857 smod = smod - 2d0*n2*dexpgp(1,n1,n2-1,n3,0)
8858 else if ((k.eq.3).and.(n3.ge.1))
then
8859 smod = smod - 2d0*n3*dexpgp(1,n1,n2,n3-1,0)
8862 dexpgp(0,n1,n2,n3,0) = smod/fk
8863 maxdexpgp(0,r-1,0) = maxdexpgp(0,r-1,0) + abs(dexpgp(0,n1,n2,n3,0))
8865 if (r.le.rmax+1)
then
8866 d(0,n1,n2,n3) = dexpgp(0,n1,n2,n3,0)
8873 if (r.le.rmax+1)
then
8875 derr(r-1) =
fac_gp*maxdexpgp(0,r-1,0)
8880 d00_err(r) = cij_err(r-2)/(2*r)
8882 dij_err(r-1)=max(cij_err(r-1),2*d00_err(r))/abs(fk)
8885 d00_err2(r) = cij_err2(r-2)/(2*r)
8887 dij_err2(r-1)=max(cij_err2(r-1),2*d00_err2(r))/abs(fk)
8894 gloop:
do g=1,min(gtrunc,r-1)
8900 maxdexpgp(1,rg,g) = 0d0
8906 dexpgp(n0,n1,n2,n3,g) = (2d0*
mm02*dexpgp(n0-1,n1,n2,n3,g) &
8907 -
z(1,1)*dexpgp(n0-1,n1+2,n2,n3,g-1) - 2d0*
z(2,1)*dexpgp(n0-1,n1+1,n2+1,n3,g-1) &
8908 - 2d0*
z(3,1)*dexpgp(n0-1,n1+1,n2,n3+1,g-1) -
z(2,2)*dexpgp(n0-1,n1,n2+2,n3,g-1) &
8909 - 2d0*
z(3,2)*dexpgp(n0-1,n1,n2+1,n3+1,g-1) -
z(3,3)*dexpgp(n0-1,n1,n2,n3+2,g-1)) &
8913 maxdexpgp(1,rg,g) = maxdexpgp(1,rg,g) + abs(dexpgp(n0,n1,n2,n3,g))
8916 if(n0.eq.1.and.n1.eq.0.and.n2.eq.1.and.n3.eq.0)
then
8917 write(*,*)
'CalcDgp Dexp(1,0,1,0,g)',r,rg,g,dexpgp(1,0,1,0,g)
8918 write(*,*)
'CalcDgp D(1,0,1,0)',r,rg,g,d(1,0,1,0)
8919 write(*,*)
'CalcDgp maxDexpgp(1,rg,g)',r,rg,g,maxdexpgp(1,rg,g)
8920 if(g.gt.0)
write(*,*)
'CalcDgp trunc',abs(dexpgp(n0,n1,n2,n3,g)), &
8921 truncfacexp*maxdexpgp(1,rg,g-1),truncfacexp,maxdexpgp(1,rg,g-1)
8925 if (g.eq.1.and.abs(dexpgp(1,n1,n2,n3,g)).gt. &
8926 truncfacexp*max(1/
m2scale,maxdexpgp(1,rg,g-1)) .or. &
8927 g.ge.2.and.abs(dexpgp(1,n1,n2,n3,g)).gt. &
8928 truncfacexp*maxdexpgp(1,rg,g-1))
then
8942 if (rg-n0.le.rmax)
then
8946 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) + dexpgp(n0,n1,n2,n3,g)
8956 maxdexpgp(0,rg-1,g) = 0d0
8961 smod = -
z(1,k)*dexpgp(0,n1+1,n2,n3,g-1) &
8962 -
z(2,k)*dexpgp(0,n1,n2+1,n3,g-1) &
8963 -
z(3,k)*dexpgp(0,n1,n2,n3+1,g-1)
8964 if ((k.eq.1).and.(n1.ge.1))
then
8965 smod = smod - 2d0*n1*dexpgp(1,n1-1,n2,n3,g)
8966 else if ((k.eq.2).and.(n2.ge.1))
then
8967 smod = smod - 2d0*n2*dexpgp(1,n1,n2-1,n3,g)
8968 else if ((k.eq.3).and.(n3.ge.1))
then
8969 smod = smod - 2d0*n3*dexpgp(1,n1,n2,n3-1,g)
8972 dexpgp(0,n1,n2,n3,g) = smod/fk
8974 maxdexpgp(0,rg-1,g) = maxdexpgp(0,rg-1,g) + abs(dexpgp(0,n1,n2,n3,g))
8976 if (g.eq.1.and.abs(dexpgp(0,n1,n2,n3,g)).gt. &
8977 truncfacexp*max(1/
m2scale**2,maxdexpgp(0,rg-1,g-1)) .or. &
8978 g.ge.2.and.abs(dexpgp(0,n1,n2,n3,g)).gt. &
8979 truncfacexp*maxdexpgp(0,rg-1,g-1))
then
8990 d00_err(rg) = max(d00_err(rg),max(2*abs(m02)*dij_err(rg-2),
maxz*dij_err(rg))/(4*r))
8992 dij_err(rg-1) = max(dij_err(rg-1),max(2*d00_err(rg),
maxz*dij_err(rg))/abs(fk))
8995 d00_err2(rg) = max(d00_err2(rg),max(2*abs(m02)*dij_err2(rg-2),
maxz*dij_err2(rg))/(4*r))
8997 dij_err2(rg-1) = max(dij_err2(rg-1),max(2*d00_err2(rg),
maxz*dij_err2(rg))/abs(fk))
9001 if (rg-n0.le.rmax)
then
9005 d(n0,n1,n2,n3) = d(n0,n1,n2,n3) + dexpgp(n0,n1,n2,n3,g)
9012 if ((rg.le.rmax+1))
then
9017 d(0,n1,n2,n3) = d(0,n1,n2,n3) + dexpgp(0,n1,n2,n3,g)
9018 if(abs(dexpgp(0,n1,n2,n3,g-1)).ne.0d0)
then
9020 derr(rg-1)=max(derr(rg-1),abs(dexpgp(0,n1,n2,n3,g))*min(1d0,abs(dexpgp(0,n1,n2,n3,g))/abs(dexpgp(0,n1,n2,n3,g-1))))
9022 derr(rg-1)=max(derr(rg-1),abs(dexpgp(0,n1,n2,n3,g)))
9029 if(dij_err2(rg-1).gt.3d0*derr(rg-1))
then
9031 if(dij_err(rg-1).gt.3d0*derr(rg-1))
then
9033 gtrunc = min(g,gtrunc)
9036 write(*,*)
'CalcDgp exit err',r,g,gtrunc
9045 write(*,*)
'CalcDgp D(0,0,0,0)',r,d(0,0,0,0)
9046 write(*,*)
'CalcDgp D(0,0,1,0)',r,d(0,0,1,0)
9047 write(*,*)
'CalcDgp D(1,0,0,0)',r,d(1,0,0,0)
9048 write(*,*)
'CalcDgp D(1,0,1,0)',r,d(1,0,1,0)
9049 write(*,*)
'CalcDgp D(0,0,3,0)',r,d(0,0,3,0)
9051 write(*,*)
'CalcDgp Dij_err',r,dij_err
9052 write(*,*)
'CalcDgp Dij_acc',r,dij_err/abs(d(0,0,0,0))
9054 write(*,*)
'CalcDgp err',r,derr
9055 write(*,*)
'CalcDgp acc',r,derr/abs(d(0,0,0,0))
9058 derr2 = max(derr,dij_err2(0:rmax))
9059 derr = max(derr,dij_err(0:rmax))
9064 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0)
then
9071 d(n0,n1,n2,rg-2*n0-n1-n2)=0d0
9079 d(0,n1,n2,r-n1-n2)=0d0
9085 111
format(a22,2(
'(',g24.17,
',',g24.17,
') ':))
9086 call seterrflag_coli(-5)
9087 call errout_coli(
'CalcDgp',
' exit rloop for D', &
9089 if (errorwriteflag)
then
9090 write(nerrout_coli,100)
' CalcDgp: exit rloop for D ', &
9091 ' should not appear'
9092 write(nerrout_coli,111)
' CalcDgp: p10 = ',p10
9093 write(nerrout_coli,111)
' CalcDgp: p21 = ',p21
9094 write(nerrout_coli,111)
' CalcDgp: p32 = ',p32
9095 write(nerrout_coli,111)
' CalcDgp: p30 = ',p30
9096 write(nerrout_coli,111)
' CalcDgp: p20 = ',p20
9097 write(nerrout_coli,111)
' CalcDgp: p31 = ',p31
9098 write(nerrout_coli,111)
' CalcDgp: m02 = ',m02
9099 write(nerrout_coli,111)
' CalcDgp: m12 = ',m12
9100 write(nerrout_coli,111)
' CalcDgp: m22 = ',m22
9101 write(nerrout_coli,111)
' CalcDgp: m32 = ',m32
9106 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0.and.r.ge.rmax)
then
9123 write(*,*)
'CalcDgp exp rmax+1',r,n0,n1,n2,n3, d(n0,n1,n2,n3)
9125 d(n0,n1,n2,n3) = (c_0(n0-1,n1,n2,n3) + 2*
mm02*d(n0-1,n1,n2,n3) &
9126 + 4*duv(n0,n1,n2,n3) &
9127 +
f(1)*d(n0-1,n1+1,n2,n3) +
f(2)*d(n0-1,n1,n2+1,n3) &
9128 +
f(3)*d(n0-1,n1,n2,n3+1)) / (2*(r-1))
9130 write(*,*)
'CalcDgp exp rmax+1',r,n0,n1,n2,n3, d(n0,n1,n2,n3)
9139 write(*,*)
'CalcDgp D(1,0,0,0) fin',d(1,0,0,0)
9140 write(*,*)
'CalcDgp D(1,0,1,0) fin',d(1,0,1,0)
9141 write(*,*)
'CalcDgp D(0,0,3,0) fin',d(0,0,3,0)
9142 write(*,*)
'CalcDgp D(0,1,1,1) fin',d(0,1,1,1)
9144 write(*,*)
'CalcDgp final err',derr
9145 write(*,*)
'CalcDgp final acc',derr/abs(d(0,0,0,0))
9164 subroutine calcdgpf(D,Duv,p10,p21,p32,p30,p20,p31, &
9165 m02,m12,m22,m32,rmax,ordgpf_min,ordgpf_max,id,Derr,Derr2)
9169 integer,
intent(in) :: rmax,ordgpf_min,ordgpf_max,id
9170 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
9171 double complex ::Zadj2(4)
9172 double complex,
allocatable :: Dexpgpf(:,:,:,:,:), DuvExpgpf(:,:,:,:)
9173 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
9174 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
9175 double precision,
intent(out) :: Derr(0:rmax),Derr2(0:rmax)
9176 double complex,
allocatable :: C_0(:,:,:,:), C_i(:,:,:,:), Shat(:,:,:,:,:)
9177 double complex,
allocatable :: Cuv_0(:,:,:,:), Cuv_i(:,:,:,:)
9178 double complex,
allocatable :: D_alt(:,:,:,:)
9179 double precision,
allocatable :: Cerr_i(:,:),Cerr2_i(:,:)
9180 double complex :: Smod(3), Daux, elimminf2_coli
9181 double precision,
allocatable :: D00_err(:),Dij_err(:),Cij_err(:),acc_req_Cextra(:)
9182 double precision,
allocatable :: D00_err2(:),Dij_err2(:),Cij_err2(:)
9183 double precision :: maxDexpgpf(0:1,0:rmax+2*ordgpf_min,0:ordgpf_max),truncfacexp,acc_aux
9184 double precision :: minZk
9185 integer :: rmaxC,rmaxExp,gtrunc,r,n0,n1,n2,n3,a,b,i,j,g,rg,m,n
9186 integer :: inds0(3),inds(3),inds2(2,4),at,bt,k,l,lt,ltt,nl,nlt,nltt
9187 integer :: bin,nid(0:3)
9188 logical :: errorwriteflag
9191 write(*,*)
'CalcDgpf in, ord',rmax,ordgpf_min,ordgpf_max
9195 rmaxc = rmax + 2*ordgpf_min + 1
9196 allocate(c_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
9197 allocate(cuv_0(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc))
9198 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
9199 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,3))
9200 allocate(cerr_i(0:rmaxc,0:3))
9201 allocate(cerr2_i(0:rmaxc,0:3))
9202 allocate(acc_req_cextra(0:rmaxc))
9208 if (mod(id/bin,2).eq.0)
then
9217 acc_req_cextra(0:rmax+1) = acc_req_cind
9219 if (
y_gpf.ne.0d0)
then
9221 acc_req_cextra(rmax+2*g) = acc_req_cextra(rmax+2*g-2)/
y_gpf
9222 acc_req_cextra(rmax+2*g+1) = acc_req_cextra(rmax+2*g-1)/
y_gpf
9224 acc_req_cextra(rmax+g+1) = min(acc_req_cextra(rmax+g+1),acc_aux)
9226 else if(
x_gpf.ne.0d0)
then
9228 acc_aux = acc_aux/
x_gpf
9229 acc_req_cextra(rmax+g+1) = acc_aux
9232 acc_req_cextra(rmax+2:rmax+2*ordgpf_min+1) = acc_inf
9238 write(*,*)
'CalcDgpf: accreq_Cextra',acc_req_cextra
9241 call calcc(c_0(:,0,:,:),cuv_0(:,0,:,:),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(:,0),cerr2_i(:,0),rmax,acc_req_cextra)
9242 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(:,1),cerr2_i(:,1),rmax,acc_req_cextra)
9243 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(:,2),cerr2_i(:,2),rmax,acc_req_cextra)
9244 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(:,3),cerr2_i(:,3),rmax,acc_req_cextra)
9247 write(*,*)
'CalcDgpf Cerr 0',cerr_i(:,0)
9248 write(*,*)
'CalcDgpf Cerr 1',cerr_i(:,1)
9249 write(*,*)
'CalcDgpf Cerr 2',cerr_i(:,2)
9250 write(*,*)
'CalcDgpf Cerr 3',cerr_i(:,3)
9258 n0 = (rmaxc-n1-n2-n3)
9259 c_0(0:n0,n1,n2,n3) = -c_0(0:n0,n1-1,n2,n3) &
9260 -c_0(0:n0,n1-1,n2+1,n3)-c_0(0:n0,n1-1,n2,n3+1)
9261 cuv_0(0:n0,n1,n2,n3) = -cuv_0(0:n0,n1-1,n2,n3) &
9262 -cuv_0(0:n0,n1-1,n2+1,n3)-cuv_0(0:n0,n1-1,n2,n3+1)
9271 allocate(shat(0:rmaxc,0:rmaxc,0:rmaxc,0:rmaxc,3))
9279 shat(n0,n1,n2,n3,:) = -c_0(n0,n1,n2,n3)
9282 shat(n0,n1,n2,n3,1) = shat(n0,n1,n2,n3,1) + c_i(n0,n2,n3,1)
9286 shat(n0,n1,n2,n3,2) = shat(n0,n1,n2,n3,2) + c_i(n0,n1,n3,2)
9290 shat(n0,n1,n2,n3,3) = shat(n0,n1,n2,n3,3) + c_i(n0,n1,n2,3)
9294 if(n0.eq.0.and.n1.eq.0.and.n2.eq.0.and.n3.eq.1)
then
9295 write(*,*)
'CalcDgpf 0 C_0',c_0(n0,n1,n2,n3)
9296 write(*,*)
'CalcDgpf 0 C_1',c_i(n0,n2,n3,1)
9297 write(*,*)
'CalcDgpf 0 C_2',c_i(n0,n1,n3,2)
9298 write(*,*)
'CalcDgpf 0 C_3',c_i(n0,n1,n2,3)
9299 write(*,*)
'CalcDgpf 0 Sh1',shat(n0,n1,n2,n3,1)
9300 write(*,*)
'CalcDgpf 0 Sh2',shat(n0,n1,n2,n3,2)
9301 write(*,*)
'CalcDgpf 0 Sh3',shat(n0,n1,n2,n3,3)
9312 if (maxval(abs(
z(1,1:3))).le.minzk)
then
9313 minzk = maxval(abs(
z(1,1:3)))
9319 if (maxval(abs(
z(2,1:3))).lt.minzk)
then
9320 minzk = maxval(abs(
z(2,1:3)))
9326 if (maxval(abs(
z(3,1:3))).lt.minzk)
then
9327 minzk = maxval(abs(
z(3,1:3)))
9335 write(*,*)
'CalcDgpf: Z',k, maxval(abs(
z(k,1:3)))
9341 allocate(dexpgpf(0:max(rmax/2,1),0:rmaxexp-2,0:rmaxexp-2,0:rmaxexp-2,0:ordgpf_max))
9345 allocate(duvexpgpf(0:rmaxexp,0:rmaxexp,0:rmaxexp,0:rmaxexp))
9347 duv(0:rmax,0:rmax,0:rmax,0:rmax) = duvexpgpf(0:rmax,0:rmax,0:rmax,0:rmax)
9350 allocate(d00_err(0:rmaxexp))
9351 allocate(dij_err(0:rmaxexp))
9352 allocate(cij_err(0:rmaxc))
9354 allocate(d00_err2(0:rmaxexp))
9355 allocate(dij_err2(0:rmaxexp))
9356 allocate(cij_err2(0:rmaxc))
9362 cij_err = max(cerr_i(:,0),cerr_i(:,1),cerr_i(:,2),cerr_i(:,3))
9367 cij_err2 = max(cerr2_i(:,0),cerr2_i(:,1),cerr2_i(:,2),cerr2_i(:,3))
9376 rloop:
do r=0,rmaxexp-2
9378 if (r.gt.rmax+2*gtrunc+2)
exit rloop
9385 maxdexpgpf(1,r,0)=0d0
9397 daux = shat(0,inds(1),inds(2),inds(3),k)
9399 dexpgpf(1,inds0(1),inds0(2),inds0(3),0) = daux/(2*(nl+1))
9401 maxdexpgpf(1,r,0) = maxdexpgpf(1,r,0) + abs(dexpgpf(1,inds0(1),inds0(2),inds0(3),0) )
9404 if (r+1.le.rmax)
then
9405 d(1,inds0(1),inds0(2),inds0(3)) = dexpgpf(1,inds0(1),inds0(2),inds0(3),0)
9414 maxdexpgpf(0,r,0)=0d0
9419 daux = 2d0*(4+r+r)*dexpgpf(1,n1,n2,n3,0) - 4*duvexpgpf(1,n1,n2,n3) &
9423 if(n1.eq.0.and.n2.eq.2.and.n3.eq.0)
then
9424 write(*,*)
'CalcDgpf 0 Daux',daux
9428 dexpgpf(0,n1,n2,n3,0) = daux/(2d0*m02)
9431 if(n1.eq.1.and.n2.eq.1.and.n3.eq.1)
then
9432 write(*,*)
'CalcDgpf D_0',r,dexpgpf(0,n1,n2,n3,0)
9436 maxdexpgpf(0,r,0) = maxdexpgpf(0,r,0) + abs(dexpgpf(0,n1,n2,n3,0))
9438 d(0,n1,n2,n3) = dexpgpf(0,n1,n2,n3,0)
9447 derr(r) =
fac_gpf*maxdexpgpf(0,r,0)
9451 d00_err(r+2) = cij_err(r+1)/2d0
9452 dij_err(r)=1d0/abs(m02)*max((2*r+4)*d00_err(r+2),cerr_i(r,0))
9454 d00_err2(r+2) = cij_err2(r+1)/2d0
9455 dij_err2(r)=1d0/abs(m02)*max((2*r+4)*d00_err2(r+2),cerr2_i(r,0))
9462 gloop:
do g=1,min(gtrunc,r/2)
9466 maxdexpgpf(1,rg,g) = 0d0
9476 daux = -
f(k)*dexpgpf(0,inds(1),inds(2),inds(3),g-1)
9479 daux = daux -
z(k,l)*dexpgpf(0,inds(1),inds(2),inds(3),g-1)
9482 inds(lt) = inds(lt)+1
9483 daux = daux -
z(k,lt)*dexpgpf(0,inds(1),inds(2),inds(3),g-1)
9485 inds(lt) = inds(lt)-1
9486 inds(ltt) = inds(ltt)+1
9487 daux = daux -
z(k,ltt)*dexpgpf(0,inds(1),inds(2),inds(3),g-1)
9489 dexpgpf(1,inds0(1),inds0(2),inds0(3),g) = daux/(2*(nl+1))
9491 maxdexpgpf(1,rg,g) = maxdexpgpf(1,rg,g) + abs(dexpgpf(1,inds0(1),inds0(2),inds0(3),g) )
9498 if (g.eq.1.and.abs(dexpgpf(1,inds0(1),inds0(2),inds0(3),g)).gt. &
9499 1d1*truncfacexp*max(1/
m2scale,maxdexpgpf(1,rg,g-1)) .or. &
9500 g.ge.2.and.abs(dexpgpf(1,inds0(1),inds0(2),inds0(3),g)).gt. &
9501 truncfacexp*maxdexpgpf(1,rg,g-1))
then
9504 write(*,*)
'CalcDgpf exit gloop',1,inds0(1),inds0(2),inds0(3),g, &
9505 abs(dexpgpf(1,inds0(1),inds0(2),inds0(3),g)),abs(dexpgpf(1,inds0(1),inds0(2),inds0(3),g-1)),maxdexpgpf(1,rg,g-1)
9519 if (rg+1.le.rmax)
then
9523 d(1,n1,n2,n3) = d(1,n1,n2,n3) + dexpgpf(1,n1,n2,n3,g)
9530 maxdexpgpf(0,rg,g) = 0d0
9538 daux = 2*(4+rg+rg)*dexpgpf(1,n1,n2,n3,g)
9544 daux = daux +
z(i,j)*dexpgpf(0,inds(1),inds(2),inds(3),g-1)
9550 dexpgpf(0,n1,n2,n3,g) = daux/(2*m02)
9552 maxdexpgpf(0,rg,g) = maxdexpgpf(0,rg,g) + abs(dexpgpf(0,n1,n2,n3,g))
9558 if (g.eq.1.and.abs(dexpgpf(0,n1,n2,n3,g)).gt. &
9559 truncfacexp*max(1/
m2scale**2,maxdexpgpf(0,rg,g-1)) .or. &
9560 g.ge.2.and.abs(dexpgpf(0,n1,n2,n3,g)).gt. &
9561 truncfacexp*maxdexpgpf(0,rg,g-1))
then
9564 write(*,*)
'CalcDgpf exit gloop',n1,n2,n3,g,abs(dexpgpf(0,n1,n2,n3,g)),abs(dexpgpf(0,n1,n2,n3,g-1)),maxdexpgpf(0,rg,g-1)
9578 d00_err(rg+2) = max(d00_err(rg+2), &
9579 fmax/2d0*dij_err(rg+1), &
9580 maxz/2d0*dij_err(rg+2))
9582 dij_err(rg)=max(dij_err(rg),
maxz/(2*abs(m02))*dij_err(rg+2), &
9583 (2*rg+4)/abs(m02)*d00_err(rg+2))
9586 d00_err2(rg+2) = max(d00_err2(rg+2), &
9587 fmax/2d0*dij_err2(rg+1), &
9588 maxz/2d0*dij_err2(rg+2))
9590 dij_err2(rg)=max(dij_err2(rg),
maxz/(2*abs(m02))*dij_err2(rg+2), &
9591 (2*rg+4)/abs(m02)*d00_err2(rg+2))
9594 if (rg+2.le.rmax)
then
9598 d(1,n1,n2,n3) = d(1,n1,n2,n3) + dexpgpf(1,n1,n2,n3,g)
9604 if (rg.le.rmax)
then
9609 d(0,n1,n2,n3) = d(0,n1,n2,n3) + dexpgpf(0,n1,n2,n3,g)
9610 if(abs(dexpgpf(0,n1,n2,n3,g-1)).ne.0d0)
then
9611 derr(rg)=max(derr(rg),abs(dexpgpf(0,n1,n2,n3,g))*min(1d0,abs(dexpgpf(0,n1,n2,n3,g))/abs(dexpgpf(0,n1,n2,n3,g-1))))
9613 derr(rg)=max(derr(rg),abs(dexpgpf(0,n1,n2,n3,g)))
9626 if(dij_err2(rg).gt.3d0*derr(rg))
then
9628 if(dij_err(rg).gt.3d0*derr(rg))
then
9630 gtrunc = min(g,gtrunc)
9634 write(*,*)
'CalcDgpf exit err',rg,g,gtrunc
9635 write(*,*)
'CalcDgpf exit err',dij_err(rg),derr(rg)
9645 write(*,*)
'CalcDgpf D(1,0,0,0)',r,d(1,0,0,0)
9646 write(*,*)
'CalcDgpf D(0,0,0,0)',r,d(0,0,0,0)
9647 write(*,*)
'CalcDgpf D(0,0,0,1)',r,d(0,0,0,1)
9648 if (r.ge.2.and.rmax.ge.2)
then
9649 write(*,*)
'CalcDgpf D(0,0,0,2)',r,d(0,0,0,2)
9651 if (r.ge.3.and.rmax.ge.3)
then
9652 write(*,*)
'CalcDgpf D(1,0,0,1)',r,d(1,0,0,1)
9653 write(*,*)
'CalcDgpf D(0,1,0,2)',r,d(0,1,0,2)
9654 write(*,*)
'CalcDgpf D(0,0,0,3)',r,d(0,0,0,3)
9655 write(*,*)
'CalcDgpf D(0,1,1,1)',r,d(0,1,1,1)
9656 write(*,*)
'CalcDgpf D(0,2,1,0)',r,d(0,2,1,0)
9659 write(*,*)
'CalcDgpf Dij_err',r,dij_err
9660 write(*,*)
'CalcDgpf Dij_acc',r,dij_err/abs(d(0,0,0,0))
9662 write(*,*)
'CalcDgpf err',r,g,derr
9663 write(*,*)
'CalcDgpf acc',r,g,derr/abs(d(0,0,0,0))
9666 derr2 = max(derr,dij_err2(0:rmax))
9667 derr = max(derr,dij_err(0:rmax))
9672 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0)
then
9677 d(0,n1,n2,rg-n1-n2)=0d0
9684 d(1,n1,n2,rg-2-n1-n2)=0d0
9690 111
format(a22,2(
'(',g24.17,
',',g24.17,
') ':))
9691 call seterrflag_coli(-5)
9692 call errout_coli(
'CalcDgpf',
' exit rloop for D', &
9694 if (errorwriteflag)
then
9695 write(nerrout_coli,100)
' CalcDgpf: exit rloop for D ', &
9696 ' should not appear'
9697 write(nerrout_coli,111)
' CalcDgpf: p10 = ',p10
9698 write(nerrout_coli,111)
' CalcDgpf: p21 = ',p21
9699 write(nerrout_coli,111)
' CalcDgpf: p32 = ',p32
9700 write(nerrout_coli,111)
' CalcDgpf: p30 = ',p30
9701 write(nerrout_coli,111)
' CalcDgpf: p20 = ',p20
9702 write(nerrout_coli,111)
' CalcDgpf: p31 = ',p31
9703 write(nerrout_coli,111)
' CalcDgpf: m02 = ',m02
9704 write(nerrout_coli,111)
' CalcDgpf: m12 = ',m12
9705 write(nerrout_coli,111)
' CalcDgpf: m22 = ',m22
9706 write(nerrout_coli,111)
' CalcDgpf: m32 = ',m32
9711 if(maxval(derr-acc_req_d*abs(d(0,0,0,0))).le.0d0.and.r.ge.rmax)
then
9723 do n0=2,max(rmax,r/2)
9725 do nlt=r-2*n0-nl,0,-1
9726 nltt = r-2*n0-nl-nlt
9734 daux = shat(n0-1,inds(1),inds(2),inds(3),k) &
9735 -
f(k)*d(n0-1,inds(1),inds(2),inds(3)) &
9736 -
z(k,1)*d(n0-1,inds(1)+1,inds(2),inds(3)) &
9737 -
z(k,2)*d(n0-1,inds(1),inds(2)+1,inds(3)) &
9738 -
z(k,3)*d(n0-1,inds(1),inds(2),inds(3)+1)
9740 d(n0,inds0(1),inds0(2),inds0(3)) = daux/(2*(nl+1))
9757 write(*,*)
'CalcDgpf exp rmax+1',r,n0,n1,n2,n3, d(n0,n1,n2,n3)
9759 d(n0,n1,n2,n3) = (c_0(n0-1,n1,n2,n3) + 2*
mm02*d(n0-1,n1,n2,n3) &
9760 + 4*duv(n0,n1,n2,n3) &
9761 +
f(1)*d(n0-1,n1+1,n2,n3) +
f(2)*d(n0-1,n1,n2+1,n3) &
9762 +
f(3)*d(n0-1,n1,n2,n3+1)) / (2*(r-1))
9764 write(*,*)
'CalcDgpf dir rmax+1',r,n0,n1,n2,n3, d(n0,n1,n2,n3)
9774 write(*,*)
'CalcDgpf D(1,0,0,0) fin',d(1,0,0,0)
9775 write(*,*)
'CalcDgpf D(0,0,0,2) fin',d(0,0,0,2)
9777 write(*,*)
'CalcDgpf D(1,0,1,0) fin',d(1,0,1,0)
9778 write(*,*)
'CalcDgpf D(0,1,1,1) fin',d(0,1,1,1)
9782 write(*,*)
'CalcDgpf final err',derr
9783 write(*,*)
'CalcDgpf final acc',derr/abs(d(0,0,0,0))
9799 subroutine copydimp3(D,D_alt,Derr,Derr_alt,Derr1,Derr1_alt,Derr2,Derr2_alt,Drmethod,Drmethod_alt,rmax,r_alt)
9801 integer,
intent(in) :: rmax,r_alt
9802 double complex,
intent(inout) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
9803 double precision,
intent(inout) :: Derr(0:rmax),Derr1(0:rmax),Derr2(0:rmax)
9804 integer,
intent(inout) :: Drmethod(0:rmax)
9805 double complex,
intent(in) :: D_alt(0:r_alt,0:r_alt,0:r_alt,0:r_alt)
9806 double precision,
intent(in) :: Derr_alt(0:r_alt),Derr1_alt(0:r_alt),Derr2_alt(0:r_alt)
9807 integer,
intent(in) :: Drmethod_alt(0:r_alt)
9809 integer :: r,n1,n2,n0
9812 if (derr_alt(r).lt.derr(r))
then
9813 drmethod(r)=drmethod_alt(r)
9815 derr1(r)=derr1_alt(r)
9816 derr2(r)=derr2_alt(r)
9818 forall (n1=0:2*r-n0)
9819 forall (n2=0:r-2*n0-n1)
9820 d(n0,n1,n2,r-2*n0-n1-n2) = d_alt(n0,n1,n2,r-2*n0-n1-n2)
9824 forall (n0=1:(r+1)/2)
9825 forall (n1=0:r+1-2*n0)
9826 forall (n2=0:r+1-2*n0-n1)
9827 d(n0,n1,n2,r+1-2*n0-n1-n2) = d_alt(n0,n1,n2,r+1-2*n0-n1-n2)