21 #define ALWAYSPV ! default
24 #define Cutrloop ! default
28 #define PVEST2 ! default
51 double precision ::
q2max,
m2max,
m2scale,
adetz,
adetx,
fmax,
maxzadjf,
maxzadjfd,
maxxadj,
azadjff,
maxz
53 double complex ::
mx(0:2,0:2),
mxinv(0:2,0:2)
94 subroutine calcc(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,id,Cerr1,Cerr2,rbasic,acc_req_Cextra)
96 integer,
intent(in) :: rmax, id
97 integer,
optional,
intent(in) :: rbasic
98 double precision,
optional,
intent(in) :: acc_req_Cextra(0:rmax)
99 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
100 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
101 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
102 double precision,
intent(out) :: Cerr1(0:rmax),Cerr2(0:rmax)
103 double complex,
allocatable :: Caux(:,:,:), Cuvaux(:,:,:), fct(:)
104 double precision,
allocatable :: Cerr1aux(:),Cerr2aux(:)
105 double complex :: x(6)
106 integer :: rank,switch,cnt,n0,n1,n2,r,rb
107 logical :: nocalc,wrica
110 write(*,*)
'CalcC in ',rmax,id
113 if (use_cache_system)
then
114 if ((ncache.gt.0).and.(ncache.le.ncache_max))
then
126 allocate(fct(ncoefsg(rmax,3)+ncoefsg(rmax-1,3)+2*(rmax+2)))
127 call readcache(fct,ncoefsg(rmax,3)+ncoefsg(rmax-1,3)+2*(rmax+2),x,6,1,id,3,rank,nocalc,wrica)
129 allocate(fct(ncoefsg(rmax,3)+2*(rmax+2)))
130 call readcache(fct,ncoefsg(rmax,3)+2*(rmax+2),x,6,1,id,3,rank,nocalc,wrica)
135 if(
present(rbasic))
then
140 if(int(fct(1)).lt.rb)
then
152 c(n0,n1,n2) = fct(cnt)
161 cuv(n0,n1,n2) = fct(cnt)
166 cerr1(r) = real(fct(cnt))
168 cerr2(r) = real(fct(cnt))
176 if(rank.eq.rmax)
then
178 if(
present(rbasic))
then
179 call calccred(c,cuv,p10,p21,p20,m02,m12,m22,rank,id,cerr1,cerr2,rbasic+rank-rmax,acc_req_cextra)
181 call calccred(c,cuv,p10,p21,p20,m02,m12,m22,rank,id,cerr1,cerr2)
186 if(
present(rbasic))
then
197 fct(cnt) = c(n0,n1,n2)
205 fct(cnt) = cuv(n0,n1,n2)
215 call writecache(fct,ncoefsg(rank,3)+ncoefsg(rank-1,3)+2*(rank+2),id,3,rank)
217 call writecache(fct,ncoefsg(rank,3)+2*(rank+2),id,3,rank)
226 allocate(caux(0:rank,0:rank,0:rank))
227 allocate(cuvaux(0:rank,0:rank,0:rank))
228 allocate(cerr1aux(0:rank))
229 allocate(cerr2aux(0:rank))
231 if(
present(rbasic))
then
232 call calccred(caux,cuvaux,p10,p21,p20,m02,m12,m22,rank,id,cerr1aux,cerr2aux,rbasic+rank-rmax,acc_req_cextra)
234 call calccred(caux,cuvaux,p10,p21,p20,m02,m12,m22,rank,id,cerr1aux,cerr2aux)
241 allocate(fct(ncoefsg(rank,3)+ncoefsg(rank-1,3)+2*(rank+2)))
243 allocate(fct(ncoefsg(rank,3)+2*(rank+2)))
245 if(
present(rbasic))
then
246 fct(cnt) = rbasic+rank-rmax
256 fct(cnt) = caux(n0,n1,n2)
264 fct(cnt) = cuvaux(n0,n1,n2)
268 fct(cnt) = cerr1aux(r)
270 fct(cnt) = cerr2aux(r)
274 call writecache(fct,ncoefsg(rank,3)+ncoefsg(rank-1,3)+2*(rank+2),id,3,rank)
276 call writecache(fct,ncoefsg(rank,3)+2*(rank+2),id,3,rank)
281 c = caux(0:rmax,0:rmax,0:rmax)
282 cuv = cuvaux(0:rmax,0:rmax,0:rmax)
283 cerr1 = cerr1aux(0:rmax)
284 cerr2 = cerr2aux(0:rmax)
301 if(
present(rbasic))
then
302 call calccred(c,cuv,p10,p21,p20,m02,m12,m22,rmax,id,cerr1,cerr2,rbasic,acc_req_cextra)
304 call calccred(c,cuv,p10,p21,p20,m02,m12,m22,rmax,id,cerr1,cerr2)
321 subroutine calccred(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,id,Cerr1,Cerr2,rbasic,acc_req_Cextra)
325 integer,
intent(in) :: rmax,id
326 integer,
intent(in),
optional :: rbasic
329 double precision,
intent(in),
optional :: acc_req_Cextra(0:rmax)
330 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
331 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
332 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
333 double precision,
intent(out) ::Cerr1(0:rmax),Cerr2(0:rmax)
334 double complex :: C_alt(0:rmax,0:rmax,0:rmax)
335 double complex :: Cuv_alt(0:rmax,0:rmax,0:rmax)
336 double precision :: Cerr(0:rmax),Cerr_alt(0:rmax),Cerr1_alt(0:rmax),Cerr2_alt(0:rmax)
337 double precision :: C0est,Ctyp
339 double complex :: C0_coli
341 double complex :: C0_coli
343 double complex :: chdet
345 double complex :: elimminf2_coli
346 integer :: r,rid,n0,n1,n2,g,gy,gp,gr,gpf,i,rdef,iexp
347 logical :: use_pv,use_pv2,use_g,use_gy,use_gp,use_gr,use_gpf,use_pvs
349 integer :: r_alt,Crmethod(0:rmax),Crmethod_alt(0:rmax),CrCalc(0:rmax),CCalc
350 double precision :: acc_pv_alt, acc_pv2_alt, acc_Cr_alt
355 double precision :: err_pv(0:rmax),err_pv2(0:rmax),err_g(0:rmax),err_gy(0:rmax), &
356 err_gp(0:rmax),err_gr(0:rmax),err_gpf(0:rmax)
357 double precision :: err_pvs(0:rmax)
358 double precision :: h_pv,w_pv,v_pv,z_pv,h_pv2,w_pv2,v_pv2,z_pv2,hw_pv2
359 double precision :: h_pvs,w_pvs,v_pvs,z_pvs
360 double precision :: x_g,u_g,z_g,err_g_B(0:rmax),err_g_exp
361 double precision :: x_gy,y_gy,v_gy,v1_gy,b_gy,err_gy_B(0:rmax),err_gy_exp
362 double precision :: w_gp,v_gp,z_gp,err_gp_B(0:rmax),err_gp_exp
363 double precision :: x_gr,y_gr,y1_gr,a_gr,err_gr_B(0:rmax),err_gr_exp
364 double precision :: x_gpf,y_gpf,v_gpf,v1_gpf,b_gpf,err_gpf_B(0:rmax),err_gpf_exp
365 double precision :: err_B,err_C0,err_C(0:rmax),err_inf,err_req_Cr(0:rmax),acc_req_Cr(0:rmax),acc_C(0:rmax)
366 double precision :: checkest,norm,Cscale
367 logical :: lerr_C0,errorwriteflag
369 character(len=*),
parameter :: fmt1 =
"(A7,'dcmplx(',d25.18,' , ',d25.18,' )')"
370 character(len=*),
parameter :: fmt10 =
"(A17,'(',d25.18,' , ',d25.18,' )')"
371 #ifdef CritPointsCOLI
372 integer,
parameter :: MaxCritPointC=50
374 integer,
parameter :: MaxCritPointC=0
376 integer,
save :: CritPointCntC
378 data critpointcntc /0/
386 write(*,*)
'CalcCred in ',rmax,id,p10,p21,p20
389 write(*,*)
'CalcCred in ',rmax,id
392 if (
present(rbasic))
then
399 write(*,*)
'CalcCred rdef ',rdef
402 if (
present(acc_req_cextra))
then
403 acc_req_cr = acc_req_cextra
405 acc_req_cr = acc_req_c
409 write(*,*)
'CalcCred acc_req_Cr ',acc_req_cr
413 mm02 = elimminf2_coli(m02)
414 mm12 = elimminf2_coli(m12)
415 mm22 = elimminf2_coli(m22)
416 q10 = elimminf2_coli(p10)
417 q21 = elimminf2_coli(p21)
418 q20 = elimminf2_coli(p20)
431 maxz = maxval(abs(
z))
435 if (
detz.ne.0d0)
then
449 write(*,*)
'Z(1) ',
z(1,1),
z(1,2),
z(2,1),
z(2,2),
detz
450 write(*,*)
'Zadj(1) ',
zadj(1,1),
zadj(1,2),
zadj(2,1),
zadj(2,2),
detz
468 fmax = max(abs(
f(1)),abs(
f(2)))
481 mx(1:2,1:2) =
z(1:2,1:2)
485 if (
detx.ne.0d0.and.
maxz.ne.0d0)
then
539 write(*,*)
'maxZ ',
maxz,
z
543 write(*,*)
'CalcCred Zadjf ',
zadjf
544 write(*,*)
'CalcCred Xadj ',
xadj(1:2,1:2)
572 cscale = max(abs(p10),abs(p21),abs(p20),abs(m02), \
575 c0est = max(abs(c0_coli(p10,p21,p20,m02,m12,m22)),1d0/cscale)
579 if(cscale.ne.0d0)
then
597 write(*,*)
'CalcCred C0 = ',c0_coli(p10,p21,p20,m02,m12,m22)
598 if(
adetz.ne.0d0)
then
599 write(*,*)
'CalcCred C0est = ',c0est,1d0/sqrt(
adetz)
601 write(*,*)
'CalcCred C0est = ',c0est
605 err_inf = acc_inf*c0est
607 err_req_cr = acc_req_cr * c0est
619 if (
adetz.ne.0d0)
then
621 err_c0 = acc_def_c0*max( c0est, 1d0/sqrt(
adetz) )
623 err_c0 = acc_def_c0 * c0est
657 write(*,*)
'CalcCred w_pv',w_pv,v_pv,z_pv,h_pv
660 if (mod(rdef,2).eq.1)
then
661 err_pv(rdef) = max( w_pv**((rdef-1)/2) * v_pv * err_c0, &
662 max(w_pv**((rdef-1)/2),1d0) * z_pv * err_b )
665 write(*,*)
'CalcCred err_pv cont', w_pv**((rdef-1)/2)* v_pv* err_c0, &
666 w_pv**((rdef-1)/2) * z_pv * err_b, err_c0,err_b
667 write(*,*)
'CalcCred err_pv cont', w_pv**((rdef-1)/2),v_pv, err_c0
671 err_pv(rdef) = max( w_pv**(rdef/2) * err_c0, &
672 max(w_pv**(rdef/2-1) * v_pv, 1d0) * z_pv * err_b )
675 write(*,*)
'CalcCred w_pv', w_pv,err_c0,sqrt(w_pv)
676 write(*,*)
'CalcCred err_pv cont', w_pv**(rdef/2) * err_c0, &
677 w_pv**(rdef/2-1) * v_pv * z_pv * err_b, z_pv * err_b, err_c0,err_b
713 if (mod(rdef,2).eq.1)
then
720 err_pv2(rdef) = max( err_c0 * max(hw_pv2**rdef,hw_pv2*v_pv2**((rdef-1)/2) ), &
721 err_b * z_pv2 * max(w_pv2*hw_pv2**(rdef),hw_pv2, &
722 w_pv2*hw_pv2*v_pv2**((rdef-1)/2), &
723 hw_pv2*v_pv2**((rdef-1)/2),w_pv2*hw_pv2, &
724 v_pv2**((rdef+1)/2),v_pv2 ) )
737 err_pv2(rdef) = max( err_c0 * max(hw_pv2**rdef,v_pv2**(rdef/2)), &
738 err_b * z_pv2 * max(w_pv2*hw_pv2**(rdef),hw_pv2, &
739 w_pv2*v_pv2**(rdef/2), w_pv2*hw_pv2, &
740 v_pv2**(rdef/2),v_pv2) )
746 err_pv(rdef) = err_pv(rdef)/impest_c
747 err_pv2(rdef) = err_pv2(rdef)/impest_c
759 write(*,*)
'CalcCred: err_pv',err_pv(rdef),err_pv2(rdef),err_req_cr(rdef)
760 write(*,*)
'CalcCred: acc_pv',err_pv(rdef)/c0est,err_pv2(rdef)/c0est,acc_req_c
768 if(use_pv.or.use_pv2)
then
770 if (min(err_pv(rdef),err_pv2(rdef)).le.err_req_cr(rdef))
then
772 if (err_pv(rdef).le.err_pv2(rdef))
then
775 write(*,*)
'CalcCred: call Cpv 1 ',rmax,id,err_pv(rdef)
779 call calccpv1(c,cuv,p10,p21,p20,m02,m12,m22,rmax,id,cerr1,cerr2)
786 crcalc(0:rmax)=crcalc(0:rmax)+1
791 checkest=cerr(rdef)/err_pv(rdef)
792 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
793 write(*,*)
'CalcCred: estimate err_pv imprecise',err_pv(rdef),cerr(rdef)
801 write(*,*)
'CalcCred: call Cpv2 1',rdef,id,err_pv2(rdef)
805 call calccpv2(c,cuv,p10,p21,p20,m02,m12,m22,rmax,id,cerr1,cerr2)
812 crcalc(0:rmax)=crcalc(0:rmax)+2
817 checkest=cerr(rdef)/err_pv2(rdef)
818 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
819 write(*,*)
'CalcCred: estimate err_pv2 imprecise',err_pv2(rdef),cerr(rdef)
830 err_c0 = acc_def_c0*max( abs(c(0,0,0)), 1d0/sqrt(
adetz) )
836 ctyp = max(abs(c(0,0,0)),abs(c(0,1,0)),abs(c(0,0,1)))
840 if(ctyp.eq.0d0) ctyp = c0est
841 err_req_cr = acc_req_cr * ctyp
845 write(*,*)
'CalcCred C0est after PV=',abs(c(0,0,0)),ctyp
846 write(*,*)
'CalcCred Cerr after PV =',cerr
847 write(*,*)
'CalcCred Cacc after PV =',cerr/ctyp
848 write(*,*)
'CalcCred err_req =',err_req_cr
849 write(*,*)
'CalcCred Cerr1-req=',cerr1-err_req_cr
850 write(*,*)
'CalcCred max(Cerr-req)',maxval(cerr-err_req_cr)
851 write(*,*)
'CalcCred max(Cerr1-req)',maxval(cerr1-err_req_cr)
855 if (maxval(cerr1-err_req_cr).lt.0)
then
893 err_g_b(rdef) = err_b * u_g**rdef * z_g
894 err_g_exp = u_g**(rdef-1) * ctyp
907 write(*,*)
'CalcCred: after Gram pars',use_g,
fac_g,x_g,u_g,z_g,err_g_b(rdef),err_g_exp
909 write(*,*)
'CalcCred: after Gram pars',use_g,err_g_exp
920 v1_gy = max(1d0,v_gy)
921 fac_gy = max(x_gy,y_gy)*v1_gy
932 err_gy_b(rdef) = err_b * b_gy*v1_gy
933 err_gy_exp = 1d0 * ctyp
946 write(*,*)
'CalcCred: after GramCay pars',use_gy,
fac_gy,x_gy,y_gy,v_gy,b_gy,err_gy_b(rdef),err_gy_exp
948 write(*,*)
'CalcCred: after GramCay pars',use_gy,err_gy_exp
969 err_gp_b(rdef) = err_b * z_gp*v_gp**rdef
970 err_gp_exp = v_gp**(rdef-1) * ctyp
983 write(*,*)
'CalcCred: after Mom pars',use_gp,
fac_gp,w_gp,v_gp,z_gp,err_gp_b(rdef),err_gp_exp
985 write(*,*)
'CalcCred: after Mom pars',use_gp,err_gp_exp
994 y1_gr = max(1d0,y_gr)
999 if (
fac_gr.ge.1.or.2*rmax.gt.rmax_b)
then
1001 err_gr_exp = err_inf
1004 err_gr_b(rdef) = err_b * a_gr
1005 err_gr_exp = y1_gr * ctyp
1011 err_gr_exp = err_inf
1018 write(*,*)
'CalcCred: after revGram pars',use_gr,
fac_gr,x_gr,y_gr,y1_gr,a_gr,err_gr_b(rdef),err_gr_exp
1020 write(*,*)
'CalcCred: after revGram pars',use_gr,err_gr_exp
1027 if (abs(m02).gt.
m2scale*dprec_cll)
then
1028 x_gpf =
fmax/abs(m02)
1029 y_gpf =
maxz/abs(m02)
1031 v1_gpf = max(1d0,v_gpf)
1032 fac_gpf = max(x_gpf,y_gpf)*v1_gpf
1037 err_gpf_exp = err_inf
1041 b_gpf = 1d0/abs(m02)
1042 err_gpf_b(rdef) = err_b * b_gpf*v1_gpf
1043 err_gpf_exp = 1d0 * ctyp
1049 err_gpf_exp = err_inf
1056 write(*,*)
'CalcCred: after pf pars',use_gpf,
fac_gpf,x_gpf,y_gpf,v_gpf,b_gpf,err_gpf_b(rdef),err_gpf_exp
1058 write(*,*)
'CalcCred: after pf pars',use_gpf,err_gpf_exp
1064 if(use_pv.or.use_pv2.or.use_g.or.use_gy.or.use_gp.or.use_gr.or.use_gpf.eqv..false.)
then
1065 call seterrflag_coli(-6)
1066 call errout_coli(
'CalcCred',
' no reduction method works', &
1069 if (errorwriteflag)
then
1070 write(nerrout_coli,fmt10)
' CalcCred: p10 = ',p10
1071 write(nerrout_coli,fmt10)
' CalcCred: p21 = ',p21
1072 write(nerrout_coli,fmt10)
' CalcCred: p20 = ',p20
1073 write(nerrout_coli,fmt10)
' CalcCred: m02 = ',m02
1074 write(nerrout_coli,fmt10)
' CalcCred: m12 = ',m12
1075 write(nerrout_coli,fmt10)
' CalcCred: m22 = ',m22
1083 write(*,*)
'CalcCred: exit'
1102 if (err_g_exp.gt.err_g_b(rdef))
then
1104 err_g_exp = err_g_exp*
fac_g
1105 err_g(rdef) = max(err_g_exp,err_g_b(rdef))
1106 if(err_g(rdef).lt.err_req_cr(rdef))
then
1109 g = min(max(g+2,2*g),rmax_c-rmax)
1114 write(*,*)
'CalcCred i g',i,g,err_g_exp,err_g_b(rdef),err_g(rdef)
1120 if (mod(i,2).eq.1)
then
1122 if (err_gy_exp.gt.err_gy_b(rdef))
then
1124 err_gy_exp = err_gy_exp*
fac_gy
1125 err_gy(rdef) = max(err_gy_exp, err_gy_b(rdef))
1126 if(err_gy(rdef).lt.err_req_cr(rdef))
then
1129 gy = min(max(gy+4,2*gy),(rmax_c-rmax)/2)
1134 write(*,*)
'CalcCred i gy',i,gy,err_gy_exp,err_gy_b(rdef),err_gy(rdef)
1142 if (err_gp_exp.gt.err_gp_b(rdef))
then
1144 err_gp_exp = err_gp_exp*
fac_gp
1145 err_gp(rdef) = max(err_gp_exp,err_gp_b(rdef))
1146 if(err_gp(rdef).lt.err_req_cr(rdef))
then
1149 gp = min(max(gp+2,2*gp),rmax_c-rmax)
1154 write(*,*)
'CalcCred i gp',i,gp,err_gp_exp,err_gp_b(rdef),err_gp(rdef)
1160 if (mod(i,2).eq.1)
then
1165 write(*,*)
'CalcCred: it gr',use_gr,err_gr_exp,err_gr_b(rdef),err_gr(rdef), &
1169 if (err_gr_exp.gt.err_gr_b(rdef))
then
1171 err_gr_exp = err_gr_exp*
fac_gr
1172 err_gr(rdef) = max(err_gr_exp, err_gr_b(rdef))
1173 if(err_gr(rdef).lt.err_req_cr(rdef))
then
1178 gr = min(max(gr+4,2*gr),rmax_c-rmax,max(0,(rmax_b-2*rmax)/2))
1184 write(*,*)
'CalcCred: it gr',i,gr, err_gr_exp,err_gr_b(rdef) ,err_gr(rdef)
1190 if (mod(i,2).eq.1)
then
1192 if (err_gpf_exp.gt.err_gpf_b(rdef))
then
1194 err_gpf_exp = err_gpf_exp*
fac_gpf
1195 err_gpf(rdef) = max(err_gpf_exp, err_gpf_b(rdef))
1196 if(err_gpf(rdef).lt.err_req_cr(rdef))
then
1199 gpf = min(max(gpf+4,2*gpf),(rmax_c-rmax)/2)
1204 write(*,*)
'CalcCred i gpf',i,gpf,err_gpf_exp,err_gpf_b(rdef),err_gpf(rdef),err_req_cr(rdef)
1214 err_g(rdef) = err_g(rdef)/impest_c
1215 err_gy(rdef) = err_gy(rdef)/impest_c
1216 err_gp(rdef) = err_gp(rdef)/impest_c
1217 err_gr(rdef) = err_gr(rdef)/impest_c
1218 err_gpf(rdef)= err_gpf(rdef)/impest_c
1221 write(*,*)
'iexp=',iexp
1223 write(*,*)
'errexp=',err_g_exp,err_gy_exp,err_gp_exp,err_gr_exp,err_gpf_exp,err_req_cr(rdef)
1224 write(*,*)
'errexptot=',i
1225 write(*,*)
'g: errexptot =',g,err_g(rdef)
1226 write(*,*)
'gy: errexptot =',gy,err_gy(rdef)
1227 write(*,*)
'gp: errexptot =',gp,err_gp(rdef)
1228 write(*,*)
'gr: errexptot =',gr,err_gr(rdef)
1229 write(*,*)
'gpf: errexptot=',gpf,err_gpf(rdef)
1230 write(*,*)
'errexptot=',i,g,err_g(rdef),gy,err_gy(rdef),gp,err_gp(rdef),gr,err_gr(rdef),gpf,err_gpf(rdef)
1231 write(*,*)
'accexptot=',i,g,err_g(rdef)/ctyp,gy,err_gy(rdef)/ctyp,gp,err_gp(rdef)/ctyp, &
1232 gr,err_gr(rdef)/ctyp,gpf,err_gpf(rdef)/ctyp
1247 call calccgn(c_alt,cuv,p10,p21,p20,m02,m12,m22,rmax,g,g,id,cerr1_alt,acc_req_cr,cerr2_alt)
1249 cerr_alt = cerr2_alt
1251 cerr_alt = cerr1_alt
1254 crcalc(0:rmax)=crcalc(0:rmax)+4
1256 crmethod_alt(0:rmax)=4
1259 checkest=cerr_alt(rdef)/err_g(rdef)
1260 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
1261 write(*,*)
'CalcCred: estimate err_g imprecise ',err_g(rdef),cerr_alt(rdef)
1267 call copycimp3(c,c_alt,cerr,cerr_alt,cerr1,cerr1_alt,cerr2,cerr2_alt,crmethod,crmethod_alt,rmax,rmax)
1270 write(*,*)
'CalcCred Cerr after exp =',cerr
1271 write(*,*)
'CalcCred Cacc=',cerr/ctyp
1272 write(*,*)
'CalcCred method=',crmethod
1277 call calccg(c_alt,cuv,p10,p21,p20,m02,m12,m22,rmax,g,g,id,cerr1_alt,acc_req_cr,cerr2_alt)
1279 cerr_alt = cerr2_alt
1281 cerr_alt = cerr1_alt
1284 crcalc(0:rmax)=crcalc(0:rmax)+4
1286 crmethod_alt(0:rmax)=4
1289 checkest=cerr_alt(rdef)/err_g(rdef)
1290 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
1291 write(*,*)
'CalcCred: estimate err_g imprecise ',err_g(rdef),cerr_alt(rdef)
1297 call copycimp3(c,c_alt,cerr,cerr_alt,cerr1,cerr1_alt,cerr2,cerr2_alt,crmethod,crmethod_alt,rmax,rmax)
1300 write(*,*)
'CalcCred Cerr after exp =',cerr
1301 write(*,*)
'CalcCred Cacc=',cerr/ctyp
1302 write(*,*)
'CalcCred method=',crmethod
1307 call calccgy(c_alt,cuv,p10,p21,p20,m02,m12,m22,rmax,gy,gy,id,cerr1_alt,acc_req_cr,cerr2_alt)
1309 cerr_alt = cerr2_alt
1311 cerr_alt = cerr1_alt
1314 crcalc(0:rmax)=crcalc(0:rmax)+8
1316 crmethod_alt(0:rmax)=8
1319 checkest=cerr_alt(rdef)/err_gy(rdef)
1320 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
1321 write(*,*)
'CalcCred: estimate err_gy imprecise',err_gy(rdef),cerr_alt(rdef),checkest
1327 call copycimp3(c,c_alt,cerr,cerr_alt,cerr1,cerr1_alt,cerr2,cerr2_alt,crmethod,crmethod_alt,rmax,rmax)
1330 write(*,*)
'CalcCred Cerr after exp =',cerr
1331 write(*,*)
'CalcCred Cacc=',cerr/ctyp
1332 write(*,*)
'CalcCred method=',crmethod
1336 call calccgp(c_alt,cuv,p10,p21,p20,m02,m12,m22,rmax,gp,gp,id,cerr1_alt,acc_req_cr,cerr2_alt)
1338 cerr_alt = cerr2_alt
1340 cerr_alt = cerr1_alt
1343 crcalc(0:rmax)=crcalc(0:rmax)+16
1345 crmethod_alt(0:rmax)=16
1348 checkest=cerr_alt(rdef)/err_gp(rdef)
1349 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
1350 write(*,*)
'CalcCred: estimate err_gp imprecise',err_gp(rdef),cerr_alt(rdef)
1356 call copycimp3(c,c_alt,cerr,cerr_alt,cerr1,cerr1_alt,cerr2,cerr2_alt,crmethod,crmethod_alt,rmax,rmax)
1359 write(*,*)
'CalcCred Cerr after exp =',cerr
1360 write(*,*)
'CalcCred Cacc=',cerr/ctyp
1361 write(*,*)
'CalcCred method=',crmethod
1365 call calccgr(c_alt,cuv,p10,p21,p20,m02,m12,m22,rmax,gr,gr,id,cerr1_alt,acc_req_cr,cerr2_alt)
1367 cerr_alt = cerr2_alt
1369 cerr_alt = cerr1_alt
1372 crcalc(0:rmax)=crcalc(0:rmax)+32
1374 crmethod_alt(0:rmax)=32
1377 checkest=cerr_alt(rdef)/err_gr(rdef)
1378 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
1379 write(*,*)
'CalcCred: estimate err_gr imprecise',err_gr(rdef),cerr_alt(rdef)
1385 call copycimp3(c,c_alt,cerr,cerr_alt,cerr1,cerr1_alt,cerr2,cerr2_alt,crmethod,crmethod_alt,rmax,rmax)
1388 write(*,*)
'CalcCred Cerr after exp =',cerr
1389 write(*,*)
'CalcCred Cacc=',cerr/ctyp
1390 write(*,*)
'CalcCred method=',crmethod
1395 call calccgpf(c_alt,cuv,p10,p21,p20,m02,m12,m22,rmax,gpf,gpf,id,cerr1_alt,acc_req_cr,cerr2_alt)
1397 cerr_alt = cerr2_alt
1399 cerr_alt = cerr1_alt
1402 crcalc(0:rmax)=crcalc(0:rmax)+8
1404 crmethod_alt(0:rmax)=8
1407 checkest=cerr_alt(rdef)/err_gpf(rdef)
1408 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
1409 write(*,*)
'CalcCred: estimate err_gpf imprecise',err_gpf(rdef),cerr_alt(rdef),checkest
1415 call copycimp3(c,c_alt,cerr,cerr_alt,cerr1,cerr1_alt,cerr2,cerr2_alt,crmethod,crmethod_alt,rmax,rmax)
1418 write(*,*)
'CalcCred Cerr after exp =',cerr
1419 write(*,*)
'CalcCred Cacc=',cerr/ctyp
1420 write(*,*)
'CalcCred method=',crmethod
1428 if(.not.lerr_c0.and.iexp.ne.0)
then
1430 err_c0 = acc_def_c0*max( abs(c(0,0,0)), 1d0/sqrt(
adetz) )
1439 ctyp = max(abs(c(0,0,0)),abs(c(0,1,0)),abs(c(0,0,1)))
1441 ctyp = abs(c(0,0,0))
1443 err_req_cr = acc_req_cr * ctyp
1445 if (maxval(cerr1-err_req_cr).lt.0)
then
1452 write(*,*)
'CalcCred no optimal method bef shift',ctyp,acc_req_cr
1453 write(*,*)
'err_req=',err_req_cr(rdef),rdef
1454 write(*,*)
'err_est=',err_pv(rdef),err_pv2(rdef),err_g(rdef) &
1455 ,err_gy(rdef),err_gp(rdef),err_gr(rdef),err_gpf(rdef)
1467 write(*,*)
'CalcCred try pv shift i = ',i
1516 #ifdef Credtestshift
1558 #ifdef Credtestshift
1604 write(*,*)
'CalcCred w_pvs',w_pvs,v_pvs,z_pv,h_pv
1607 if (mod(rdef,2).eq.1)
then
1608 err_pvs(rdef) = max( w_pvs**((rdef-1)/2) * v_pvs * err_c0, &
1609 max(w_pvs**((rdef-1)/2),1d0) * z_pv * err_b )
1612 write(*,*)
'CalcCred err_pvs', w_pvs**((rdef-1)/2)* v_pvs* err_c0, &
1613 w_pvs**((rdef-1)/2) * z_pv * err_b, err_c0,err_b
1614 write(*,*)
'CalcCred err_pvs', w_pvs**((rdef-1)/2),v_pvs, err_c0
1618 err_pvs(rdef) = max( w_pvs**(rdef/2) * err_c0, &
1619 max(w_pvs**(rdef/2-1) * v_pvs, 1d0) * z_pv * err_b )
1622 write(*,*)
'CalcCred w_pvs', w_pvs,err_c0,sqrt(w_pvs)
1625 write(*,*)
'CalcCred err_pvs', w_pvs**(rdef/2) * err_c0, &
1626 w_pvs**(rdef/2-1) * v_pvs * z_pv * err_b, z_pv * err_b, err_c0,err_b
1634 write(*,*)
'CalcCred use_pvs',use_pvs,err_pvs(rdef).lt.err_pv(rdef),i
1635 write(*,*)
'CalcCred err_pvs',err_pvs(rdef),err_pv(rdef),i
1638 if(use_pvs.and.err_pvs(rdef).lt. min(err_pv(rdef),err_pv2(rdef),err_g(rdef) &
1639 ,err_gy(rdef),err_gp(rdef),err_gr(rdef),err_gpf(rdef)) )
then
1642 write(*,*)
'CalcCred: call Cpvs 1 ',rmax,id,err_pvs(rdef)
1647 call calccpvshift(c_alt,cuv,p10,p20,p21,m12,m02,m22,rmax,id,cerr1_alt,cerr2_alt)
1654 c_alt(0:n0,n1,n2) = -c_alt(0:n0,n1-1,n2)-c_alt(0:n0,n1,n2)-c_alt(0:n0,n1-1,n2+1)
1660 elseif (i.eq.2)
then
1661 call calccpvshift(c_alt,cuv,p21,p10,p20,m22,m12,m02,rmax,id,cerr1_alt,cerr2_alt)
1667 c_alt(0:n0,n1,n2) = -c_alt(0:n0,n1,n2-1)-c_alt(0:n0,n1,n2)-c_alt(0:n0,n1+1,n2-1)
1676 cerr_alt = cerr2_alt
1678 cerr_alt = cerr1_alt
1683 crmethod_alt(0:rmax)=1
1684 if (cerr_alt(rmax).lt.cerr(rmax))
then
1689 checkest=cerr_alt(rdef)/err_pvs(rdef)
1690 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
1691 write(*,*)
'CalcCred: estimate err_pvs imprecise ',err_pvs(rdef),cerr_alt(rdef)
1697 call copycimp3(c,c_alt,cerr,cerr_alt,cerr1,cerr1_alt,cerr2,cerr2_alt,crmethod,crmethod_alt,rmax,rmax)
1700 ctyp = max(abs(c(0,0,0)),abs(c(0,1,0)),abs(c(0,0,1)))
1702 ctyp = abs(c(0,0,0))
1704 err_req_cr = acc_req_cr * ctyp
1707 write(*,*)
'CalcCred C0est after PVS=',abs(c(0,0,0)),ctyp
1708 write(*,*)
'CalcCred Cerr after PVS =',cerr
1709 write(*,*)
'CalcCred Cacc after PVS =',cerr/ctyp
1710 write(*,*)
'CalcCred err_req =',err_req_cr
1711 write(*,*)
'CalcCred Cerr1-req=',cerr1-err_req_cr
1712 write(*,*)
'CalcCred max(Cerr-req)',maxval(cerr-err_req_cr)
1713 write(*,*)
'CalcCred max(Cerr1-req)',maxval(cerr1-err_req_cr)
1717 if (maxval(cerr1-err_req_cr).lt.0)
then
1727 write(*,*)
'CalcCred no optimal method aft shift',ctyp,acc_req_cr
1728 write(*,*)
'err_req=',err_req_cr(rdef),rdef
1729 write(*,*)
'err_est=',err_pv(rdef),err_pv2(rdef),err_g(rdef) &
1730 ,err_gy(rdef),err_gp(rdef),err_gr(rdef),err_gpf(rdef)
1743 if(use_pv.and.mod(crcalc(r),2).ne.1)
then
1749 if (mod(r,2).eq.1)
then
1750 err_pv(r) = max( w_pv**((r-1)/2) * v_pv * err_c0, &
1751 max(w_pv**((r-1)/2),1d0) * z_pv * err_b )
1756 else if (r.ne.0)
then
1757 err_pv(r) = max( w_pv**(r/2) * err_c0, &
1758 max(w_pv**(r/2-1) * v_pv , 1d0) * z_pv * err_b )
1769 err_pv(r) = err_pv(r)/impest_c
1772 if (use_pv2.and.mod(crcalc(r),4)-mod(crcalc(r),2).ne.2)
then
1775 if (mod(r,2).eq.1)
then
1781 err_pv2(r) = max( err_c0 * max(hw_pv2**r,hw_pv2*v_pv2**((r-1)/2) ), &
1782 err_b * z_pv2 * max(w_pv2*hw_pv2**(r),hw_pv2, &
1783 hw_pv2*v_pv2**((r-1)/2),w_pv2*hw_pv2, &
1784 w_pv2*hw_pv2*v_pv2**((r-1)/2), &
1785 v_pv2**((r+1)/2),v_pv2) )
1799 err_pv2(r) = max( err_c0 * max(hw_pv2**r,v_pv2**(r/2)), &
1800 err_b * z_pv2 * max(w_pv2*hw_pv2**(r),hw_pv2, &
1801 hw_pv2*v_pv2**(r/2), w_pv2*hw_pv2, &
1802 v_pv2**(r/2),v_pv2) )
1805 err_pv2(r) = err_inf
1808 err_pv2(r) = err_pv2(r)/impest_c
1811 if (use_g.and.mod(crcalc(r),8)-mod(crcalc(r),4).ne.4)
then
1813 err_g_b(r) = err_b * u_g**r * z_g
1814 err_g_exp = u_g**(r-1) * ctyp
1820 err_g_exp = err_g_exp*
fac_g
1821 err_g(r) = max(err_g_exp,err_g_b(r))
1825 if (err_g_exp.lt.err_g_b(r).or.err_g(r).lt.err_req_cr(r))
exit
1828 g = min(max(g+2,2*g),rmax_c-r)
1830 err_g(r) = err_g(r)/impest_c
1833 if (use_gy.and.mod(crcalc(r),16)-mod(crcalc(r),8).ne.8)
then
1835 err_gy_b(r) = err_b * b_gy*v1_gy
1836 err_gy_exp = 1d0 * ctyp
1841 if (mod(i,2).eq.1)
then
1843 err_gy_exp = err_gy_exp*
fac_gy
1844 err_gy(r) = max(err_gy_exp, err_gy_b(r))
1845 if (err_gy_exp.lt.err_gy_b(r).or.err_gy(r).lt.err_req_cr(r))
exit
1849 gy = min(max(gy+4,2*gy),(rmax_c-r)/2)
1851 err_gy(r) = err_gy(r)/impest_c
1854 if (use_gp.and.mod(crcalc(r),32)-mod(crcalc(r),16).ne.16)
then
1856 err_gp_b(r) = err_b * z_gp*v_gp**r
1857 err_gp_exp = v_gp**(r-1) * ctyp
1862 err_gp_exp = err_gp_exp*
fac_gp
1863 err_gp(r) = max(err_gp_exp,err_gp_b(r))
1864 if (err_gp_exp.lt.err_gp_b(r).or.err_gp(r).lt.err_req_cr(r))
exit
1867 gp = min(max(gp+2,2*gp),rmax_c-r)
1869 err_gp(r) = err_gp(r)/impest_c
1872 if (mod(crcalc(r),64)-mod(crcalc(r),32).ne.32.and.use_gr)
then
1874 err_gr_b(r) = err_b * a_gr
1875 err_gr_exp = y1_gr * ctyp
1880 if (mod(i,2).eq.1)
then
1882 err_gr_exp = err_gr_exp*
fac_gr
1883 err_gr(r) = max(err_gr_exp,err_gr_b(r))
1886 write(*,*)
'CalcCgr err_gr',i,gr,err_gr_exp,err_gr_b(r),err_gr(r),err_req_cr(r)
1889 if (err_gr_exp.lt.err_gr_b(r).or.err_gr(r).lt.err_req_cr(r))
exit
1895 gr = min(max(gr+2,2*gr),rmax_c-r,max(0,(rmax_b-2*r)/2))
1897 err_gr(r) = err_gr(r)/impest_c
1901 if (use_gpf.and.mod(crcalc(r),128)-mod(crcalc(r),64).ne.64)
then
1903 err_gpf_b(r) = err_b * b_gpf*v1_gpf
1904 err_gpf_exp = 1d0 * ctyp
1909 if (mod(i,2).eq.1)
then
1911 err_gpf_exp = err_gpf_exp*
fac_gpf
1912 err_gpf(r) = max(err_gpf_exp, err_gpf_b(r))
1913 if (err_gpf_exp.lt.err_gpf_b(r).or.err_gpf(r).lt.err_req_cr(r))
exit
1917 gpf = min(max(gpf+4,2*gpf),(rmax_c-r)/2)
1919 err_gpf(r) = err_gpf(r)/impest_c
1924 write(*,*)
'CalcCred: bef final loop expansion depth',r,g,gy,gp,gr,gpf
1925 write(*,*)
'CalcCred: bef final loop err methods',r,err_pv(r),err_pv2(r) &
1926 ,err_g(r),err_gy(r),err_gp(r),err_gr(r),err_gp(r),err_gpf(r)
1927 write(*,*)
'CalcCred: bef final loop acc methods',r,err_pv(r)/ctyp,err_pv2(r)/ctyp, &
1928 err_g(r)/ctyp,err_gy(r)/ctyp,err_gp(r)/ctyp,err_gr(r)/ctyp,err_gpf(r)/ctyp
1929 write(*,*)
'CalcCred: bef final loop',r,crcalc(r),crmethod(r)
1934 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)) &
1935 .and.min(err_pv(r),err_pv2(r)).lt.err_inf)
then
1937 if (use_pv.and.err_pv(r).le.err_pv2(r).and.mod(crcalc(r),2).ne.1)
then
1951 write(*,*)
'CalcCred: call Cpv 2',r,id
1954 call calccpv1(c_alt,cuv,p10,p21,p20,m02,m12,m22,r,id,cerr1_alt,cerr2_alt)
1956 call calccpv1(c_alt(0:r,0:r,0:r),cuv_alt(0:r,0:r,0:r), &
1957 p10,p21,p20,m02,m12,m22,r,id,cerr1_alt(0:r),cerr2_alt(0:r))
1960 cerr_alt = cerr2_alt
1962 cerr_alt = cerr1_alt
1965 crcalc(0:r)=crcalc(0:r)+1
1968 checkest=cerr_alt(r)/(err_pv(r)*abs(c_alt(0,0,0)))
1971 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
1972 write(*,*)
'CalcCred: estimate acc_pv imprecise',err_pv(r),cerr_alt(r)/abs(c_alt(0,0,0))
1976 err_pv(0:r)=cerr_alt(0:r)
1978 call copycimp3(c,c_alt(0:r,0:r,0:r),cerr,cerr_alt(0:r),cerr1,cerr1_alt(0:r), &
1979 cerr2,cerr2_alt(0:r),crmethod,crmethod_alt(0:r),rmax,r)
1982 ctyp = max(abs(c(0,0,0)),abs(c(0,1,0)),abs(c(0,0,1)))
1984 ctyp = abs(c(0,0,0))
1986 err_req_cr = acc_req_cr * ctyp
1989 write(*,*)
'CalcCred: after pv 2nd try Cmethod',crmethod
1990 write(*,*)
'CalcCred: after pv 2nd try Cerr(r)',cerr
1991 write(*,*)
'CalcCred: after pv 2nd try Cacc(r)',cerr/ctyp
1994 if(checkest.gt.impest_c.and.mode_coli.lt.1)
goto 100
1996 elseif (use_pv2.and.err_pv2(r).le.err_pv(r).and.mod(crcalc(r),4)-mod(crcalc(r),2).ne.2)
then
2010 call calccpv2(c_alt,cuv,p10,p21,p20,m02,m12,m22,r,id,cerr1_alt,cerr2_alt)
2012 call calccpv2(c_alt(0:r,0:r,0:r),cuv_alt(0:r,0:r,0:r),p10,p21,p20,m02,m12,m22,r,id,cerr1_alt(0:r),cerr2_alt(0:r))
2015 cerr_alt = cerr2_alt
2017 cerr_alt = cerr1_alt
2020 crcalc(0:r)=crcalc(0:r)+2
2023 checkest=cerr_alt(r)/(err_pv(r)*abs(c_alt(0,0,0)))
2026 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
2027 write(*,*)
'CalcCred: estimate err_pv2 imprecise',err_pv2(r),cerr_alt(r)
2031 err_pv2(0:r)=cerr_alt(0:r)
2033 call copycimp3(c,c_alt(0:r,0:r,0:r),cerr,cerr_alt(0:r),cerr1,cerr1_alt(0:r), &
2034 cerr2,cerr2_alt(0:r),crmethod,crmethod_alt(0:r),rmax,r)
2037 ctyp = max(abs(c(0,0,0)),abs(c(0,1,0)),abs(c(0,0,1)))
2039 ctyp = abs(c(0,0,0))
2041 err_req_cr = acc_req_cr * ctyp
2044 write(*,*)
'CalcCred: after pv 2nd try Cmethod',crmethod
2045 write(*,*)
'CalcCred: after pv 2nd try Cerr(r)',cerr
2046 write(*,*)
'CalcCred: after pv 2nd try Cacc(r)',cerr/ctyp
2049 if(checkest.gt.impest_c.and.mode_coli.lt.1)
goto 100
2056 write(*,*)
'CalcCred: explore exps once more'
2059 if (use_g.and.err_g(r).le.min(err_gy(r),err_gp(r),err_gr(r),err_gpf(r)) &
2060 .and.mod(crcalc(r),8)-mod(crcalc(r),4).ne.4)
then
2074 call calccg(c_alt,cuv,p10,p21,p20,m02,m12,m22,r,g,rmax_c,id,cerr1_alt,acc_req_cr,cerr2_alt)
2076 call calccg(c_alt(0:r,0:r,0:r),cuv_alt(0:r,0:r,0:r), &
2077 p10,p21,p20,m02,m12,m22,r,g,rmax_c,id,cerr1_alt(0:r),acc_req_cr(0:r),cerr2_alt(0:r))
2080 cerr_alt = cerr2_alt
2082 cerr_alt = cerr1_alt
2085 crcalc(0:r)=crcalc(0:r)+4
2088 checkest=cerr_alt(r)/err_g(r)
2091 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
2092 write(*,*)
'CalcCred: estimate err_g imprecise ',err_g(r),cerr_alt(r)
2096 err_g(0:r)=cerr_alt(0:r)
2098 call copycimp3(c,c_alt(0:r,0:r,0:r),cerr,cerr_alt(0:r),cerr1,cerr1_alt(0:r), &
2099 cerr2,cerr2_alt(0:r),crmethod,crmethod_alt(0:r),rmax,r)
2102 ctyp = max(abs(c(0,0,0)),abs(c(0,1,0)),abs(c(0,0,1)))
2104 ctyp = abs(c(0,0,0))
2106 err_req_cr = acc_req_cr * ctyp
2109 write(*,*)
'CalcCred: after exps 2nd try Cmethod',crmethod
2110 write(*,*)
'CalcCred: after exps 2nd try Cerr(r)',cerr
2111 write(*,*)
'CalcCred: after exps 2nd try Cacc(r)',cerr/ctyp
2114 if(checkest.gt.impest_c.and.mode_coli.lt.1)
goto 100
2116 else if (use_gy.and.err_gy(r).le.min(err_g(r),err_gp(r),err_gr(r),err_gpf(r)) &
2117 .and.mod(crcalc(r),16)-mod(crcalc(r),8).ne.8)
then
2130 call calccgy(c_alt,cuv,p10,p21,p20,m02,m12,m22,r,gy,(rmax_c)/2,id,cerr1_alt,acc_req_cr,cerr2_alt)
2132 call calccgy(c_alt(0:r,0:r,0:r),cuv_alt(0:r,0:r,0:r), &
2133 p10,p21,p20,m02,m12,m22,r,gy,(rmax_c)/2,id,cerr1_alt(0:r),acc_req_cr(0:r),cerr2_alt(0:r))
2136 cerr_alt = cerr2_alt
2138 cerr_alt = cerr1_alt
2141 crcalc(0:r)=crcalc(0:r)+8
2144 checkest=cerr_alt(r)/err_gy(r)
2147 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
2148 write(*,*)
'CalcCred: estimate err_gy imprecise',err_gy(r),cerr_alt(r)
2151 err_gy(0:r)=cerr_alt(0:r)
2153 call copycimp3(c,c_alt(0:r,0:r,0:r),cerr,cerr_alt(0:r),cerr1,cerr1_alt(0:r), &
2154 cerr2,cerr2_alt(0:r),crmethod,crmethod_alt(0:r),rmax,r)
2157 ctyp = max(abs(c(0,0,0)),abs(c(0,1,0)),abs(c(0,0,1)))
2159 ctyp = abs(c(0,0,0))
2161 err_req_cr = acc_req_cr * ctyp
2164 write(*,*)
'CalcCred: after exps 2nd try Cmethod',crmethod
2165 write(*,*)
'CalcCred: after exps 2nd try Cerr(r)',cerr
2166 write(*,*)
'CalcCred: after exps 2nd try Cacc(r)',cerr/ctyp
2169 if(checkest.gt.impest_c.and.mode_coli.lt.1)
goto 100
2171 elseif (use_gp.and.err_gp(r).le.min(err_g(r),err_gy(r),err_gr(r),err_gpf(r)) &
2172 .and.mod(crcalc(r),32)-mod(crcalc(r),16).ne.16)
then
2186 call calccgp(c_alt,cuv,p10,p21,p20,m02,m12,m22,r,gp,rmax_c,id,cerr1_alt,acc_req_cr,cerr2_alt)
2188 call calccgp(c_alt(0:r,0:r,0:r),cuv_alt(0:r,0:r,0:r), &
2189 p10,p21,p20,m02,m12,m22,r,gp,rmax_c,id,cerr1_alt(0:r),acc_req_cr(0:r),cerr2_alt(0:r))
2192 cerr_alt = cerr2_alt
2194 cerr_alt = cerr1_alt
2197 crcalc(0:r)=crcalc(0:r)+16
2199 crmethod_alt(0:r)=16
2200 checkest=cerr_alt(r)/err_gp(r)
2203 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
2204 write(*,*)
'CalcCred: estimate err_gp imprecise',err_gp(r),cerr_alt(r)
2208 err_gp(0:r)=cerr_alt(0:r)
2210 call copycimp3(c,c_alt(0:r,0:r,0:r),cerr,cerr_alt(0:r),cerr1,cerr1_alt(0:r), &
2211 cerr2,cerr2_alt(0:r),crmethod,crmethod_alt(0:r),rmax,r)
2214 ctyp = max(abs(c(0,0,0)),abs(c(0,1,0)),abs(c(0,0,1)))
2216 ctyp = abs(c(0,0,0))
2218 err_req_cr = acc_req_cr * ctyp
2221 write(*,*)
'CalcCred: after exps 2nd try Cmethod',crmethod
2222 write(*,*)
'CalcCred: after exps 2nd try Cerr(r)',cerr
2223 write(*,*)
'CalcCred: after exps 2nd try Cacc(r)',cerr/ctyp
2226 if(checkest.gt.impest_c.and.mode_coli.lt.1)
goto 100
2228 elseif (use_gr.and.err_gr(r).le.min(err_g(r),err_gy(r),err_gp(r),err_gpf(r)) &
2229 .and.mod(crcalc(r),64)-mod(crcalc(r),32).ne.32)
then
2242 call calccgr(c_alt,cuv,p10,p21,p20,m02,m12,m22,r,gr,rmax_c,id,cerr1_alt,acc_req_cr,cerr2_alt)
2244 call calccgr(c_alt(0:r,0:r,0:r),cuv_alt(0:r,0:r,0:r), &
2245 p10,p21,p20,m02,m12,m22,r,gr,rmax_c,id,cerr1_alt(0:r),acc_req_cr(0:r),cerr2_alt(0:r))
2248 cerr_alt = cerr2_alt
2250 cerr_alt = cerr1_alt
2253 crcalc(0:r)=crcalc(0:r)+32
2255 crmethod_alt(0:r)=32
2256 checkest=cerr_alt(r)/err_gr(r)
2259 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
2260 write(*,*)
'CalcCred: estimate err_gr imprecise',err_gr(r),cerr_alt(r)
2264 err_gr(0:r)=cerr_alt(0:r)
2266 call copycimp3(c,c_alt(0:r,0:r,0:r),cerr,cerr_alt(0:r),cerr1,cerr1_alt(0:r), &
2267 cerr2,cerr2_alt(0:r),crmethod,crmethod_alt(0:r),rmax,r)
2270 ctyp = max(abs(c(0,0,0)),abs(c(0,1,0)),abs(c(0,0,1)))
2272 ctyp = abs(c(0,0,0))
2274 err_req_cr = acc_req_cr * ctyp
2277 write(*,*)
'CalcCred: after exps 2nd try Cmethod',crmethod
2278 write(*,*)
'CalcCred: after exps 2nd try Cerr(r)',cerr
2279 write(*,*)
'CalcCred: after exps 2nd try Cacc(r)',cerr/ctyp
2282 if(checkest.gt.impest_c.and.mode_coli.lt.1)
goto 100
2284 else if (use_gpf.and.err_gpf(r).le.min(err_g(r),err_gy(r),err_gp(r),err_gr(r)) &
2285 .and.mod(crcalc(r),128)-mod(crcalc(r),64).ne.64)
then
2298 call calccgpf(c_alt,cuv,p10,p21,p20,m02,m12,m22,r,gpf,(rmax_c)/2,id,cerr1_alt,acc_req_cr,cerr2_alt)
2300 call calccgpf(c_alt(0:r,0:r,0:r),cuv_alt(0:r,0:r,0:r), &
2301 p10,p21,p20,m02,m12,m22,r,gpf,(rmax_c)/2,id,cerr1_alt(0:r),acc_req_cr(0:r),cerr2_alt(0:r))
2304 cerr_alt = cerr2_alt
2306 cerr_alt = cerr1_alt
2309 crcalc(0:r)=crcalc(0:r)+64
2311 crmethod_alt(0:r)=64
2312 checkest=cerr_alt(r)/err_gpf(r)
2315 if(checkest.gt.1d2*impest_c.or.checkest.lt.1d-2*impest_c)
then
2316 write(*,*)
'CalcCred: estimate err_gpf imprecise',err_gpf(r),cerr_alt(r)
2319 err_gpf(0:r)=cerr_alt(0:r)
2321 call copycimp3(c,c_alt(0:r,0:r,0:r),cerr,cerr_alt(0:r),cerr1,cerr1_alt(0:r), &
2322 cerr2,cerr2_alt(0:r),crmethod,crmethod_alt(0:r),rmax,r)
2325 ctyp = max(abs(c(0,0,0)),abs(c(0,1,0)),abs(c(0,0,1)))
2327 ctyp = abs(c(0,0,0))
2329 err_req_cr = acc_req_cr * ctyp
2332 write(*,*)
'CalcCred: after exps 2nd try Cmethod',crmethod
2333 write(*,*)
'CalcCred: after exps 2nd try Cerr(r)',cerr
2334 write(*,*)
'CalcCred: after exps 2nd try Cacc(r)',cerr/ctyp
2337 if(checkest.gt.impest_c.and.mode_coli.lt.1)
goto 100
2346 if(.not.lerr_c0)
then
2348 err_c0 = acc_def_c0*max( abs(c(0,0,0)), 1d0/sqrt(
adetz) )
2356 write(*,*)
'CalcCred: final loop err methods',r,err_pv(r),err_pv2(r), &
2357 err_g(r),err_gy(r),err_gp(r),err_gr(r),err_gpf(r)
2358 write(*,*)
'CalcCred: final loop acc methods',r,err_pv(r)/ctyp,err_pv2(r)/ctyp, &
2359 err_g(r)/ctyp,err_gy(r)/ctyp,err_gp(r)/ctyp,err_gr(r)/ctyp,err_gpf(r)/ctyp
2360 write(*,*)
'CalcCred: final loop',r,crcalc(r),crmethod(r)
2364 norm = abs(c(0,0,0))
2368 norm = max(norm,abs(c(0,n1,n2)))
2371 acc_c(0:rdef) = cerr(0:rdef)/norm
2375 if (maxval(acc_c(0:rdef)-sqrt(reqacc_coli)).gt.0)
then
2379 if (maxval(acc_c(0:rdef)-reqacc_coli).gt.0)
then
2384 write(*,*)
'CalcCred final acc_C=',cerr/norm,critacc_coli
2385 write(*,*)
'CalcCred final method C=',crmethod
2388 if (maxval(acc_c(0:rdef)-critacc_coli).gt.0)
then
2393 write(*,*)
'CritPoint C',critacc_coli,acc_c
2394 write(*,*)
'CritPoint C',maxval(acc_c(0:rdef)-critacc_coli),maxval(acc_c(0:rdef)),rdef,acc_c(rdef)
2401 #ifdef CritPointsCOLI
2402 critpointcntc = critpointcntc + 1
2404 if (critpointcntc.le.maxcritpointc.and.monitoring)
then
2405 call critpointsout_coli(
'C_coli',acc_c(rdef))
2406 write(ncpout_coli,*)
'arguments of CalcCred_coli:'
2407 write(ncpout_coli,*)
'rank = ', rmax
2408 if(
present(rbasic))
write(ncpout_coli,*)
'rbas = ', rbasic
2409 write(ncpout_coli,fmt1)
'p10 = ', p10
2410 write(ncpout_coli,fmt1)
'p21 = ', p21
2411 write(ncpout_coli,fmt1)
'p20 = ', p20
2412 write(ncpout_coli,fmt1)
'm02 = ', m02
2413 write(ncpout_coli,fmt1)
'm12 = ', m12
2414 write(ncpout_coli,fmt1)
'm22 = ', m22
2415 if (critpointcntc.eq.maxcritpointc)
then
2416 write(ncpout_coli,*)
2417 write(ncpout_coli,*)
2418 write(ncpout_coli,*)
2419 write(ncpout_coli,*)
'***********************************************************'
2420 write(ncpout_coli,*)
2421 write(ncpout_coli,*)
' Further output of bad C functions will be suppressed '
2428 write(*,*)
'CalcCred out'
2441 subroutine calccuv(Cuv,Buv_0,m02,f,rmax,id)
2443 integer,
intent(in) :: rmax,id
2444 double complex,
intent(in) :: m02,f(2)
2446 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
2447 double complex,
intent(in) :: Buv_0(0:rmax-1,0:rmax-1,0:rmax-1)
2448 integer :: r,n0,n1,n2,r0
2458 do n0=max(1,r-rmax),r/2
2462 cuv(n0,n1,n2) = (buv_0(n0-1,n1,n2) + 2*m02*cuv(n0-1,n1,n2) &
2463 + f(1)*cuv(n0-1,n1+1,n2) &
2464 + f(2)*cuv(n0-1,n1,n2+1)) / (2*r)
2486 subroutine calccpv1(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,id,Cerr,Cerr2)
2490 integer,
intent(in) :: rmax,id
2491 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
2492 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
2493 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
2494 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
2498 double complex :: B_0(0:rmax-1,0:rmax-1,0:rmax-1), Buv_0(0:rmax-1,0:rmax-1,0:rmax-1)
2499 double complex :: B_i(0:rmax-1,0:rmax-1,2), Buv_i(0:rmax-1,0:rmax-1,2)
2500 double complex :: C_alt(0:rmax,0:rmax,0:rmax)
2501 double complex :: Smod(2)
2502 double complex :: C0_coli, elimminf2_coli
2505 double precision :: C00_err(0:rmax),Cij_err(0:rmax)
2506 double precision :: C00_err2(0:rmax),Cij_err2(0:rmax)
2507 double precision :: B_err,B_max
2508 integer :: rmaxB,r,n0,n1,n2,nn0,nn1,nn2,i,j
2509 integer :: bin,k,nid(0:2)
2514 write(*,*)
'CalcCpv1 in ',rmax,id
2517 write(*,*)
'CalcCpv1 in ',rmax,id
2521 c(0,0,0) = c0_coli(p10,p21,p20,m02,m12,m22)
2525 cerr(0) = acc_def_c0*max(1d0/sqrt(
adetz),abs(c(0,0,0)))
2526 cerr2(0) = acc_def_c0*max(1d0/sqrt(
adetz),abs(c(0,0,0)))
2528 if (rmax.eq.0)
return
2548 if (mod(id/bin,2).eq.0)
then
2555 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
2556 call calcb(b_i(:,:,1),buv_i(:,:,1),p20,m02,m22,rmaxb,nid(1))
2557 call calcb(b_i(:,:,2),buv_i(:,:,2),p10,m02,m12,rmaxb,nid(2))
2564 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
2565 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
2566 b_max = max(b_max,abs(b_0(0,n1,n2)))
2572 b_max=max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
2584 cij_err(0) = cerr(0)
2585 b_err = acc_def_b*b_max
2589 cij_err2(0) = cerr2(0)
2594 write(*,*)
'CalcDpv1 Cij_err(0)=',cij_err(0)
2609 c(n0,n1,n2) = + 4*cuv(n0,n1,n2) +
detx/
detz*c(n0-1,n1,n2)
2610 c(n0,n1,n2) = c(n0,n1,n2) &
2612 ) /
detz * b_0(n0-1,n1,n2)
2616 c(n0,n1,n2) = c(n0,n1,n2) &
2619 c(n0,n1,n2) = c(n0,n1,n2) &
2623 c(n0,n1,n2) = c(n0,n1,n2) &
2626 c(n0,n1,n2) = c(n0,n1,n2) &
2630 c(n0,n1,n2) = c(n0,n1,n2) / (2*r)
2668 smod(1) = smod(1) - 2d0*nn1*c(n0+1,nn1-1,nn2)
2670 smod(1) = smod(1) + b_i(n0,nn2,1)
2674 smod(2) = smod(2) - 2d0*nn2*c(n0+1,nn1,nn2-1)
2676 smod(2) = smod(2) + b_i(n0,nn1,2)
2679 c(n0,n1,n2) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
2680 -
zadjs(j)*b_0(n0,nn1,nn2) &
2685 write(*,*)
'Ca(0,n1,n2)=',n1,n2,c(0,n1,n2),nn1,nn2,j
2686 write(*,*)
'Ca(0,n1,n2)=',
zadj(1,j),smod(1),
zadj(2,j),smod(2)
2687 write(*,*)
'Ca(0,n1,n2)=',
zadjs(j),b_0(n0,nn1,nn2),
zadjf(j),c(n0,nn1,nn2)
2688 write(*,*)
'Ca(0,n1,n2)=',
zadj(1,j)*smod(1),
zadj(2,j)*smod(2)
2689 write(*,*)
'Ca(0,n1,n2)=',-
zadjs(j)*b_0(n0,nn1,nn2),-
zadjf(j)*c(n0,nn1,nn2)
2713 smod(1) = smod(1) - 2d0*nn1*c(n0+1,nn1-1,nn2)
2715 smod(1) = smod(1) + b_i(n0,nn2,1)
2719 smod(2) = smod(2) - 2d0*nn2*c(n0+1,nn1,nn2-1)
2721 smod(2) = smod(2) + b_i(n0,nn1,2)
2724 c_alt(n0,n1,n2) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
2725 -
zadjs(j)*b_0(n0,nn1,nn2) &
2728 cerr(r)=max(cerr(r),abs(c(n0,n1,n2)-c_alt(n0,n1,n2)))
2729 cerr2(r)=max(cerr2(r),abs(c(n0,n1,n2)-c_alt(n0,n1,n2)))
2733 write(*,*)
'Cb(0,n1,n2)=',n1,n2,c_alt(0,n1,n2),nn1,nn2,j
2734 write(*,*)
'Cb(0,n1,n2)=',
zadj(1,j),smod(1),
zadj(2,j),smod(2)
2735 write(*,*)
'Cb(0,n1,n2)=',
zadjs(j),b_0(n0,nn1,nn2),
zadjf(j),c(n0,nn1,nn2)
2744 c00_err(r) = max(2*abs(m02)*cij_err(r-2), b_err, &
2757 cij_err(r) = max(
maxzadjf*cij_err(r-1), &
2761 c00_err2(r) = max(2*abs(m02)*cij_err2(r-2), b_err, &
2781 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
2782 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
2789 write(*,*)
'CalcCpv1 Cerrsym',cerr
2790 write(*,*)
'CalcCpv1 Caccsym',cerr/abs(c(0,0,0))
2792 write(*,*)
'CalcCpv1 Cijerr',cij_err(1:rmax)
2793 write(*,*)
'CalcCpv1 Cijacc',cij_err(1:rmax)/abs(c(0,0,0))
2796 cerr2 = max(cerr2,cij_err2(0:rmax))
2797 cerr = max(cerr,cij_err(0:rmax))
2800 write(*,*)
'CalcCpv1 Cerr',cerr
2801 write(*,*)
'CalcCpv1 Cacc',cerr/abs(c(0,0,0))
2825 subroutine calccpv1o(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,id,Cerr,Cerr2)
2829 integer,
intent(in) :: rmax,id
2830 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
2831 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
2832 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
2833 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
2834 double complex,
allocatable :: B_0(:,:,:), Buv_0(:,:,:)
2835 double complex,
allocatable :: B_i(:,:,:), Buv_i(:,:,:)
2836 double complex,
allocatable :: C_alt(:,:,:)
2837 double complex :: Smod(2)
2838 double complex :: C0_coli, elimminf2_coli
2839 double precision,
allocatable :: C00_err(:),Cij_err(:)
2840 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
2841 double precision :: B_err,B_max
2842 integer :: rmaxB,r,n0,n1,n2,nn0,nn1,nn2,i,j
2843 integer :: bin,k,nid(0:2)
2846 write(*,*)
'CalcCpv1o in ',rmax,id
2849 write(*,*)
'CalcCpv1o in ',rmax,id
2853 c(0,0,0) = c0_coli(p10,p21,p20,m02,m12,m22)
2857 cerr(0) = acc_def_c0*max(1d0/sqrt(
adetz),abs(c(0,0,0)))
2858 cerr2(0) = acc_def_c0*max(1d0/sqrt(
adetz),abs(c(0,0,0)))
2860 if (rmax.eq.0)
return
2865 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
2866 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
2867 allocate(b_i(0:rmaxb,0:rmaxb,2))
2868 allocate(buv_i(0:rmaxb,0:rmaxb,2))
2871 allocate(c00_err(0:rmax))
2872 allocate(cij_err(0:rmax))
2873 allocate(c00_err2(0:rmax))
2874 allocate(cij_err2(0:rmax))
2880 if (mod(id/bin,2).eq.0)
then
2887 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
2888 call calcb(b_i(:,:,1),buv_i(:,:,1),p20,m02,m22,rmaxb,nid(1))
2889 call calcb(b_i(:,:,2),buv_i(:,:,2),p10,m02,m12,rmaxb,nid(2))
2896 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
2897 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
2898 b_max = max(b_max,abs(b_0(0,n1,n2)))
2904 b_max=max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
2948 cij_err(0) = cerr(0)
2949 b_err = acc_def_b*b_max
2953 cij_err2(0) = cerr2(0)
2957 allocate(c_alt(0:rmax,0:rmax,0:rmax))
2966 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
2967 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
2988 smod(i) = -b_0(n0,nn1,nn2)-
f(i)*c(n0,nn1,nn2)
2992 smod(1) = smod(1) - 2d0*nn1*c(n0+1,nn1-1,nn2)
2994 smod(1) = smod(1) + b_i(n0,nn2,1)
2998 smod(2) = smod(2) - 2d0*nn2*c(n0+1,nn1,nn2-1)
3000 smod(2) = smod(2) + b_i(n0,nn1,2)
3003 c(n0,n1,n2) =
zinv(1,j)*smod(1) +
zinv(2,j)*smod(2)
3026 smod(i) = -b_0(n0,nn1,nn2)-
f(i)*c(n0,nn1,nn2)
3030 smod(1) = smod(1) - 2d0*nn1*c(n0+1,nn1-1,nn2)
3032 smod(1) = smod(1) + b_i(n0,nn2,1)
3036 smod(2) = smod(2) - 2d0*nn2*c(n0+1,nn1,nn2-1)
3038 smod(2) = smod(2) + b_i(n0,nn1,2)
3041 c_alt(n0,n1,n2) =
zinv(1,j)*smod(1) +
zinv(2,j)*smod(2)
3043 cerr(r)=max(cerr(r),abs(c(n0,n1,n2)-c_alt(n0,n1,n2)))
3044 cerr2(r)=max(cerr2(r),abs(c(n0,n1,n2)-c_alt(n0,n1,n2)))
3055 c00_err(r) = max(2*abs(m02)*cij_err(r-2), b_err, &
3067 cij_err(r) = max(
maxzadjf*cij_err(r-1), &
3071 c00_err2(r) = max(2*abs(m02)*cij_err2(r-2), b_err, &
3091 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
3092 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
3099 write(*,*)
'CalcCpv1o Cerrsym',cerr
3100 write(*,*)
'CalcCpv1o Caccsym',cerr/abs(c(0,0,0))
3102 write(*,*)
'CalcCpv1o Cijerr',cij_err(1:rmax)
3103 write(*,*)
'CalcCpv1o Cijacc',cij_err(1:rmax)/abs(c(0,0,0))
3106 cerr2 = max(cerr2,cij_err2(0:rmax))
3107 cerr = max(cerr,cij_err(0:rmax))
3110 write(*,*)
'CalcCpv1o Cerr',cerr
3111 write(*,*)
'CalcCpv1o Cacc',cerr/abs(c(0,0,0))
3129 subroutine calccpv(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,id,Cerr,Cerr2)
3133 integer,
intent(in) :: rmax,id
3134 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
3135 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
3136 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
3137 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
3138 double complex,
allocatable :: B_0(:,:,:), Buv_0(:,:,:)
3139 double complex,
allocatable :: B_i(:,:,:), Buv_i(:,:,:)
3140 double complex,
allocatable :: C_alt(:,:,:)
3141 double complex :: Smod(2)
3142 double complex :: C0_coli, elimminf2_coli
3143 double precision,
allocatable :: C00_err(:),Cij_err(:)
3144 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
3145 double precision :: B_err,B_max
3146 integer :: rmaxB,r,n0,n1,n2,nn0,nn1,nn2,i,j
3147 integer :: bin,k,nid(0:2)
3150 write(*,*)
'CalcCpv in ', id
3153 write(*,*)
'CalcCpv !n ', id
3157 c(0,0,0) = c0_coli(p10,p21,p20,m02,m12,m22)
3161 cerr(0) = acc_def_c0*max(1d0/sqrt(
adetz),abs(c(0,0,0)))
3162 cerr2(0) = acc_def_c0*max(1d0/sqrt(
adetz),abs(c(0,0,0)))
3166 if (rmax.eq.0)
return
3171 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
3172 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
3173 allocate(b_i(0:rmaxb,0:rmaxb,2))
3174 allocate(buv_i(0:rmaxb,0:rmaxb,2))
3177 allocate(c00_err(0:rmax))
3178 allocate(cij_err(0:rmax))
3179 allocate(c00_err2(0:rmax))
3180 allocate(cij_err2(0:rmax))
3186 if (mod(id/bin,2).eq.0)
then
3193 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
3194 call calcb(b_i(:,:,1),buv_i(:,:,1),p20,m02,m22,rmaxb,nid(1))
3195 call calcb(b_i(:,:,2),buv_i(:,:,2),p10,m02,m12,rmaxb,nid(2))
3202 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
3203 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
3204 b_max = max(b_max,abs(b_0(0,n1,n2)))
3210 b_max=max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
3254 cij_err(0) = cerr(0)
3255 b_err = acc_def_b*b_max
3259 cij_err2(0) = cerr2(0)
3263 allocate(c_alt(0:rmax,0:rmax,0:rmax))
3268 if (mod(r,2).eq.0)
then
3271 c(n0,0,0) = (b_0(n0-1,0,0) + 2*
mm02*c(n0-1,0,0) + 4*cuv(n0,0,0) &
3272 +
f(1)*c(n0-1,1,0) +
f(2)*c(n0-1,0,1)) / (2*r)
3291 smod(i) = -b_0(n0,nn1,nn2)-
f(i)*c(n0,nn1,nn2)
3295 smod(1) = smod(1) - 2d0*nn1*c(n0+1,nn1-1,nn2)
3297 smod(1) = smod(1) + b_i(n0,nn2,1)
3301 smod(2) = smod(2) - 2d0*nn2*c(n0+1,nn1,nn2-1)
3303 smod(2) = smod(2) + b_i(n0,nn1,2)
3306 c(n0,n1,n2) =
zinv(1,j)*smod(1) +
zinv(2,j)*smod(2)
3323 if (n1.ge.1.and.n2.ge.1)
then
3330 smod(i) = -b_0(n0,nn1,nn2)-
f(i)*c(n0,nn1,nn2)
3334 smod(1) = smod(1) - 2d0*nn1*c(n0+1,nn1-1,nn2)
3336 smod(1) = smod(1) + b_i(n0,nn2,1)
3340 smod(2) = smod(2) - 2d0*nn2*c(n0+1,nn1,nn2-1)
3342 smod(2) = smod(2) + b_i(n0,nn1,2)
3345 c_alt(n0,n1,n2) =
zinv(1,j)*smod(1) +
zinv(2,j)*smod(2)
3347 cerr(r)=max(cerr(r),abs(c(n0,n1,n2)-c_alt(n0,n1,n2)))
3348 cerr2(r)=max(cerr2(r),abs(c(n0,n1,n2)-c_alt(n0,n1,n2)))
3354 c00_err(r) = max(abs(m02)*cij_err(r-2), b_err, &
3366 cij_err(r) = max(
maxzadjf*cij_err(r-1), &
3370 c00_err2(r) = max(abs(m02)*cij_err(r-2), b_err, &
3383 max(c00_err2(r),b_err))/sqrt(
adetz)
3398 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
3399 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
3405 write(*,*)
'CalcCpv Cerrsym',cerr
3406 write(*,*)
'CalcCpv Caccsym',cerr/abs(c(0,0,0))
3408 write(*,*)
'CalcCpv Cijerr',cij_err(1:rmax)
3409 write(*,*)
'CalcCpv Cijacc',cij_err(1:rmax)/abs(c(0,0,0))
3412 cerr2 = max(cerr2,cij_err2(0:rmax))
3413 cerr = max(cerr,cij_err(0:rmax))
3416 write(*,*)
'CalcCpv Cerr',cerr
3417 write(*,*)
'CalcCpv Cacc',cerr/abs(c(0,0,0))
3433 subroutine calccpv2(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,id,Cerr,Cerr2)
3437 integer,
intent(in) :: rmax,id
3438 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
3439 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
3440 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
3441 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
3442 double complex,
allocatable :: B_0(:,:,:), B_i(:,:,:)
3443 double complex,
allocatable :: Buv_0(:,:,:), Buv_i(:,:,:)
3444 double complex,
allocatable :: C_alt(:,:,:)
3445 double complex :: C0_coli, elimminf2_coli
3447 double complex :: Caux(1:rmax,0:rmax-1,0:rmax-1), Smod(2)
3448 double complex :: chdet
3449 double precision,
allocatable :: C00_err(:),Cij_err(:)
3450 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
3451 double precision :: B_err,B_max
3452 integer :: rmaxB,r,n0,n1,n2,k
3453 integer :: bin,nid(0:3)
3456 write(*,*)
'CalcCpv2 in '
3459 write(*,*)
'CalcCpv2 in '
3463 c(0,0,0) = c0_coli(p10,p21,p20,m02,m12,m22)
3467 cerr(0) = acc_def_c0*max( abs(c(0,0,0)), 1d0/sqrt(
adetz) )
3468 cerr2(0) = acc_def_c0*max( abs(c(0,0,0)), 1d0/sqrt(
adetz) )
3472 if (rmax.eq.0)
return
3478 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
3479 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
3480 allocate(b_i(0:rmaxb,0:rmaxb,2))
3481 allocate(buv_i(0:rmaxb,0:rmaxb,2))
3484 allocate(c00_err(0:rmax+1))
3485 allocate(cij_err(0:rmax))
3486 allocate(c00_err2(0:rmax+1))
3487 allocate(cij_err2(0:rmax))
3493 if (mod(id/bin,2).eq.0)
then
3500 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
3501 call calcb(b_i(:,:,1),buv_i(:,:,1),p20,m02,m22,rmaxb,nid(1))
3502 call calcb(b_i(:,:,2),buv_i(:,:,2),p10,m02,m12,rmaxb,nid(2))
3509 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
3510 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
3511 b_max = max(b_max,abs(b_0(0,n1,n2)))
3514 b_max=max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
3539 cij_err(0) = cerr(0)
3540 b_err = acc_def_b*b_max
3544 cij_err2(0) = cerr2(0)
3548 allocate(c_alt(0:rmax,0:rmax,0:rmax))
3558 smod(k) = -b_0(n0-1,n1,n2)
3562 smod(1) = smod(1) - 2*n1*c(n0,n1-1,n2)
3564 smod(1) = smod(1) + b_i(n0-1,n2,1)
3568 smod(2) = smod(2) - 2*n2*c(n0,n1,n2-1)
3570 smod(2) = smod(2) + b_i(n0-1,n1,2)
3573 caux(n0,n1,n2) = (c(n0-1,n1,n2) -
mxinv(1,0)*smod(1) &
3584 c(n0,n1,n2) = (caux(n0,n1,n2) + 4d0*cuv(n0,n1,n2) &
3585 + b_0(n0-1,n1,n2))/r/2d0
3600 smod(k) = -b_0(0,n1,n2)
3604 smod(1) = smod(1) - 2*n1*c(1,n1-1,n2)
3606 smod(1) = smod(1) + b_i(0,n2,1)
3610 smod(2) = smod(2) - 2*n2*c(1,n1,n2-1)
3612 smod(2) = smod(2) + b_i(0,n1,2)
3615 caux(1,n1,n2) = (c(0,n1,n2) -
mxinv(1,0)*smod(1) &
3618 c(0,n1+1,n2) =
mxinv(0,1)*caux(1,n1,n2) &
3620 c_alt(0,n1,n2+1) =
mxinv(0,2)*caux(1,n1,n2) &
3624 c(0,0,r) = c_alt(0,0,r)
3626 cerr(r)=max(cerr(r),abs(c(0,n1,n2+1)-c_alt(0,n1,n2+1)))
3627 cerr2(r)=max(cerr2(r),abs(c(0,n1,n2+1)-c_alt(0,n1,n2+1)))
3632 c00_err(r+1) = max(b_err,
adetx/
adetz*cij_err(r-1), &
3638 cij_err(r) = max(
maxzadjf*max(2*(r+1)*c00_err(r+1),b_err), &
3641 c00_err2(r+1) = max(b_err,
adetx/
adetz*cij_err2(r-1), &
3647 cij_err2(r) = max(
maxzadjf*max(2*(r+1)*c00_err2(r+1),b_err), &
3664 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
3665 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
3667 write(*,*)
'C1(n0+1)',n0,n1,n2
3668 write(*,*)
'C1(n0+1)',(b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
3669 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
3676 do n0=max(2,r-rmax),r/2
3681 smod(k) = -b_0(n0-1,n1,n2)
3685 smod(1) = smod(1) - 2*n1*c(n0,n1-1,n2)
3687 smod(1) = smod(1) + b_i(n0-1,n2,1)
3691 smod(2) = smod(2) - 2*n2*c(n0,n1,n2-1)
3693 smod(2) = smod(2) + b_i(n0-1,n1,2)
3696 caux(n0,n1,n2) = (c(n0-1,n1,n2) -
mxinv(1,0)*smod(1) &
3707 c(n0,n1,n2) = (caux(n0,n1,n2) + 4d0*cuv(n0,n1,n2) &
3708 + b_0(n0-1,n1,n2))/r/2d0
3711 write(*,*)
'C2(n0+1)',n0,n1,n2
3712 write(*,*)
'C2(n0+1)',(caux(n0,n1,n2) + 4d0*cuv(n0,n1,n2) &
3713 + b_0(n0-1,n1,n2))/r/2d0
3722 write(*,*)
'CalcCpv2 Cerrsym',cerr
3723 write(*,*)
'CalcCpv2 Caccsym',cerr/abs(c(0,0,0))
3725 write(*,*)
'CalcCpv2 Cijerr',cij_err
3726 write(*,*)
'CalcCpv2 Cijacc',cij_err/abs(c(0,0,0))
3729 cerr2 = max(cerr2,cij_err2(0:rmax))
3730 cerr = max(cerr,cij_err(0:rmax))
3733 write(*,*)
'CalcCpv2 Cerr',cerr
3734 write(*,*)
'CalcCpv2 Cacc',cerr/abs(c(0,0,0))
3754 subroutine calccpvshift(Cshift,Cuvshift,p10shift,p21shift,p20shift,m02shift,m12shift,m22shift,rmax,id,Cerr,Cerr2)
3758 integer,
intent(in) :: rmax,id
3759 double complex,
intent(in) :: p10shift,p21shift,p20shift,m02shift,m12shift,m22shift
3760 double complex,
intent(out) :: Cuvshift(0:rmax,0:rmax,0:rmax)
3761 double complex,
intent(out) :: Cshift(0:rmax,0:rmax,0:rmax)
3762 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
3763 double complex,
allocatable :: B_0(:,:,:), Buv_0(:,:,:)
3764 double complex,
allocatable :: B_i(:,:,:), Buv_i(:,:,:)
3765 double complex,
allocatable :: Cshift_alt(:,:,:)
3766 double complex :: Smod(2)
3767 double complex :: C0_coli, elimminf2_coli
3768 double precision,
allocatable :: C00_err(:),Cij_err(:)
3769 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
3770 double precision :: B_err,B_max
3771 integer :: rmaxB,r,n0,n1,n2,nn0,nn1,nn2,i,j
3772 integer :: bin,k,nid(0:2)
3773 logical :: use_cache_system_save
3776 write(*,*)
'CalcCpvshift p in ',p10shift,p21shift,p20shift
3777 write(*,*)
'CalcCpvshift m in ',m02shift,m12shift,m22shift
3778 write(*,*)
'CalcCpvshift in ',rmax,id
3781 write(*,*)
'CalcCpvshift in ',rmax,id
3785 cshift(0,0,0) = c0_coli(p10shift,p21shift,p20shift,m02shift,m12shift,m22shift)
3786 cuvshift(0,0,0) = 0d0
3789 cerr(0) = acc_def_c0*max(1d0/sqrt(
adetzshift),abs(cshift(0,0,0)))
3790 cerr2(0) = acc_def_c0*max(1d0/sqrt(
adetzshift),abs(cshift(0,0,0)))
3792 if (rmax.eq.0)
return
3797 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
3798 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
3799 allocate(b_i(0:rmaxb,0:rmaxb,2))
3800 allocate(buv_i(0:rmaxb,0:rmaxb,2))
3803 allocate(c00_err(0:rmax))
3804 allocate(cij_err(0:rmax))
3805 allocate(c00_err2(0:rmax))
3806 allocate(cij_err2(0:rmax))
3812 if (mod(id/bin,2).eq.0)
then
3822 use_cache_system_save = use_cache_system
3823 use_cache_system = .false.
3824 call calcb(b_0(:,0,:),buv_0(:,0,:),p21shift,m12shift,m22shift,rmaxb,0)
3825 call calcb(b_i(:,:,1),buv_i(:,:,1),p20shift,m02shift,m22shift,rmaxb,0)
3826 call calcb(b_i(:,:,2),buv_i(:,:,2),p10shift,m02shift,m12shift,rmaxb,0)
3827 use_cache_system = use_cache_system_save
3836 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
3837 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
3838 b_max = max(b_max,abs(b_0(0,n1,n2)))
3844 b_max=max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
3855 cij_err(0) = cerr(0)
3856 b_err = acc_def_b*b_max
3860 cij_err2(0) = cerr2(0)
3865 write(*,*)
'CalcDpvshift Cij_err(0)=',cij_err(0)
3871 allocate(cshift_alt(0:rmax,0:rmax,0:rmax))
3881 cshift(n0,n1,n2) = cshift(n0,n1,n2) &
3887 cshift(n0,n1,n2) = cshift(n0,n1,n2) &
3890 cshift(n0,n1,n2) = cshift(n0,n1,n2) &
3894 cshift(n0,n1,n2) = cshift(n0,n1,n2) &
3897 cshift(n0,n1,n2) = cshift(n0,n1,n2) &
3901 cshift(n0,n1,n2) = cshift(n0,n1,n2) / (2*r)
3939 smod(1) = smod(1) - 2d0*nn1*cshift(n0+1,nn1-1,nn2)
3941 smod(1) = smod(1) + b_i(n0,nn2,1)
3945 smod(2) = smod(2) - 2d0*nn2*cshift(n0+1,nn1,nn2-1)
3947 smod(2) = smod(2) + b_i(n0,nn1,2)
3982 smod(1) = smod(1) - 2d0*nn1*cshift(n0+1,nn1-1,nn2)
3984 smod(1) = smod(1) + b_i(n0,nn2,1)
3988 smod(2) = smod(2) - 2d0*nn2*cshift(n0+1,nn1,nn2-1)
3990 smod(2) = smod(2) + b_i(n0,nn1,2)
3997 cerr(r)=max(cerr(r),abs(cshift(n0,n1,n2)-cshift_alt(n0,n1,n2)))
3998 cerr2(r)=max(cerr2(r),abs(cshift(n0,n1,n2)-cshift_alt(n0,n1,n2)))
4011 c00_err(r) = max(2*abs(m02shift)*cij_err(r-2), b_err, &
4023 c00_err2(r) = max(2*abs(m02shift)*cij_err2(r-2), b_err, &
4038 cshift(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02shift*cshift(n0-1,n1,n2) + 4*cuvshift(n0,n1,n2) &
4039 +
fshift(1)*cshift(n0-1,n1+1,n2) +
fshift(2)*cshift(n0-1,n1,n2+1)) / (2*r)
4046 write(*,*)
'CalcCpvshift Cerrsym',cerr
4047 write(*,*)
'CalcCpvshift Caccsym',cerr/abs(cshift(0,0,0))
4049 write(*,*)
'CalcCpvshift Cijerr',cij_err(1:rmax)
4050 write(*,*)
'CalcCpvshift Cijacc',cij_err(1:rmax)/abs(cshift(0,0,0))
4053 cerr2 = max(cerr2,cij_err2(0:rmax))
4054 cerr = max(cerr,cij_err(0:rmax))
4057 write(*,*)
'CalcCpvshift Cerr',cerr
4058 write(*,*)
'CalcCpvshift Cacc',cerr/abs(cshift(0,0,0))
4065 write(*,*)
'CalcCpvshift C',cshift(0,0,0)
4066 write(*,*)
'CalcCpvshift C1',cshift(0,1,0)
4067 write(*,*)
'CalcCpvshift C2',cshift(0,0,1)
4068 write(*,*)
'CalcCpvshift C11',cshift(0,2,0)
4069 write(*,*)
'CalcCpvshift C12',cshift(0,1,1)
4070 write(*,*)
'CalcCpvshift C22',cshift(0,0,2)
4088 subroutine calccgn(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,ordg_min,ordg_max,id,Cerr,acc_req_Cr,Cerr2)
4092 integer,
intent(in) :: rmax,ordg_min,ordg_max,id
4093 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
4094 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
4095 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
4096 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
4097 double precision,
intent(in) :: acc_req_Cr(0:rmax)
4098 double complex :: Xtilde,Zkl,Zadjfj,Zadj2,Zadjkl
4099 double complex,
allocatable :: Cexpg(:,:,:,:), CuvExpg(:,:,:)
4100 double complex,
allocatable :: B_0(:,:,:), B_i(:,:,:), Shat(:,:,:,:)
4101 double complex,
allocatable :: Buv_0(:,:,:), Buv_i(:,:,:)
4102 double complex :: Smod(2), Skl, CexpgAux
4103 double complex :: C0_coli, elimminf2_coli
4104 double precision,
allocatable :: C00_err(:),Cij_err(:)
4105 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
4106 double precision :: B_err,B_max
4107 double precision :: maxCexpg(0:1,0:rmax+ordg_min+1,0:ordg_max),truncfacexp
4108 integer :: rmaxB,rmaxExp,gtrunc,r,n0,n1,n2,k,l,i,j,m,n,sgn,g,rg
4109 integer :: inds0(2), inds(2), inds2(2), ktlt(2)
4110 integer :: bin,nid(0:2)
4112 double complex,
allocatable :: D_alt(:,:,:,:)
4115 write(*,*)
'CalcCgn in ',rmax,ordg_min,ordg_max
4118 write(*,*)
'CalcCgn in ',rmax,ordg_min,ordg_max
4122 rmaxb = rmax + ordg_min
4123 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
4124 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
4125 allocate(b_i(0:rmaxb,0:rmaxb,2))
4126 allocate(buv_i(0:rmaxb,0:rmaxb,2))
4133 if (mod(id/bin,2).eq.0)
then
4140 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
4141 call calcb(b_i(:,:,1),buv_i(:,:,1),p20,m02,m22,rmaxb,nid(1))
4142 call calcb(b_i(:,:,2),buv_i(:,:,2),p10,m02,m12,rmaxb,nid(2))
4149 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
4150 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
4151 b_max = max(b_max,abs(b_0(0,n1,n2)))
4155 b_max = max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
4156 b_err = acc_def_b*b_max
4193 allocate(shat(0:rmaxb,0:rmaxb,0:rmaxb,2))
4200 shat(n0,n1,n2,:) = -b_0(n0,n1,n2)
4204 shat(n0,0,k,1) = shat(n0,0,k,1) + b_i(n0,k,1)
4205 shat(n0,k,0,2) = shat(n0,k,0,2) + b_i(n0,k,2)
4243 allocate(cexpg(0:rmaxexp/2,0:rmaxexp-1,0:rmaxexp-1,0:ordg_max))
4246 allocate(cuvexpg(0:rmaxexp,0:rmaxexp,0:rmaxexp))
4248 cuv(0:rmax,0:rmax,0:rmax) = cuvexpg(0:rmax,0:rmax,0:rmax)
4251 allocate(c00_err(0:rmaxexp))
4252 allocate(cij_err(0:rmaxexp))
4253 allocate(c00_err2(0:rmaxexp))
4254 allocate(cij_err2(0:rmaxexp))
4266 write(*,*)
'CalcCgn rmax = ',rmax,rmaxexp
4267 write(*,*)
'CalcCgn Cij_err = ',cij_err
4268 write(*,*)
'CalcCgn B0 = ', b_0(0,0,0),b_i(0,0,1),b_i(0,0,2)
4280 rloop:
do r=1,rmaxexp
4286 if (r.gt.rmax+gtrunc+1)
exit rloop
4289 write(*,*)
'CalcCgn rloop',r
4307 cexpgaux = 2d0*
zadj(k,l)*b_0(n0-1,n1,n2) &
4308 + xtilde*cexpg(n0-1,n1,n2,0) &
4309 + 4d0*
zadj(k,l)*cuvexpg(n0,n1,n2)
4314 cexpgaux = cexpgaux +
zadj(i,l)*shat(n0-1,inds(1),inds(2),i)
4320 cexpgaux = cexpgaux -
zadj(k,l)*shat(n0-1,inds(1),inds(2),i)
4326 skl =
f(n)*shat(n0-1,inds0(1),inds0(2),m)
4329 if (inds(m).ge.1)
then
4331 skl = skl - 2d0*
f(n)*inds0(m)*cexpg(n0,inds(1),inds(2),0)
4332 if (inds(n).ge.1)
then
4334 skl = skl - 4d0*inds0(m)*(inds(n)+1)*cexpg(n0+1,inds(1),inds(2),0)
4338 if (inds(n).ge.1)
then
4340 skl = skl + 2d0*inds0(n)*shat(n0,inds(1),inds(2),m) &
4341 - 2d0*
f(m)*inds0(n)*cexpg(n0,inds(1),inds(2),0)
4344 cexpgaux = cexpgaux - zadj2*skl
4346 cexpg(n0,n1,n2,0) = cexpgaux/(2d0*zadjkl)/(2d0*(r-n0)+1)
4349 maxcexpg(1,r,0) = maxcexpg(1,r,0) + abs(cexpg(n0,n1,n2,0) )
4352 if (r-n0.le.rmax)
then
4353 c(n0,n1,n2) = cexpg(n0,n1,n2,0)
4362 maxcexpg(0,r-1,0)=0d0
4366 smod = shat(0,n1,n2,:)
4368 smod(1) = smod(1) - 2d0*n1*cexpg(1,n1-1,n2,0)
4371 smod(2) = smod(2) - 2d0*n2*cexpg(1,n1,n2-1,0)
4374 cexpg(0,n1,n2,0) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2))/zadjfj
4375 maxcexpg(0,r-1,0) = maxcexpg(0,r-1,0) + abs(cexpg(0,n1,n2,0))
4376 if (r-n0.le.rmax+1)
then
4377 c(0,n1,n2) = cexpg(0,n1,n2,0)
4392 if(r.le.rmax+1)
then
4394 cerr(r-1) =
fac_g*maxcexpg(0,r-1,0)
4398 c00_err(r) = max(max(
maxzadj*b_err,
fmax*b_err)/abs(zadjkl),b_err) &
4400 cij_err(r-1)=
maxzadj*max(b_err,2*c00_err(r))/abs(zadjfj)
4402 c00_err2(r) = max(max(
maxzadj*b_err,
fmax*b_err)/abs(zadjkl),b_err) &
4404 cij_err2(r-1)=
maxzadj*max(b_err,2*c00_err2(r))/abs(zadjfj)
4407 write(*,*)
'CalcCgn C00_err',r,
maxzadj,
fmax,abs(zadjkl),b_err,abs(zadjfj)
4408 write(*,*)
'CalcCgn C00_err',r, c00_err(r), cij_err(r-1)
4416 gloop:
do g=1,min(gtrunc,r-1)
4424 maxcexpg(1,rg,g) = 0d0
4435 cexpgaux = xtilde*cexpg(n0-1,n1,n2,g) &
4436 -
detz*cexpg(n0-1,inds(1),inds(2),g-1)
4445 if (inds(m).ge.1)
then
4447 skl = skl - 2d0*
f(n)*inds0(m)*cexpg(n0,inds(1),inds(2),g)
4448 if (inds(n).ge.1)
then
4450 skl = skl - 4d0*inds0(m)*(inds(n)+1)*cexpg(n0+1,inds(1),inds(2),g)
4454 if (inds(n).ge.1)
then
4456 skl = skl - 2d0*
f(m)*inds0(n)*cexpg(n0,inds(1),inds(2),g)
4459 cexpgaux = cexpgaux - zadj2*skl
4461 cexpg(n0,n1,n2,g) = cexpgaux/(2d0*zadjkl)/(2d0*(rg-n0)+1)
4465 maxcexpg(1,rg,g) = maxcexpg(1,rg,g) + abs(cexpg(n0,n1,n2,g))
4467 if (g.eq.1.and.abs(cexpg(n0,n1,n2,g)).gt. &
4468 truncfacexp*max(1d0,maxcexpg(1,rg,g-1)) .or. &
4469 g.ge.2.and.abs(cexpg(n0,n1,n2,g)).gt. &
4470 truncfacexp*maxcexpg(1,rg,g-1))
then
4473 write(*,*)
'CalcCgn exit gloop',n0,n1,n2,g,abs(cexpg(n0,n1,n2,g)),maxcexpg(1,rg,g-1),truncfacexp
4486 if (rg-n0.le.rmax)
then
4489 c(n0,n1,n2) = c(n0,n1,n2) + cexpg(n0,n1,n2,g)
4502 maxcexpg(0,rg-1,g) = 0d0
4508 smod(1) = smod(1) - 2d0*n1*cexpg(1,n1-1,n2,g)
4511 smod(2) = smod(2) - 2d0*n2*cexpg(1,n1,n2-1,g)
4517 cexpg(0,n1,n2,g) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
4518 -
detz*cexpg(0,inds(1),inds(2),g-1))/zadjfj
4520 maxcexpg(0,rg-1,g) = maxcexpg(0,rg-1,g) + abs(cexpg(0,n1,n2,g))
4530 if (g.eq.1.and.abs(cexpg(0,n1,n2,g)).gt. &
4531 truncfacexp*max(1/
m2max,maxcexpg(0,rg-1,g-1)) .or. &
4532 g.ge.2.and.abs(cexpg(0,n1,n2,g)).gt. &
4533 truncfacexp*maxcexpg(0,rg-1,g-1))
then
4536 write(*,*)
'CalcCgn exit gloop',0,n1,n2,g,abs(cexpg(0,n1,n2,g)),maxcexpg(0,rg-1,g-1),truncfacexp
4546 c00_err(rg) = max(c00_err(rg), &
4547 max( abs(m02)*cij_err(rg-2), &
4548 max(
adetz*cij_err(rg),
fmax**2*cij_err(rg-2),
fmax*c00_err(rg-1))/abs(zadjkl) ) &
4551 cij_err(rg-1) = max(cij_err(rg-1),max(2*
maxzadj*c00_err(rg),
adetz*cij_err(rg))/abs(zadjfj) )
4554 c00_err2(rg) = max(c00_err2(rg), &
4555 max( abs(m02)*cij_err2(rg-2), &
4556 max(
adetz*cij_err2(rg),
fmax**2*cij_err2(rg-2),
fmax*c00_err2(rg-1))/abs(zadjkl) ) &
4559 cij_err2(rg-1) = max(cij_err2(rg-1),max(2*
maxzadj*c00_err2(rg),
adetz*cij_err2(rg))/abs(zadjfj) )
4566 if (rg-n0.le.rmax)
then
4569 c(n0,n1,n2) = c(n0,n1,n2) + cexpg(n0,n1,n2,g)
4576 if ((rg.le.rmax+1))
then
4580 c(0,n1,n2) = c(0,n1,n2) + cexpg(0,n1,n2,g)
4581 if(abs(cexpg(0,n1,n2,g-1)).ne.0d0)
then
4583 cerr(rg-1)=max(cerr(rg-1),abs(cexpg(0,n1,n2,g))*min(1d0,abs(cexpg(0,n1,n2,g))/abs(cexpg(0,n1,n2,g-1))))
4585 cerr(rg-1)=max(cerr(rg-1),abs(cexpg(0,n1,n2,g)))
4594 if(cij_err(rg-1).gt.cerr(rg-1))
then
4595 gtrunc = min(g,gtrunc)
4598 write(*,*)
'CalcCgn exit err',r,g,gtrunc &
4599 ,cij_err(rg-1),cerr(rg-1)
4608 write(*,*)
'CalcCgn C(0,0,0) = ',r,c(0,0,0)
4610 write(*,*)
'CalcCgn C(1,0,0) = ',r,c(1,0,0)
4611 write(*,*)
'CalcCgn C(0,1,0) = ',r,c(0,1,0)
4612 write(*,*)
'CalcCgn C(0,0,1) = ',r,c(0,0,1)
4614 if(r.gt.2.and.rmax.ge.2)
then
4615 write(*,*)
'CalcCgn C(0,2,0) = ',r,c(0,2,0)
4617 write(*,*)
'CalcCgn C(0,0,2) = ',r,c(0,0,2)
4619 if(r.gt.3.and.rmax.ge.3)
then
4620 write(*,*)
'CalcCgn C(1,0,1) = ',r,c(1,0,1)
4621 write(*,*)
'CalcCgn C(1,1,0) = ',r,c(1,1,0)
4622 write(*,*)
'CalcCgn C(1,0,1) = ',r,c(1,0,1)
4624 write(*,*)
'CalcCgn C(0,3,0) = ',r,c(0,3,0)
4625 write(*,*)
'CalcCgn C(0,2,1) = ',r,c(0,2,1)
4626 write(*,*)
'CalcCgn C(0,0,3) = ',r,c(0,0,3)
4628 write(*,*)
'CalcCgn Cij_err',r,cij_err
4629 write(*,*)
'CalcCgn Cij_acc',r,cij_err/abs(c(0,0,0))
4631 write(*,*)
'CalcCgn err',r,cerr
4632 write(*,*)
'CalcCgn acc',r,cerr/abs(c(0,0,0))
4635 cerr2 = max(cerr,cij_err2(0:rmax))
4636 cerr = max(cerr,cij_err(0:rmax))
4645 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0)
then
4652 c(n0,n1,rg-2*n0-n1)=0d0
4662 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0.and.r.gt.rmax)
then
4674 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
4675 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
4687 write(*,*)
'CalcCgn final err',cerr
4688 write(*,*)
'CalcCgn final acc',cerr/abs(c(0,0,0))
4693 write(*,*)
'CalcCgn rmax',rmax
4697 write(*,*)
'CalcCgn out',r,n0,n1,r-2*n0-n1,c(n0,n1,r-2*n0-n1)
4716 subroutine calccg(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,ordg_min,ordg_max,id,Cerr,acc_req_Cr,Cerr2)
4720 integer,
intent(in) :: rmax,ordg_min,ordg_max,id
4721 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
4722 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
4723 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
4724 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
4725 double precision,
intent(in) :: acc_req_Cr(0:rmax)
4726 double complex :: Xtilde,Zkl,Zadjfj
4727 double complex,
allocatable :: Cexpg(:,:,:,:), CuvExpg(:,:,:)
4728 double complex,
allocatable :: B_0(:,:,:), B_i(:,:,:), Shat(:,:,:,:)
4729 double complex,
allocatable :: Buv_0(:,:,:), Buv_i(:,:,:)
4730 double complex :: Smod(2), Skl
4731 double complex :: C0_coli, elimminf2_coli
4732 double precision,
allocatable :: C00_err(:),Cij_err(:)
4733 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
4734 double precision :: B_err,B_max
4735 double precision :: maxCexpg(0:1,0:rmax+ordg_min+1,0:ordg_max),truncfacexp
4736 integer :: rmaxB,rmaxExp,gtrunc,r,n0,n1,n2,k,l,j,sgn,g,rg,mr
4737 integer :: inds0(2), inds(2), ktlt(2)
4738 integer :: bin,nid(0:2)
4741 write(*,*)
'CalcCg in ',rmax,ordg_min,ordg_max,id
4744 write(*,*)
'CalcCg in ',rmax,ordg_min,ordg_max,id
4750 rmaxb = rmax + ordg_min
4751 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
4752 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
4753 allocate(b_i(0:rmaxb,0:rmaxb,2))
4754 allocate(buv_i(0:rmaxb,0:rmaxb,2))
4760 if (mod(id/bin,2).eq.0)
then
4767 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
4768 call calcb(b_i(:,:,1),buv_i(:,:,1),p20,m02,m22,rmaxb,nid(1))
4769 call calcb(b_i(:,:,2),buv_i(:,:,2),p10,m02,m12,rmaxb,nid(2))
4776 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
4777 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
4778 b_max = max(b_max,abs(b_0(0,n1,n2)))
4782 b_max = max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
4783 b_err = acc_def_b*b_max
4820 allocate(shat(0:rmaxb,0:rmaxb,0:rmaxb,2))
4827 shat(n0,n1,n2,:) = -b_0(n0,n1,n2)
4831 shat(n0,0,k,1) = shat(n0,0,k,1) + b_i(n0,k,1)
4832 shat(n0,k,0,2) = shat(n0,k,0,2) + b_i(n0,k,2)
4845 if (abs(
z(1,1)).ge.abs(
z(2,2)))
then
4846 if (abs(
z(1,1)).ge.abs(
z(1,2)))
then
4858 if (abs(
z(2,2)).ge.abs(
z(1,2)))
then
4874 xtilde =
xadj(3-k,3-l)
4876 xtilde = -
xadj(3-k,3-l)
4885 allocate(cexpg(0:rmaxexp/2,0:rmaxexp-1,0:rmaxexp-1,0:ordg_max))
4889 allocate(cuvexpg(0:rmaxexp,0:rmaxexp,0:rmaxexp))
4891 cuv(0:rmax,0:rmax,0:rmax) = cuvexpg(0:rmax,0:rmax,0:rmax)
4894 allocate(c00_err(0:rmaxexp))
4895 allocate(cij_err(0:rmaxexp))
4896 allocate(c00_err2(0:rmaxexp))
4897 allocate(cij_err2(0:rmaxexp))
4915 rloop:
do r=1,rmaxexp
4917 if (r.gt.rmax+gtrunc+1)
exit rloop
4937 if (inds(k).ge.1)
then
4939 skl = skl - 2d0*
f(l)*inds0(k)*cexpg(n0,inds(1),inds(2),0)
4940 if (inds(l).ge.1)
then
4942 skl = skl - 4d0*inds0(k)*(inds(l)+1)*cexpg(n0+1,inds(1),inds(2),0)
4946 if (inds(l).ge.1)
then
4948 skl = skl + 2d0*inds0(l)*shat(n0,inds(1),inds(2),k) &
4949 - 2d0*
f(k)*inds0(l)*cexpg(n0,inds(1),inds(2),0)
4952 cexpg(n0,n1,n2,0) = (2d0*zkl*b_0(n0-1,n1,n2) + xtilde*cexpg(n0-1,n1,n2,0) &
4953 -
z(1,k)*shat(n0-1,n1+1,n2,l) -
z(2,k)*shat(n0-1,n1,n2+1,l) &
4954 +
f(l)*shat(n0-1,n1,n2,k) + 4d0*zkl*cuvexpg(n0,n1,n2) + skl) &
4955 /(2d0*zkl)/(2d0*(r-n0)+1d0)
4958 maxcexpg(1,r,0) = maxcexpg(1,r,0) + abs(cexpg(n0,n1,n2,0))
4961 if (r-n0.le.rmax)
then
4962 c(n0,n1,n2) = cexpg(n0,n1,n2,0)
4971 maxcexpg(0,r-1,0)=0d0
4975 smod = shat(0,n1,n2,:)
4977 smod(1) = smod(1) - 2d0*n1*cexpg(1,n1-1,n2,0)
4980 smod(2) = smod(2) - 2d0*n2*cexpg(1,n1,n2-1,0)
4983 cexpg(0,n1,n2,0) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2))/zadjfj
4985 maxcexpg(0,r-1,0) = maxcexpg(0,r-1,0) + abs(cexpg(0,n1,n2,0))
4986 if (r-n0.le.rmax+1)
then
4987 c(0,n1,n2) = cexpg(0,n1,n2,0)
4992 if(r.le.rmax+1)
then
4994 cerr(r-1) =
fac_g*maxcexpg(0,r-1,0)
4998 c00_err(r) = max(max(
maxzadj*b_err,
fmax*b_err)/abs(zkl),b_err) &
5000 cij_err(r-1)=
maxzadj*max(b_err,2*c00_err(r))/abs(zadjfj)
5002 c00_err2(r) = max(max(
maxzadj*b_err,
fmax*b_err)/abs(zkl),b_err) &
5004 cij_err2(r-1)=
maxzadj*max(b_err,2*c00_err2(r))/abs(zadjfj)
5015 gloop:
do g=1,min(gtrunc,r-1)
5022 maxcexpg(1,rg,g) = 0d0
5031 if (inds(k).ge.1)
then
5033 skl = skl - 2d0*
f(l)*inds0(k)*cexpg(n0,inds(1),inds(2),g)
5034 if (inds(l).ge.1)
then
5036 skl = skl - 4d0*inds0(k)*(inds(l)+1)*cexpg(n0+1,inds(1),inds(2),g)
5040 if (inds(l).ge.1)
then
5042 skl = skl - 2d0*
f(k)*inds0(l)*cexpg(n0,inds(1),inds(2),g)
5046 cexpg(n0,n1,n2,g) = (xtilde*cexpg(n0-1,n1,n2,g) + skl &
5047 -
detz*sgn*cexpg(n0-1,inds(1),inds(2),g-1)) &
5048 /(2d0*zkl)/(2d0*(rg-n0)+1d0)
5050 maxcexpg(1,rg,g) = maxcexpg(1,rg,g) + abs(cexpg(n0,n1,n2,g))
5052 if (g.eq.1.and.abs(cexpg(n0,n1,n2,g)).gt. &
5053 truncfacexp*max(1d0,maxcexpg(1,rg,g-1)).or. &
5054 g.ge.2.and.abs(cexpg(n0,n1,n2,g)).gt. &
5055 truncfacexp*maxcexpg(1,rg,g-1))
then
5058 write(*,*)
'CalcCg exit gloop',n0,n1,n2,g,abs(cexpg(n0,n1,n2,g)),maxcexpg(1,rg,g-1)
5073 if (rg-n0.le.rmax)
then
5076 c(n0,n1,n2) = c(n0,n1,n2) + cexpg(n0,n1,n2,g)
5085 maxcexpg(0,rg-1,g) = 0d0
5091 smod(1) = smod(1) - 2d0*n1*cexpg(1,n1-1,n2,g)
5094 smod(2) = smod(2) - 2d0*n2*cexpg(1,n1,n2-1,g)
5101 cexpg(0,n1,n2,g) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
5102 -
detz*cexpg(0,inds(1),inds(2),g-1))/zadjfj
5104 maxcexpg(0,rg-1,g) = maxcexpg(0,rg-1,g) + abs(cexpg(0,n1,n2,g))
5106 if (g.eq.1.and.abs(cexpg(0,n1,n2,g)).gt. &
5107 truncfacexp*max(1/
m2max,maxcexpg(0,rg-1,g-1)) .or. &
5108 g.ge.2.and.abs(cexpg(0,n1,n2,g)).gt. &
5109 truncfacexp*maxcexpg(0,rg-1,g-1))
then
5112 write(*,*)
'CalcCg exit gloop',0,n1,n2,g,abs(cexpg(0,n1,n2,g)),maxcexpg(0,rg-1,g-1)
5113 write(*,*)
'CalcCg exit gloop',abs(cexpg(0,n1,n2,g)).gt.truncfacexp*maxcexpg(0,rg-1,g-1),truncfacexp
5114 write(*,*)
'CalcCg exit gloop',
zadj(1,j)*smod(1)/zadjfj ,
zadj(2,j)*smod(2)/zadjfj, &
5115 -
detz*cexpg(0,inds(1),inds(2),g-1)/zadjfj
5121 write(*,*)
'CalcCg exit gloop',rmax,g,rmaxexp
5142 c00_err(rg) = max(c00_err(rg), &
5143 max(
adetz*cij_err(rg),abs(xtilde)*cij_err(rg-2),
fmax*c00_err(rg-1))/abs(zkl) &
5146 cij_err(rg-1) = max(cij_err(rg-1),max(2*
maxzadj*c00_err(rg),
adetz*cij_err(rg))/abs(zadjfj) )
5149 c00_err2(rg) = max(c00_err2(rg), &
5150 max(
adetz*cij_err2(rg),abs(xtilde)*cij_err2(rg-2),
fmax*c00_err2(rg-1))/abs(zkl) &
5153 cij_err2(rg-1) = max(cij_err2(rg-1),max(2*
maxzadj*c00_err2(rg),
adetz*cij_err2(rg))/abs(zadjfj) )
5161 if (rg-n0.le.rmax)
then
5164 c(n0,n1,n2) = c(n0,n1,n2) + cexpg(n0,n1,n2,g)
5170 if ((rg.le.rmax+1))
then
5174 c(0,n1,n2) = c(0,n1,n2) + cexpg(0,n1,n2,g)
5175 if(abs(cexpg(0,n1,n2,g-1)).ne.0d0)
then
5177 cerr(rg-1)=max(cerr(rg-1),abs(cexpg(0,n1,n2,g))*min(1d0,abs(cexpg(0,n1,n2,g))/abs(cexpg(0,n1,n2,g-1))))
5179 cerr(rg-1)=max(cerr(rg-1),abs(cexpg(0,n1,n2,g)))
5188 if(cij_err(rg-1).gt.cerr(rg-1))
then
5189 gtrunc = min(g,gtrunc)
5192 write(*,*)
'CalcCg exit err',r,g,gtrunc
5202 write(*,*)
'CalcCg C(0,0,0) = ',r,c(0,0,0)
5203 write(*,*)
'CalcCg C(2,0,0) = ',r,c(1,0,0)
5204 write(*,*)
'CalcCg C(0,1,0) = ',r,c(0,1,0)
5205 write(*,*)
'CalcCg C(0,0,1) = ',r,c(0,0,1)
5206 if(r.ge.5.and.rmax.ge.5)
then
5207 write(*,*)
'CalcCg C(2,1,0) = ',r,c(2,1,0)
5212 write(*,*)
'CalcCg Cerr r =',r,cerr
5213 write(*,*)
'CalcCg Cij_err =',r,cij_err
5216 cerr2 = max(cerr,cij_err2(0:rmax))
5217 cerr = max(cerr,cij_err(0:rmax))
5220 write(*,*)
'CalcCg Cerr =',r,cerr,maxval(cerr)
5247 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0)
then
5254 c(n0,n1,rg-2*n0-n1)=0d0
5264 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0.and.r.gt.rmax)
then
5278 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
5279 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
5303 write(*,*)
'CalcCg final err',cerr
5304 write(*,*)
'CalcCg final acc',cerr/abs(c(0,0,0))
5308 write(*,*)
'CalcCg rmax',rmax
5312 write(*,*)
'CalcCg out ',r,n0,n1,r-2*n0-n1
5313 write(*,*)
'CalcCg out C',c(n0,n1,r-2*n0-n1)
5330 subroutine calccgr(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,ordgr_min,ordgr_max,id,Cerr,acc_req_Cr,Cerr2)
5334 integer,
intent(in) :: rmax,ordgr_min,ordgr_max,id
5335 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
5336 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
5337 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
5338 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
5339 double precision,
intent(in) :: acc_req_Cr(0:rmax)
5340 double complex,
allocatable :: B_0(:,:,:), B_i(:,:,:), Shat(:,:,:,:)
5341 double complex,
allocatable :: Buv_0(:,:,:), Buv_i(:,:,:)
5342 double precision :: B_err,B_max
5343 double complex :: Zadjfj,Zadj2(2,2), Zadjkl, Zadj2f(2,2,2)
5344 double complex,
allocatable :: Cexpgr(:,:,:,:), CuvExpgr(:,:,:)
5345 double complex :: Smod(2), Skl, Caux
5346 double complex :: elimminf2_coli
5347 double precision,
allocatable :: C00_err(:),Cij_err(:)
5348 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
5349 double precision :: maxZadj2f
5350 double precision :: maxCexpgr(0:1,0:2*(rmax+ordgr_min),0:ordgr_max),truncfacexp
5351 integer :: rmaxB,rmaxExp,gtrunc,r,n0,n1,n2,k,l,i,j,m,n,g,rg,lt,ltt,nn,nntt
5352 integer :: inds0(2), inds1(2), inds(2)
5353 integer :: bin,nid(0:2)
5356 write(*,*)
'CalcCgr in ',rmax,ordgr_min,ordgr_max
5357 write(*,*)
'CalcCgr in, f ',
f
5360 write(*,*)
'CalcCgr in ',rmax,ordgr_min,ordgr_max
5364 rmaxb = 2*rmax + 2*ordgr_min
5365 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
5366 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
5367 allocate(b_i(0:rmaxb,0:rmaxb,2))
5368 allocate(buv_i(0:rmaxb,0:rmaxb,2))
5374 if (mod(id/bin,2).eq.0)
then
5381 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
5382 call calcb(b_i(:,:,1),buv_i(:,:,1),p20,m02,m22,rmaxb,nid(1))
5383 call calcb(b_i(:,:,2),buv_i(:,:,2),p10,m02,m12,rmaxb,nid(2))
5390 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
5391 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
5394 b_max = max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
5395 b_err = acc_def_b*b_max
5445 allocate(shat(0:rmaxb,0:rmaxb,0:rmaxb,2))
5452 shat(n0,n1,n2,:) = -b_0(n0,n1,n2)
5456 shat(n0,0,k,1) = shat(n0,0,k,1) + b_i(n0,k,1)
5457 shat(n0,k,0,2) = shat(n0,k,0,2) + b_i(n0,k,2)
5470 zadj2f(1,2,1) = -
f(2)
5471 zadj2f(1,2,2) =
f(1)
5476 if (abs(zadj2f(1,2,1)).gt.maxzadj2f)
then
5477 maxzadj2f = abs(zadj2f(1,2,1))
5484 if (abs(zadj2f(1,2,2)).gt.maxzadj2f)
then
5485 maxzadj2f = abs(zadj2f(1,2,2))
5494 write(*,*)
'CalcCgr maxZadj2f ',maxzadj2f,maxval(abs(zadj2f(:,:,:)))
5495 write(*,*)
'CalcCgr Zadj2f ',zadj2f
5496 write(*,*)
'CalcCgr Zadj2f ',zadj2f(1,2,1),zadj2f(1,2,1)
5497 write(*,*)
'CalcCgr f ',
f
5504 write(*,*)
'CalcCgr k,n,nt,l',k,n,l,m
5505 write(*,*)
'CalcCgr pars', maxzadj2f,zadj2f(k,n,l),
zadj(k,l),
maxzadj
5506 write(*,*)
'CalcCgr pars', abs(
zadjf(l)),abs(
detz)
5507 write(*,*)
'CalcCgr pars', abs(
zadjf(l)/ maxzadj2f),abs(
detz/maxzadj2f)
5515 allocate(cexpgr(0:rmaxexp/2,0:rmaxexp,0:rmaxexp,0:ordgr_max))
5518 allocate(cuvexpgr(0:(rmaxexp+1),0:rmaxexp+1,0:rmaxexp+1))
5520 cuv(0:rmax,0:rmax,0:rmax) = cuvexpgr(0:rmax,0:rmax,0:rmax)
5523 allocate(c00_err(0:rmaxexp))
5524 allocate(cij_err(0:rmaxexp))
5525 allocate(c00_err2(0:rmaxexp))
5526 allocate(cij_err2(0:rmaxexp))
5546 rloop:
do r=0,rmaxexp/2
5552 if (r.gt.rmax+gtrunc)
exit rloop
5555 write(*,*)
'CalcCgr rloop',r,rmaxexp,rmaxb
5565 maxcexpgr(1,r,0)=0d0
5571 write(*,*)
'CalcCgr rloop',n0,nn,zadj2f(k,n,l)
5578 write(*,*)
'CalcCgr inds0',n0,inds0
5585 write(*,*)
'CalcCgr inds1',n0,inds1
5588 caux = -
zadj(k,l)*b_0(n0-1,inds1(1),inds1(2))
5602 write(*,*)
'CalcCgr Caux 1c',-
zadj(k,l)*b_0(n0-1,inds1(1),inds1(2))
5603 write(*,*)
'CalcCgr Caux 1s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5609 caux = caux -
zadj(i,l)*shat(n0-1,inds(1),inds(2),i)
5611 write(*,*)
'CalcCgr Caux 2ci', -
zadj(i,l)*shat(n0-1,inds(1),inds(2),i)
5616 write(*,*)
'CalcCgr Caux 2s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5622 caux = caux +
zadj(k,l)*shat(n0-1,inds(1),inds(2),i)
5624 write(*,*)
'CalcCgr Caux 3ci',
zadj(k,l)*shat(n0-1,inds(1),inds(2),i)
5630 write(*,*)
'CalcCgr Caux 3s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5633 caux = caux + 2*(nn+1) *zadj2(n ,m )*shat(n0,inds0(1),inds0(2),m)
5637 write(*,*)
'CalcCgr Caux 4ca', 2*(nn+1) *zadj2(n ,m )*shat(n0,inds0(1),inds0(2),m)
5638 write(*,*)
'CalcCgr Caux 4s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5645 if (inds(n).gt.1)
then
5647 caux = caux - 4*(nn+1)*nn * zadj2(n,m ) * cexpgr(n0+1,inds(1),inds(2),0)
5649 write(*,*)
'CalcCgr Caux 6c',4*(nn+1)*nn* zadj2(n,m ) *cexpgr(n0+1,inds(1),inds(2),0)
5650 write(*,*)
'CalcCgr Caux 6s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5654 if (inds(n).gt.0.and.inds(m).gt.0)
then
5657 caux = caux - 4*(nn+1)*(inds(m)+1)* zadj2(n,m ) * cexpgr(n0+1,inds(1),inds(2),0)
5659 write(*,*)
'CalcCgr Caux 6c',-4*(nn+1)*(inds(m)+1)* zadj2(n,m ) *cexpgr(n0+1,inds(1),inds(2),0)
5660 write(*,*)
'CalcCgr Caux 6s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5665 cexpgr(n0,inds0(1),inds0(2),0) = caux/(2*(nn+1)* zadj2f(k,n,l))
5668 maxcexpgr(1,r,0) = maxcexpgr(1,r,0) + abs(cexpgr(n0,inds0(1),inds0(2),0) )
5673 c(n0,inds0(1),inds0(2)) = cexpgr(n0,inds0(1),inds0(2),0)
5682 maxcexpgr(0,r,0)=0d0
5686 smod = shat(0,n1,n2,:)
5688 smod(1) = smod(1) - 2d0*n1*cexpgr(1,n1-1,n2,0)
5691 smod(2) = smod(2) - 2d0*n2*cexpgr(1,n1,n2-1,0)
5694 cexpgr(0,n1,n2,0) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
5696 maxcexpgr(0,r,0) = maxcexpgr(0,r,0) + abs(cexpgr(0,n1,n2,0))
5698 c(0,n1,n2) = cexpgr(0,n1,n2,0)
5704 write(*,*)
'CalcCgr C(0,n1,n2,0)=',n1,n2,c(0,n1,n2)
5707 if(n0.eq.0.and.n1.eq.0.and.n2.eq.3)
then
5708 write(*,*)
'C(0,0,3)= ',0,c(n0,n1,n2)
5720 cerr(r) =
fac_gr*maxcexpgr(0,r,0)
5725 c00_err(r+1) =
maxzadj*b_err/(2*maxzadj2f)
5727 cij_err(r)=
maxzadj*max(b_err,2*c00_err(r+1))/abs(zadjfj)
5730 c00_err2(r+1) =
maxzadj*b_err/(2*maxzadj2f)
5732 cij_err2(r)=
maxzadj*max(b_err,2*c00_err2(r+1))/abs(zadjfj)
5739 gloop:
do g=1,min(gtrunc,r)
5743 write(*,*)
'CalcCgr: gloop ',r,rg,g
5749 maxcexpgr(1,rg,g) = 0d0
5760 write(*,*)
'CalcCgr Caux r inds=',n0,inds0
5763 caux = 2*
zadj(k,l) * (2+rg-n0) * cexpgr(n0,inds1(1),inds1(2),g-1)
5766 write(*,*)
'CalcCgr Caux r1c',2*
zadj(k,l)*(2+rg-n0)* cexpgr(n0,inds1(1),inds1(2),g-1)
5767 write(*,*)
'CalcCgr Caux r1c',2*
zadj(k,l)*(2+rg-n0),cexpgr(n0,inds1(1),inds1(2),g-1) &
5768 ,n0,inds1(1),inds1(2)
5769 write(*,*)
'CalcCgr Caux r1s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5774 inds(k) = inds(k) + 1
5775 inds(l) = inds(l) + 1
5776 caux = caux +
detz * cexpgr(n0-1,inds(1),inds(2),g-2)
5779 write(*,*)
'CalcCgr Caux r2c',
detz * cexpgr(n0-1,inds(1),inds(2),g-2)
5780 write(*,*)
'CalcCgr Caux r2s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5785 inds(k) = inds(k) + 1
5786 caux = caux +
zadjf(l) * cexpgr(n0-1,inds(1),inds(2),g-1)
5789 write(*,*)
'CalcCgr Caux r3c',
zadjf(l)* cexpgr(n0-1,inds(1),inds(2),g-1)
5790 write(*,*)
'CalcCgr Caux r3c',
zadjf(l),cexpgr(n0-1,inds(1),inds(2),g-1),n0-1,inds(1),inds(2)
5791 write(*,*)
'CalcCgr Caux r3s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5798 if (inds(n).gt.1)
then
5800 caux = caux - 4*(nn+1)*nn * zadj2(n,m ) * cexpgr(n0+1,inds(1),inds(2),g)
5802 write(*,*)
'CalcCgr Caux r6c',4*(nn+1)*nn* zadj2(n,m ) *cexpgr(n0+1,inds(1),inds(2),g)
5803 write(*,*)
'CalcCgr Caux r6s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5807 if (inds(n).gt.0.and.inds(m).gt.0)
then
5810 caux = caux - 4*(nn+1)*(inds(m)+1)* zadj2(n,m ) * cexpgr(n0+1,inds(1),inds(2),g)
5812 write(*,*)
'CalcCgr Caux r6c',4*(nn+1)*(inds(m)+1)* zadj2(n,m ) *cexpgr(n0+1,inds(1),inds(2),g)
5813 write(*,*)
'CalcCgr Caux r6c',n,m,nn,4*(nn+1)*(inds(m)+1),zadj2(n,m ),cexpgr(n0+1,inds(1),inds(2),g)
5814 write(*,*)
'CalcCgr Caux r6s',caux,caux/(2*(nn+1)* zadj2f(k,n,l))
5820 cexpgr(n0,inds0(1),inds0(2),g) = caux/(2*(nn+1)* zadj2f(k,n,l))
5823 maxcexpgr(1,rg,g) = maxcexpgr(1,rg,g) + abs(cexpgr(n0,inds0(1),inds0(2),g))
5825 if (g.eq.1.and.abs(cexpgr(n0,inds0(1),inds0(2),g)).gt. &
5826 truncfacexp*max(1d0,maxcexpgr(1,rg,g-1)) .or. &
5827 g.ge.2.and.abs(cexpgr(n0,inds0(1),inds0(2),g)).gt. &
5828 truncfacexp*maxcexpgr(1,rg,g-1))
then
5831 write(*,*)
'CalcCgr exit gloop',n0,inds0(1),inds0(2),g,rg, &
5832 abs(cexpgr(n0,inds0(1),inds0(2),g)),maxcexpgr(1,rg,g-1),truncfacexp
5844 if (rg.le.rmax)
then
5847 if (rg.le.rmax)
then
5850 c(n0,n1,n2) = c(n0,n1,n2) + cexpgr(n0,n1,n2,g)
5864 maxcexpgr(0,rg,g) = 0d0
5870 smod(1) = smod(1) - 2d0*n1*cexpgr(1,n1-1,n2,g)
5873 smod(2) = smod(2) - 2d0*n2*cexpgr(1,n1,n2-1,g)
5879 cexpgr(0,n1,n2,g) = (
zadj(1,j)*smod(1) +
zadj(2,j)*smod(2) &
5880 -
detz*cexpgr(0,inds(1),inds(2),g-1))/zadjfj
5882 maxcexpgr(0,rg,g) = maxcexpgr(0,rg,g) + abs(cexpgr(0,n1,n2,g))
5892 if (g.eq.1.and.abs(cexpgr(0,n1,n2,g)).gt. &
5893 truncfacexp*max(1d0/
m2scale,maxcexpgr(0,rg,g-1)).or. &
5894 g.ge.2.and.abs(cexpgr(0,n1,n2,g)).gt. &
5895 truncfacexp*maxcexpgr(0,rg,g-1))
then
5898 write(*,*)
'CalcCgr exit gloop',0,n1,n2,g,abs(cexpgr(0,n1,n2,g)),maxcexpgr(0,rg,g-1),truncfacexp
5908 c00_err(rg+1) = max( c00_err(rg+1), &
5909 max(
maxzadj*(2+rg-2*n0)*c00_err(rg+2), &
5910 abs(
detz)*cij_err(rg+2), &
5914 cij_err(rg)=max(cij_err(rg), &
5915 max(2*
maxzadj*c00_err(rg+1),abs(
detz)*cij_err(rg))/abs(zadjfj) )
5918 c00_err2(rg+1) = max( c00_err2(rg+1), &
5919 max(
maxzadj*(2+rg-2*n0)*c00_err2(rg+2), &
5920 abs(
detz)*cij_err2(rg+2), &
5924 cij_err2(rg)=max(cij_err2(rg), &
5925 max(2*
maxzadj*c00_err2(rg+1),abs(
detz)*cij_err2(rg))/abs(zadjfj) )
5928 if (rg.le.rmax)
then
5931 if (rg.le.rmax)
then
5934 c(n0,n1,n2) = c(n0,n1,n2) + cexpgr(n0,n1,n2,g)
5941 if (rg.le.rmax)
then
5945 c(0,n1,n2) = c(0,n1,n2) + cexpgr(0,n1,n2,g)
5946 if(abs(cexpgr(0,n1,n2,g-1)).ne.0d0)
then
5948 cerr(rg)=max(cerr(rg),abs(cexpgr(0,n1,n2,g))*min(1d0,abs(cexpgr(0,n1,n2,g))/abs(cexpgr(0,n1,n2,g-1))))
5950 cerr(rg)=max(cerr(rg),abs(cexpgr(0,n1,n2,g)))
5954 write(*,*)
'CalcCgr Cerr calc',rg,cerr(rg),n1,n2,abs(cexpgr(0,n1,n2,g)),abs(cexpgr(0,n1,n2,g-1))
5960 if(cij_err(rg).gt.3d0*cerr(rg))
then
5961 gtrunc = min(g,gtrunc)
5964 write(*,*)
'CalcCgr exit err',r,rg,g,gtrunc,cij_err(rg),cerr(rg)
5974 write(*,*)
'CalcCgr C(0,0,0) = ',r,c(0,0,0)
5976 write(*,*)
'CalcCgr C(1,0,0) = ',r,c(1,0,0)
5977 write(*,*)
'CalcCgr C(0,1,0) = ',r,c(0,1,0)
5978 write(*,*)
'CalcCgr C(0,0,1) = ',r,c(0,0,1)
5979 write(*,*)
'CalcCgr C(0,0,0) = ',r,c(0,0,0)
5981 if(r.ge.2.and.rmax.ge.2)
then
5982 write(*,*)
'CalcCgr C(1,1,0) = ',r,c(1,1,0)
5983 write(*,*)
'CalcCgr C(1,0,1) = ',r,c(1,0,1)
5984 write(*,*)
'CalcCgr C(1,0,0) = ',r,c(1,0,0)
5987 write(*,*)
'CalcCgr C(0,0,2) = ',r,c(0,0,2)
5989 if(r.ge.3.and.rmax.ge.2)
then
5992 write(*,*)
'CalcCgr C(1,0,2) = ',r,c(1,0,2)
5993 write(*,*)
'CalcCgr C(0,3,0) = ',r,c(0,3,0)
5994 write(*,*)
'CalcCgr C(0,2,1) = ',r,c(0,2,1)
5995 write(*,*)
'CalcCgr C(0,0,3) = ',r,c(0,0,3)
5996 write(*,*)
'CalcCgr C(0,1,1) = ',r,c(0,1,1)
5997 write(*,*)
'CalcCgr C(0,0,2) = ',r,c(0,0,2)
5999 write(*,*)
'CalcCgr Cij_err',r,cij_err
6000 write(*,*)
'CalcCgr Cij_acc',r,cij_err/abs(c(0,0,0))
6002 write(*,*)
'CalcCgr err',r,cerr
6003 write(*,*)
'CalcCgr acc',r,cerr/abs(c(0,0,0))
6006 cerr2 = max(cerr,cij_err2(0:rmax))
6007 cerr = max(cerr,cij_err(0:rmax))
6016 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0)
then
6020 c(n0,n1,rg-n0-n1)=0d0
6025 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0.and.r.ge.rmax)
then
6038 write(*,*)
'CalcCgr final err',cerr
6039 write(*,*)
'CalcCgr final acc',cerr /abs(c(0,0,0))
6044 write(*,*)
'CalcCgr rmax',rmax
6048 write(*,*)
'CalcCgr out',r,n0,n1,r-2*n0-n1,c(n0,n1,r-2*n0-n1)
6069 subroutine calccgy(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,ordgy_min,ordgy_max,id,Cerr,acc_req_Cr,Cerr2)
6073 integer,
intent(in) :: rmax,ordgy_min,ordgy_max,id
6074 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
6075 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
6076 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
6077 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
6078 double precision,
intent(in) :: acc_req_Cr(0:rmax)
6079 double complex,
allocatable :: Cexpgy(:,:,:,:), CuvExpgy(:,:,:)
6080 double complex,
allocatable :: B_0(:,:,:), B_i(:,:,:), Shat(:,:,:,:)
6081 double complex,
allocatable :: Buv_0(:,:,:), Buv_i(:,:,:)
6082 double complex :: Smod, Caux, Zadj2f
6083 double complex :: C0_coli, elimminf2_coli
6084 double precision,
allocatable :: C00_err(:),Cij_err(:)
6085 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
6086 double precision :: B_err,B_max,aZadj2f
6087 double precision :: maxCexpgy(0:1,0:rmax+2*ordgy_min,0:ordgy_max),truncfacexp
6088 integer :: rmaxB,rmaxExp,gtrunc,r,n0,n1,n2,i,j,jt,g,rg
6089 integer :: inds0(2),inds(2),k,l,lt,nl,nlt
6090 integer :: bin,nid(0:2)
6093 write(*,*)
'CalcCgy in ',rmax,ordgy_min,ordgy_max,id
6094 write(*,*)
'CalcCgy in ',p10,p21,p20,m02,m12,m22
6097 write(*,*)
'CalcCgy in ',rmax,ordgy_min,ordgy_max,id
6102 rmaxb = rmax + 2*ordgy_min + 1
6103 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
6104 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
6105 allocate(b_i(0:rmaxb,0:rmaxb,2))
6106 allocate(buv_i(0:rmaxb,0:rmaxb,2))
6112 if (mod(id/bin,2).eq.0)
then
6119 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
6120 call calcb(b_i(:,:,1),buv_i(:,:,1),p20,m02,m22,rmaxb,nid(1))
6121 call calcb(b_i(:,:,2),buv_i(:,:,2),p10,m02,m12,rmaxb,nid(2))
6128 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
6129 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
6132 b_max = max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
6133 b_err = acc_def_b*b_max
6146 if (abs(
detz).lt.abs(4d0*
q10*
q20 +
z(2,1)*
z(2,1))*1d-4)
then
6171 allocate(shat(0:rmaxb,0:rmaxb,0:rmaxb,2))
6178 shat(n0,n1,n2,:) = -b_0(n0,n1,n2)
6182 shat(n0,0,k,1) = shat(n0,0,k,1) + b_i(n0,k,1)
6183 shat(n0,k,0,2) = shat(n0,k,0,2) + b_i(n0,k,2)
6189 if (abs(
xadj(1,1)).ge.abs(
xadj(2,2)))
then
6190 if (abs(
xadj(1,1)).ge.abs(
xadj(1,2)))
then
6202 if (abs(
xadj(2,2)).ge.abs(
xadj(1,2)))
then
6214 azadj2f = abs(zadj2f)
6216 if (abs(
zadj(1,1)).ge.abs(
zadj(2,2)))
then
6217 if (abs(
zadj(1,1)).ge.abs(
zadj(1,2)))
then
6227 if (abs(
zadj(2,2)).ge.abs(
zadj(1,2)))
then
6239 write(*,*)
'CalcCgy: Zadj',k,l,
zadj(k,l)
6240 write(*,*)
'CalcCgy: Xadj',i,j,
xadj(i,j)
6248 allocate(cexpgy(0:max(rmax/2,1),0:rmaxexp-2,0:rmaxexp-2,0:ordgy_max))
6251 allocate(cuvexpgy(0:rmaxexp,0:rmaxexp,0:rmaxexp))
6253 cuv(0:rmax,0:rmax,0:rmax) = cuvexpgy(0:rmax,0:rmax,0:rmax)
6256 allocate(c00_err(0:rmaxexp))
6257 allocate(cij_err(0:rmaxexp))
6258 allocate(c00_err2(0:rmaxexp))
6259 allocate(cij_err2(0:rmaxexp))
6280 write(*,*)
'CalcCgy gtrunc orig=',gtrunc
6281 write(*,*)
'CalcCgy rmaxExp-2=',rmaxexp-2
6286 rloop:
do r=0,rmaxexp-2
6289 write(*,*)
'CalcCgy rloop=',r,rmaxexp-2,rmax+2*gtrunc+2
6290 write(*,*)
'CalcCgy rloop=',rmax,gtrunc
6293 if (r.gt.rmax+2*gtrunc+2)
exit rloop
6300 maxcexpgy(1,r,0)=0d0
6308 caux =
zadj(k,1)*shat(0,inds(1),inds(2),1) &
6309 +
zadj(k,2)*shat(0,inds(1),inds(2),2)
6313 caux = caux - 2*nlt*
zadj(k,lt)*cexpgy(1,inds(1),inds(2),0)
6316 cexpgy(1,inds0(1),inds0(2),0) = caux/(2*(nl+1)*
zadj(k,l))
6317 maxcexpgy(1,r,0) = maxcexpgy(1,r,0) + abs(cexpgy(1,inds0(1),inds0(2),0) )
6319 if (r+1.le.rmax)
then
6320 c(1,inds0(1),inds0(2)) = cexpgy(1,inds0(1),inds0(2),0)
6326 maxcexpgy(0,r,0)=0d0
6332 caux = (2*(2+r)*cexpgy(1,n1,n2,0) - 4*cuvexpgy(1,n1,n2) &
6333 - b_0(0,n1,n2))*
zadj(i,j)
6337 smod = shat(0,n1,n2,jt)
6339 if (inds(jt).ge.1)
then
6340 inds(jt) = inds(jt)-1
6341 smod = smod - 2d0*(inds(jt)+1)*cexpgy(1,inds(1),inds(2),0)
6344 caux = caux + zadj2f*smod
6348 cexpgy(0,n1,n2,0) = caux/
xadj(i,j)
6349 maxcexpgy(0,r,0) = maxcexpgy(0,r,0) + abs(cexpgy(0,n1,n2,0))
6351 c(0,n1,n2) = cexpgy(0,n1,n2,0)
6358 cerr(r) =
fac_gy*maxcexpgy(0,r,0)
6365 c00_err(r+2) = b_err /2d0
6366 cij_err(r)=max(abs(
zadj(i,j))/abs(
xadj(i,j))*max(b_err,2*(r+2)*c00_err(r+2)), &
6367 fmax/abs(
xadj(i,j))*max(b_err,2*c00_err(r+1)))
6369 c00_err2(r+2) = b_err /2d0
6370 cij_err2(r)=max(abs(
zadj(i,j))/abs(
xadj(i,j))*max(b_err,2*(r+2)*c00_err2(r+2)), &
6371 fmax/abs(
xadj(i,j))*max(b_err,2*c00_err2(r+1)))
6379 gloop:
do g=1,min(gtrunc,r/2)
6385 maxcexpgy(1,rg,g) = 0d0
6393 caux = -
zadjf(k)*cexpgy(0,inds(1),inds(2),g-1)
6396 caux = caux -
detz*cexpgy(0,inds(1),inds(2),g-1)
6401 caux = caux - 2*nlt*
zadj(k,lt)*cexpgy(1,inds(1),inds(2),g)
6404 cexpgy(1,inds0(1),inds0(2),g) = caux/(2*(nl+1)*
zadj(k,l))
6405 maxcexpgy(1,rg,g) = maxcexpgy(1,rg,g) + abs(cexpgy(1,inds0(1),inds0(2),g) )
6407 if (g.eq.1.and.abs(cexpgy(1,inds0(1),inds0(2),g)).gt. &
6408 truncfacexp*max(1d0,maxcexpgy(1,rg,g-1)) .or. &
6409 g.ge.2.and.abs(cexpgy(1,inds0(1),inds0(2),g)).gt. &
6410 truncfacexp*maxcexpgy(1,rg,g-1))
then
6412 write(*,*)
'CalcCgy exit gloop',n1,n2,g,abs(cexpgy(1,inds0(1),inds0(2),g)),maxcexpgy(1,rg,g-1)
6413 write(*,*)
'CalcCgy exit gloop',g,rg,inds0(1),inds0(2)
6425 if (rg+1.le.rmax)
then
6430 c(1,inds0(1),inds0(2)) = c(1,inds0(1),inds0(2)) &
6431 + cexpgy(1,inds0(1),inds0(2),g)
6437 maxcexpgy(0,rg,g) = 0d0
6443 caux = 2*(2+rg)*cexpgy(1,n1,n2,g)*
zadj(i,j)
6447 if (inds0(jt).ge.1)
then
6449 inds(jt) = inds(jt)-1
6450 caux = caux - 2d0*zadj2f*inds0(jt)*cexpgy(1,inds(1),inds(2),g)
6455 inds0(i) = inds0(i)+1
6456 caux = caux -
zadjf(j)*cexpgy(0,inds0(1),inds0(2),g-1)
6460 cexpgy(0,n1,n2,g) = caux/
xadj(i,j)
6464 maxcexpgy(0,rg,g) = maxcexpgy(0,rg,g) + abs(cexpgy(0,n1,n2,g))
6466 if (g.eq.1.and.abs(cexpgy(0,n1,n2,g)).gt. &
6467 truncfacexp*max(1d0/
m2scale,maxcexpgy(0,rg,g-1)).or. &
6468 g.ge.2.and.abs(cexpgy(0,n1,n2,g)).gt. &
6469 truncfacexp*maxcexpgy(0,rg,g-1))
then
6472 write(*,*)
'CalcCgy exit gloop',n1,n2,g,rg
6473 write(*,*)
'CalcCgy exit gloop',abs(cexpgy(0,n1,n2,g)),maxcexpgy(0,rg,g-1),1d0/
m2scale
6474 write(*,*)
'CalcCgy exit gloop',truncfacexp
6490 c00_err(rg+2) =max(c00_err(rg+2), &
6491 max(abs(
zadjf(k))/2d0*cij_err(rg+1), &
6492 abs(
detz)/2d0*cij_err(rg+2))/abs(
zadj(k,l)))
6496 write(*,*)
'CalcCgy test2',rg,i,j,cij_err(rg)
6497 write(*,*)
'CalcCgy test2',rg,cij_err(rg+1),c00_err(rg+1)
6498 write(*,*)
'CalcCgy test2',rg,abs(
zadj(i,j))
6499 write(*,*)
'CalcCgy test2',rg,abs(zadj2f)
6500 write(*,*)
'CalcCgy test2',rg,abs(
zadjf(j))
6501 write(*,*)
'CalcCgy test2',rg,abs(
xadj(i,j))
6504 cij_err(rg)= max( cij_err(rg), &
6505 max(2*(rg+2)*abs(
zadj(i,j))*c00_err(rg+2), &
6506 2*abs(zadj2f)*c00_err(rg+1), &
6507 abs(
zadjf(j))*cij_err(rg+1))/abs(
xadj(i,j)))
6510 c00_err2(rg+2) =max(c00_err2(rg+2), &
6511 max(abs(
zadjf(k))/2d0*cij_err2(rg+1), &
6512 abs(
detz)/2d0*cij_err2(rg+2))/abs(
zadj(k,l)))
6515 cij_err2(rg)= max( cij_err2(rg), &
6516 max(2*(rg+2)*abs(
zadj(i,j))*c00_err2(rg+2), &
6517 2*abs(zadj2f)*c00_err2(rg+1), &
6518 abs(
zadjf(j))*cij_err2(rg+1))/abs(
xadj(i,j)))
6522 if (rg+1.le.rmax)
then
6527 c(1,inds0(1),inds0(2)) = c(1,inds0(1),inds0(2)) &
6528 + cexpgy(1,inds0(1),inds0(2),g)
6533 if ((rg.le.rmax))
then
6537 c(0,n1,n2) = c(0,n1,n2) + cexpgy(0,n1,n2,g)
6540 write(*,*)
'CalcCgy test1',rg,n1,n2,cerr(rg)
6541 write(*,*)
'CalcCgy test1',cexpgy(0,n1,n2,g)
6542 write(*,*)
'CalcCgy test1',cexpgy(0,n1,n2,g-1)
6545 if(abs(cexpgy(0,n1,n2,g-1)).ne.0d0)
then
6546 cerr(rg)=max(cerr(rg),abs(cexpgy(0,n1,n2,g))*min(1d0,abs(cexpgy(0,n1,n2,g))/abs(cexpgy(0,n1,n2,g-1))))
6548 cerr(rg)=max(cerr(rg),abs(cexpgy(0,n1,n2,g)))
6552 write(*,*)
'CalcCgy test1',cerr(rg)
6558 if(cij_err(rg).gt.cerr(rg))
then
6559 gtrunc = min(g,gtrunc)
6563 write(*,*)
'CalcCgy adjust gtrunc',r,g,gtrunc
6576 write(*,*)
'CalcCgy Cerr r =',r
6577 write(*,*)
'CalcCgy Cerr r =',r,cerr
6578 write(*,*)
'CalcCgy Cacc r =',r,cerr/abs(c(0,0,0))
6579 write(*,*)
'CalcCgy Cij_err =',r,cij_err
6580 write(*,*)
'CalcCgy C(0,0,0)=',r,c(0,0,0)
6581 if(rmax.ge.1.and.r.ge.1)
then
6582 write(*,*)
'CalcCgy C(0,1,0)=',r,c(0,1,0)
6583 if(rmax.ge.2.and.r.ge.2)
then
6584 write(*,*)
'CalcCgy C(0,1,1)=',r,c(0,1,1)
6585 if(rmax.ge.3.and.r.ge.3)
then
6586 write(*,*)
'CalcCgy C(0,1,2)=',r,c(0,1,2)
6588 write(*,*)
'CalcCgy C(0,0,4)=',r,c(0,0,4)
6595 cerr2 = max(cerr,cij_err2(0:rmax))
6596 cerr = max(cerr,cij_err(0:rmax))
6599 write(*,*)
'CalcCgy Cerr =',r,cerr,maxval(cerr)
6605 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0)
then
6617 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0.and.r.ge.rmax)
then
6621 write(*,*)
'CalcCgy exit rloop',r,cerr,maxval(cerr)
6642 caux =
zadj(k,1)*shat(n0-1,inds(1),inds(2),1) &
6643 +
zadj(k,2)*shat(n0-1,inds(1),inds(2),2) &
6644 -
zadjf(k)*c(n0-1,inds(1),inds(2))
6647 caux = caux -
detz*c(n0-1,inds(1),inds(2))
6652 caux = caux - 2*nlt*
zadj(k,lt)*c(n0,inds(1),inds(2))
6655 c(n0,inds0(1),inds0(2)) = caux/(2*(nl+1)*
zadj(k,l))
6666 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
6667 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
6673 write(*,*)
'CalcCgy final err',cerr
6674 write(*,*)
'CalcCgy final acc',cerr/abs(c(0,0,0))
6679 write(*,*)
'CalcCgy rmax',rmax
6683 write(*,*)
'CalcCgy out',r,n0,n1,r-2*n0-n1,c(n0,n1,r-2*n0-n1)
6705 subroutine calccgyo(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,ordgy_min,ordgy_max,id,Cerr,acc_req_Cr,Cerr2)
6709 integer,
intent(in) :: rmax,ordgy_min,ordgy_max,id
6710 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
6711 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
6712 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
6713 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
6714 double precision,
intent(in) :: acc_req_Cr(0:rmax)
6715 double complex,
allocatable :: Cexpgy(:,:,:,:), CuvExpgy(:,:,:)
6716 double complex,
allocatable :: B_0(:,:,:), B_i(:,:,:), Shat(:,:,:,:)
6717 double complex,
allocatable :: Buv_0(:,:,:), Buv_i(:,:,:)
6718 double complex :: Smod, Caux
6719 double complex :: C0_coli, elimminf2_coli
6720 double precision,
allocatable :: C00_err(:),Cij_err(:)
6721 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
6722 double precision :: B_err,B_max
6723 double precision :: maxCexpgy(0:1,0:rmax+2*ordgy_min,0:ordgy_max),truncfacexp
6724 integer :: rmaxB,rmaxExp,gtrunc,r,n0,n1,n2,a,b,j,sgnab,g,rg
6725 integer :: inds0(2),inds(2),at,bt,k,l,lt,nl,nlt
6726 integer :: bin,nid(0:2)
6729 write(*,*)
'CalcCgy in ',rmax,ordgy_min,ordgy_max,id
6732 write(*,*)
'CalcCgy in ',rmax,ordgy_min,ordgy_max,id
6737 rmaxb = rmax + 2*ordgy_min + 1
6738 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
6739 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
6740 allocate(b_i(0:rmaxb,0:rmaxb,2))
6741 allocate(buv_i(0:rmaxb,0:rmaxb,2))
6747 if (mod(id/bin,2).eq.0)
then
6754 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
6755 call calcb(b_i(:,:,1),buv_i(:,:,1),p20,m02,m22,rmaxb,nid(1))
6756 call calcb(b_i(:,:,2),buv_i(:,:,2),p10,m02,m12,rmaxb,nid(2))
6763 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
6764 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
6767 b_max = max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
6768 b_err = acc_def_b*b_max
6781 if (abs(
detz/( 4d0*
q10*
q20 +
z(2,1)*
z(2,1))).lt.1d-4)
then
6806 allocate(shat(0:rmaxb,0:rmaxb,0:rmaxb,2))
6813 shat(n0,n1,n2,:) = -b_0(n0,n1,n2)
6817 shat(n0,0,k,1) = shat(n0,0,k,1) + b_i(n0,k,1)
6818 shat(n0,k,0,2) = shat(n0,k,0,2) + b_i(n0,k,2)
6824 if (abs(
xadj(1,1)).ge.abs(
xadj(2,2)))
then
6825 if (abs(
xadj(1,1)).ge.abs(
xadj(1,2)))
then
6839 if (abs(
xadj(2,2)).ge.abs(
xadj(1,2)))
then
6854 if (abs(
zadj(1,1)).ge.abs(
zadj(2,2)))
then
6855 if (abs(
zadj(1,1)).ge.abs(
zadj(1,2)))
then
6865 if (abs(
zadj(2,2)).ge.abs(
zadj(1,2)))
then
6878 allocate(cexpgy(0:max(rmax/2,1),0:rmaxexp-2,0:rmaxexp-2,0:ordgy_max))
6881 allocate(cuvexpgy(0:rmaxexp,0:rmaxexp,0:rmaxexp))
6883 cuv(0:rmax,0:rmax,0:rmax) = cuvexpgy(0:rmax,0:rmax,0:rmax)
6886 allocate(c00_err(0:rmaxexp))
6887 allocate(cij_err(0:rmaxexp))
6888 allocate(c00_err2(0:rmaxexp))
6889 allocate(cij_err2(0:rmaxexp))
6912 rloop:
do r=0,rmaxexp-2
6914 if (r.gt.rmax+2*gtrunc+2)
exit rloop
6921 maxcexpgy(1,r,0)=0d0
6929 caux =
zadj(k,1)*shat(0,inds(1),inds(2),1) &
6930 +
zadj(k,2)*shat(0,inds(1),inds(2),2)
6934 caux = caux - 2*nlt*
zadj(k,lt)*cexpgy(1,inds(1),inds(2),0)
6937 cexpgy(1,inds0(1),inds0(2),0) = caux/(2*(nl+1)*
zadj(k,l))
6938 maxcexpgy(1,r,0) = maxcexpgy(1,r,0) + abs(cexpgy(1,inds0(1),inds0(2),0) )
6939 if (r+2.le.rmax)
then
6940 c(1,inds0(1),inds0(2)) = cexpgy(1,inds0(1),inds0(2),0)
6946 maxcexpgy(0,r,0)=0d0
6952 caux = (2*(2+r)*cexpgy(1,n1,n2,0) - 4*cuvexpgy(1,n1,n2) &
6953 - b_0(0,n1,n2))*
z(a,b)
6957 smod = shat(0,n1,n2,a)
6959 if (inds(a).ge.1)
then
6961 smod = smod - 2d0*(inds(a)+1)*cexpgy(1,inds(1),inds(2),0)
6965 caux = caux -
f(b)*smod
6969 cexpgy(0,n1,n2,0) = caux/
xadj(a,b)
6970 maxcexpgy(0,r,0) = maxcexpgy(0,r,0) + abs(cexpgy(0,n1,n2,0))
6972 c(0,n1,n2) = cexpgy(0,n1,n2,0)
6980 cerr(r-1) =
fac_gy*maxcexpgy(0,r,0)
6984 c00_err(r+2) = b_err
6985 cij_err(r)=max(abs(
zadj(a,b))/abs(
xadj(a,b))*max(b_err,c00_err(r+2)), &
6986 fmax/abs(
xadj(a,b))*max(b_err,c00_err(r+1)))
6988 c00_err2(r+2) = b_err
6989 cij_err2(r)=max(abs(
zadj(a,b))/abs(
xadj(a,b))*max(b_err,c00_err(r+2)), &
6990 fmax/abs(
xadj(a,b))*max(b_err,c00_err2(r+1)))
6998 gloop:
do g=1,min(gtrunc,r/2)
7004 maxcexpgy(1,rg,g) = 0d0
7012 caux = -
zadjf(k)*cexpgy(0,inds(1),inds(2),g-1)
7015 caux = caux -
detz*cexpgy(0,inds(1),inds(2),g-1)
7020 caux = caux - 2*nlt*
zadj(k,lt)*cexpgy(1,inds(1),inds(2),g)
7023 cexpgy(1,inds0(1),inds0(2),g) = caux/(2*(nl+1)*
zadj(k,l))
7024 maxcexpgy(1,rg,g) = maxcexpgy(1,rg,g) + abs(cexpgy(1,inds0(1),inds0(2),g) )
7027 if (g.eq.1.and.abs(cexpgy(1,inds0(1),inds0(2),g)).gt. &
7028 truncfacexp*max(1d0,maxcexpgy(1,rg,g-1)) .or. &
7029 g.ge.2.and.abs(cexpgy(1,inds0(1),inds0(2),g)).gt. &
7030 truncfacexp*maxcexpgy(1,rg,g-1))
then
7033 write(*,*)
'CalcCgy cycle loop',n1,n2,g,abs(cexpgy(1,inds0(1),inds0(2),g)),maxcexpgy(1,rg,g-1)
7046 if (rg+2.le.rmax)
then
7051 c(1,inds0(1),inds0(2)) = c(1,inds0(1),inds0(2)) &
7052 + cexpgy(1,inds0(1),inds0(2),g)
7058 maxcexpgy(0,rg,g) = 0d0
7064 caux = 2*(2+rg)*cexpgy(1,n1,n2,g)*
z(a,b)
7068 if (inds0(a).ge.1)
then
7071 caux = caux + 2d0*
f(b)*inds0(a)*cexpgy(1,inds(1),inds(2),g)
7076 inds0(at) = inds0(at)+1
7077 caux = caux - sgnab*
zadjf(bt)*cexpgy(0,inds0(1),inds0(2),g-1)
7081 cexpgy(0,n1,n2,g) = caux/
xadj(a,b)
7085 maxcexpgy(0,rg,g) = maxcexpgy(0,rg,g) + abs(cexpgy(0,n1,n2,g))
7087 if (g.eq.1.and.abs(cexpgy(0,n1,n2,g)).gt. &
7088 truncfacexp*max(1d0/
m2scale,maxcexpgy(0,rg,g-1)).or. &
7089 g.ge.2.and.abs(cexpgy(0,n1,n2,g)).gt. &
7090 truncfacexp*maxcexpgy(0,rg,g-1))
then
7093 write(*,*)
'CalcCgy cycle loop',n1,n2,g,abs(cexpgy(0,n1,n2,g)),maxcexpgy(0,rg,g-1)
7110 c00_err(rg+2) =max(c00_err(rg+2), &
7111 max(abs(
zadjf(k))*cij_err(rg+1),abs(
detz)*cij_err(rg+2))/abs(
zadj(k,l)))
7113 cij_err(rg)= max( cij_err(rg), &
7114 max(abs(
z(a,b))*c00_err(rg+2),abs(
f(b))*c00_err(rg+1), &
7115 abs(
zadjf(b))*cij_err(rg+1))/abs(
xadj(a,b)))
7118 c00_err2(rg+2) =max(c00_err2(rg+2), &
7119 max(abs(
zadjf(k))*cij_err2(rg+1),abs(
detz)*cij_err2(rg+2))/abs(
zadj(k,l)))
7121 cij_err2(rg)= max( cij_err2(rg), &
7122 max(abs(
z(a,b))*c00_err2(rg+2),abs(
f(b))*c00_err2(rg+1), &
7123 abs(
zadjf(b))*cij_err2(rg+1))/abs(
xadj(a,b)))
7126 if (rg+2.le.rmax)
then
7131 c(1,inds0(1),inds0(2)) = c(1,inds0(1),inds0(2)) &
7132 + cexpgy(1,inds0(1),inds0(2),g)
7137 if ((rg.le.rmax))
then
7141 c(0,n1,n2) = c(0,n1,n2) + cexpgy(0,n1,n2,g)
7142 if(abs(cexpgy(0,n1,n2,g-1)).ne.0d0)
then
7143 cerr(rg)=max(cerr(rg),abs(cexpgy(0,n1,n2,g))*min(1d0,abs(cexpgy(0,n1,n2,g))/abs(cexpgy(0,n1,n2,g-1))))
7145 cerr(rg)=max(cerr(rg),abs(cexpgy(0,n1,n2,g)))
7150 if(cij_err(rg).gt.cerr(rg))
then
7151 gtrunc = min(g,gtrunc)
7155 write(*,*)
'CalcCgy exit err',r,g,gtrunc
7168 write(*,*)
'CalcCgy Cerr r =',r
7169 write(*,*)
'CalcCgy Cerr r =',r,cerr
7170 write(*,*)
'CalcCgy Cacc r =',r,cerr/abs(c(0,0,0))
7171 write(*,*)
'CalcCgy Cij_err =',r,cij_err
7174 cerr2 = max(cerr,cij_err2(0:rmax))
7175 cerr = max(cerr,cij_err(0:rmax))
7178 write(*,*)
'CalcCgy Cerr =',r,cerr,maxval(cerr)
7184 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0)
then
7196 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0.and.r.ge.rmax)
then
7214 caux =
zadj(k,1)*shat(n0-1,inds(1),inds(2),1) &
7215 +
zadj(k,2)*shat(n0-1,inds(1),inds(2),2) &
7216 -
zadjf(k)*c(n0-1,inds(1),inds(2))
7219 caux = caux -
detz*c(n0-1,inds(1),inds(2))
7224 caux = caux - 2*nlt*
zadj(k,lt)*c(n0,inds(1),inds(2))
7227 c(n0,inds0(1),inds0(2)) = caux/(2*(nl+1)*
zadj(k,l))
7238 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
7239 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
7245 write(*,*)
'CalcCgy final err',cerr
7246 write(*,*)
'CalcCgy final acc',cerr/abs(c(0,0,0))
7251 write(*,*)
'CalcCgyo rmax',rmax
7255 write(*,*)
'CalcCgyo out',r,n0,n1,r-2*n0-n1,c(n0,n1,r-2*n0-n1)
7273 subroutine calccgp(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,ordgp_min,ordgp_max,id,Cerr,acc_req_Cr,Cerr2)
7277 integer,
intent(in) :: rmax,ordgp_min,ordgp_max,id
7278 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
7279 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
7280 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
7281 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
7282 double precision,
intent(in) :: acc_req_Cr(0:rmax)
7283 double complex,
allocatable :: Cexpgp(:,:,:,:), CuvExpgp(:,:,:)
7284 double complex,
allocatable :: B_0(:,:,:), B_k(:,:), Shat(:,:,:)
7285 double complex,
allocatable :: Buv_0(:,:,:), Buv_k(:,:)
7286 double complex :: Smod, fk, elimminf2_coli
7287 double precision,
allocatable :: C00_err(:),Cij_err(:)
7288 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
7289 double precision :: B_err,B_max
7290 double precision :: maxCexpgp(0:1,0:rmax+ordgp_min+1,0:ordgp_max),truncfacexp
7291 integer :: rmaxB,rmaxExp,gtrunc,r,n0,n1,n2,k,l,g,rg
7292 integer :: bin,nid(0:2),i
7295 write(*,*)
'CalcCgp in ',rmax,ordgp_min,ordgp_max,id
7298 write(*,*)
'CalcCgp in ',rmax,ordgp_min,ordgp_max,id
7320 if (abs(
f(1)).ge.abs(
f(2)))
then
7329 rmaxb = rmax + ordgp_min
7330 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
7331 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
7332 allocate(b_k(0:rmaxb,0:rmaxb))
7333 allocate(buv_k(0:rmaxb,0:rmaxb))
7339 if (mod(id/bin,2).eq.0)
then
7346 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
7348 call calcb(b_k(:,:),buv_k(:,:),p20,m02,m22,rmaxb,nid(1))
7350 call calcb(b_k(:,:),buv_k(:,:),p10,m02,m12,rmaxb,nid(2))
7358 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
7359 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
7360 b_max = max(b_max,abs(b_0(0,n1,n2)))
7363 b_max = max(b_max,maxval(abs(b_k(0,0:rmaxb))))
7364 b_err = acc_def_b*b_max
7367 write(*,*)
'CalcCgp B_max ', b_max
7371 allocate(shat(0:rmaxb,0:rmaxb,0:rmaxb))
7378 shat(n0,n1,n2) = -b_0(n0,n1,n2)
7383 shat(n0,0,l) = shat(n0,0,l) + b_k(n0,l)
7385 shat(n0,l,0) = shat(n0,l,0) + b_k(n0,l)
7394 allocate(cexpgp(0:rmaxexp/2,0:rmaxexp-1,0:rmaxexp-1,0:ordgp_max))
7398 allocate(cuvexpgp(0:rmaxexp,0:rmaxexp,0:rmaxexp))
7400 cuv(0:rmax,0:rmax,0:rmax) = cuvexpgp(0:rmax,0:rmax,0:rmax)
7403 allocate(c00_err(0:rmaxexp))
7404 allocate(cij_err(0:rmaxexp))
7405 allocate(c00_err2(0:rmaxexp))
7406 allocate(cij_err2(0:rmaxexp))
7427 write(*,*)
'CalcCgp rmaxExp',rmaxexp,rmax,gtrunc
7431 rloop:
do r=1,rmaxexp
7434 write(*,*)
'CalcCgp r',r,rmax+gtrunc+1
7438 if (r.gt.rmax+gtrunc+1)
exit rloop
7447 maxcexpgp(1,r,0)=0d0
7452 cexpgp(n0,n1,n2,0) = (2d0*cuvexpgp(n0,n1,n2) + b_0(n0-1,n1,n2) &
7453 +
mm02*cexpgp(n0-1,n1,n2,0))/((r-n0)+1d0)/2d0
7456 maxcexpgp(1,r,0) = maxcexpgp(1,r,0) + abs(cexpgp(n0,n1,n2,0) )
7459 if (r-n0.le.rmax)
then
7460 c(n0,n1,n2) = cexpgp(n0,n1,n2,0)
7469 maxcexpgp(0,r-1,0)=0d0
7473 smod = shat(0,n1,n2)
7474 if ((k.eq.1).and.(n1.ge.1))
then
7475 smod = smod - 2d0*n1*cexpgp(1,n1-1,n2,0)
7476 else if ((k.eq.2).and.(n2.ge.1))
then
7477 smod = smod - 2d0*n2*cexpgp(1,n1,n2-1,0)
7480 cexpgp(0,n1,n2,0) = smod/fk
7481 maxcexpgp(0,r-1,0) = maxcexpgp(0,r-1,0) + abs(cexpgp(0,n1,n2,0))
7483 if (r.le.rmax+1)
then
7484 c(0,n1,n2) = cexpgp(0,n1,n2,0)
7489 if (r.le.rmax+1)
then
7491 cerr(r-1) =
fac_gp*maxcexpgp(0,r-1,0)
7496 c00_err(r) = b_err/(2*r)
7498 cij_err(r-1) = b_err/abs(fk)
7501 c00_err2(r) = b_err/(2*r)
7503 cij_err2(r-1) = b_err/abs(fk)
7506 write(*,*)
'CalcCgp B_err',b_err,abs(fk), cij_err(r-1),r
7514 gloop:
do g=1,min(gtrunc,r-1)
7520 maxcexpgp(1,rg,g) = 0d0
7525 cexpgp(n0,n1,n2,g) = (2d0*
mm02*cexpgp(n0-1,n1,n2,g) &
7526 -
z(1,1)*cexpgp(n0-1,n1+2,n2,g-1) - 2d0*
z(2,1)*cexpgp(n0-1,n1+1,n2+1,g-1) &
7527 -
z(2,2)*cexpgp(n0-1,n1,n2+2,g-1))/((rg-n0)+1d0)/4d0
7530 maxcexpgp(1,rg,g) = maxcexpgp(1,rg,g) + abs(cexpgp(n0,n1,n2,g))
7533 if (g.eq.1.and.abs(cexpgp(1,n1,n2,g)).gt. &
7534 truncfacexp*max(1d0,maxcexpgp(1,rg,g-1)) .or. &
7535 g.ge.2.and.abs(cexpgp(1,n1,n2,g)).gt. &
7536 truncfacexp*maxcexpgp(1,rg,g-1))
then
7539 write(*,*)
'CalcCg exit rloop',n0,n1,n2,g,abs(cexpgp(n0,n1,n2,g)),maxcexpgp(1,rg,g-1)
7556 if (rg-n0.le.rmax)
then
7559 c(n0,n1,n2) = c(n0,n1,n2) + cexpgp(n0,n1,n2,g)
7568 maxcexpgp(0,rg-1,g) = 0d0
7572 smod = -
z(1,k)*cexpgp(0,n1+1,n2,g-1) &
7573 -
z(2,k)*cexpgp(0,n1,n2+1,g-1)
7574 if ((k.eq.1).and.(n1.ge.1))
then
7575 smod = smod - 2d0*n1*cexpgp(1,n1-1,n2,g)
7576 else if ((k.eq.2).and.(n2.ge.1))
then
7577 smod = smod - 2d0*n2*cexpgp(1,n1,n2-1,g)
7580 cexpgp(0,n1,n2,g) = smod/fk
7582 maxcexpgp(0,rg-1,g) = maxcexpgp(0,rg-1,g) + abs(cexpgp(0,n1,n2,g))
7584 if (g.eq.1.and.abs(cexpgp(0,n1,n2,g)).gt. &
7585 truncfacexp*max(1/
m2max,maxcexpgp(0,rg-1,g-1)) .or. &
7586 g.ge.2.and.abs(cexpgp(0,n1,n2,g)).gt. &
7587 truncfacexp*maxcexpgp(0,rg-1,g-1))
then
7590 write(*,*)
'CalcCgp exit gloop',0,n1,n2,g,abs(cexpgp(0,n1,n2,g)),maxcexpgp(0,rg,g-1)
7600 c00_err(rg) = max(c00_err(rg),max(2*abs(m02)*cij_err(rg-2),
maxz*cij_err(rg))/(4*r) )
7602 cij_err(rg-1) = max(cij_err(rg-1),max(2*c00_err(rg),
maxz*cij_err(rg))/abs(fk) )
7605 c00_err2(rg) = max(c00_err2(rg),max(2*abs(m02)*cij_err2(rg-2),
maxz*cij_err2(rg))/(4*r) )
7607 cij_err2(rg-1) = max(cij_err2(rg-1),max(2*c00_err2(rg),
maxz*cij_err2(rg))/abs(fk) )
7611 if (rg-n0.le.rmax)
then
7614 c(n0,n1,n2) = c(n0,n1,n2) + cexpgp(n0,n1,n2,g)
7620 if ((rg.le.rmax+1))
then
7624 c(0,n1,n2) = c(0,n1,n2) + cexpgp(0,n1,n2,g)
7626 if(abs(cexpgp(0,n1,n2,g-1)).ne.0d0)
then
7627 cerr(rg-1)=max(cerr(rg-1),abs(cexpgp(0,n1,n2,g))*min(1d0,abs(cexpgp(0,n1,n2,g))/abs(cexpgp(0,n1,n2,g-1))))
7629 cerr(rg-1)=max(cerr(rg-1),abs(cexpgp(0,n1,n2,g)))
7634 if(cij_err(rg-1).gt.cerr(rg-1))
then
7635 gtrunc = min(g,gtrunc)
7638 write(*,*)
'CalcCgp exit err',r,g,gtrunc
7648 write(*,*)
'CalcCgp Cerr r =',r,cerr
7649 write(*,*)
'CalcCgp Cacc r =',r,cerr/abs(c(0,0,0))
7650 write(*,*)
'CalcCgp Cij_err =',r,cij_err
7653 cerr2 = max(cerr,cij_err2(0:rmax))
7654 cerr = max(cerr,cij_err(0:rmax))
7657 write(*,*)
'CalcCgp Cerr =',r,cerr,maxval(cerr)
7658 write(*,*)
'CalcCgp accr =',r,acc_req_cr*abs(c(0,0,0)),maxval(acc_req_cr*abs(c(0,0,0)))
7659 write(*,*)
'CalcCgp C-ar =',r,cerr-acc_req_cr*abs(c(0,0,0)),maxval(cerr-acc_req_cr*abs(c(0,0,0)))
7665 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0)
then
7669 c(n0,n1,rg-2*n0-n1)=0d0
7679 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0.and.r.gt.rmax)
then
7695 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
7696 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
7702 write(*,*)
'CalcCgp final err',cerr
7703 write(*,*)
'CalcCgp final acc',cerr/abs(c(0,0,0))
7707 write(*,*)
'CalcCgp rmax',rmax
7711 write(*,*)
'CalcCgp out',r,n0,n1,r-2*n0-n1,c(n0,n1,r-2*n0-n1)
7740 subroutine calccgpf(C,Cuv,p10,p21,p20,m02,m12,m22,rmax,ordgpf_min,ordgpf_max,id,Cerr,acc_req_Cr,Cerr2)
7744 integer,
intent(in) :: rmax,ordgpf_min,ordgpf_max,id
7745 double complex,
intent(in) :: p10,p21,p20,m02,m12,m22
7746 double complex,
intent(out) :: C(0:rmax,0:rmax,0:rmax)
7747 double complex,
intent(out) :: Cuv(0:rmax,0:rmax,0:rmax)
7748 double precision,
intent(out) :: Cerr(0:rmax),Cerr2(0:rmax)
7749 double precision,
intent(in) :: acc_req_Cr(0:rmax)
7750 double complex,
allocatable :: Cexpgpf(:,:,:,:), CuvExpgpf(:,:,:)
7751 double complex,
allocatable :: B_0(:,:,:), B_i(:,:,:), Shat(:,:,:,:)
7752 double complex,
allocatable :: Buv_0(:,:,:), Buv_i(:,:,:)
7753 double complex :: Smod, Caux, Zadj2f
7754 double complex :: C0_coli, elimminf2_coli
7755 double precision,
allocatable :: C00_err(:),Cij_err(:)
7756 double precision,
allocatable :: C00_err2(:),Cij_err2(:)
7757 double precision :: B_err,B_max,aZadj2f
7758 double precision :: maxCexpgpf(0:1,0:rmax+2*ordgpf_min,0:ordgpf_max),truncfacexp
7759 double precision :: minZk
7760 integer :: rmaxB,rmaxExp,gtrunc,r,n0,n1,n2,i,j,jt,g,rg
7761 integer :: inds0(2),inds(2),k,l,lt,nl,nlt
7762 integer :: bin,nid(0:2)
7765 write(*,*)
'CalcCgpf in ',rmax,ordgpf_min,ordgpf_max,id
7766 write(*,*)
'CalcCgpf in ',p10,p21,p20,m02,m12,m22
7769 write(*,*)
'CalcCgpf in ',rmax,ordgpf_min,ordgpf_max,id
7774 rmaxb = rmax + 2*ordgpf_min + 1
7775 allocate(b_0(0:rmaxb,0:rmaxb,0:rmaxb))
7776 allocate(buv_0(0:rmaxb,0:rmaxb,0:rmaxb))
7777 allocate(b_i(0:rmaxb,0:rmaxb,2))
7778 allocate(buv_i(0:rmaxb,0:rmaxb,2))
7784 if (mod(id/bin,2).eq.0)
then
7791 call calcb(b_0(:,0,:),buv_0(:,0,:),p21,m12,m22,rmaxb,nid(0))
7792 call calcb(b_i(:,:,1),buv_i(:,:,1),p20,m02,m22,rmaxb,nid(1))
7793 call calcb(b_i(:,:,2),buv_i(:,:,2),p10,m02,m12,rmaxb,nid(2))
7800 b_0(0:n0,n1,n2) = -b_0(0:n0,n1-1,n2)-b_0(0:n0,n1-1,n2+1)
7801 buv_0(0:n0,n1,n2) = -buv_0(0:n0,n1-1,n2)-buv_0(0:n0,n1-1,n2+1)
7804 b_max = max(b_max,maxval(abs(b_i(0,0:rmaxb,1:2))))
7805 b_err = acc_def_b*b_max
7818 if (abs(
detz).lt.abs(4d0*
q10*
q20 +
z(2,1)*
z(2,1))*1d-4)
then
7843 allocate(shat(0:rmaxb,0:rmaxb,0:rmaxb,2))
7850 shat(n0,n1,n2,:) = -b_0(n0,n1,n2)
7854 shat(n0,0,k,1) = shat(n0,0,k,1) + b_i(n0,k,1)
7855 shat(n0,k,0,2) = shat(n0,k,0,2) + b_i(n0,k,2)
7862 if (maxval(abs(
z(1,1:2))).le.minzk)
then
7863 minzk = maxval(abs(
z(1,1:2)))
7868 if (maxval(abs(
z(2,1:2))).lt.minzk)
then
7869 minzk = maxval(abs(
z(2,1:2)))
7876 write(*,*)
'CalcCgpf: minZk',k,minzk
7883 allocate(cexpgpf(0:max(rmax/2,1),0:rmaxexp-2,0:rmaxexp-2,0:ordgpf_max))
7886 allocate(cuvexpgpf(0:rmaxexp,0:rmaxexp,0:rmaxexp))
7888 cuv(0:rmax,0:rmax,0:rmax) = cuvexpgpf(0:rmax,0:rmax,0:rmax)
7891 allocate(c00_err(0:rmaxexp))
7892 allocate(cij_err(0:rmaxexp))
7893 allocate(c00_err2(0:rmaxexp))
7894 allocate(cij_err2(0:rmaxexp))
7915 write(*,*)
'CalcCgpf: gtrunc orig=',gtrunc
7916 write(*,*)
'CalcCgpf: rmaxExp-2=',rmaxexp-2
7921 rloop:
do r=0,rmaxexp-2
7924 write(*,*)
'CalcCgpf: rloop=',r,rmaxexp-2,rmax+2*gtrunc+2
7925 write(*,*)
'CalcCgpf: rloop=',rmax,gtrunc
7928 if (r.gt.rmax+2*gtrunc+2)
exit rloop
7935 maxcexpgpf(1,r,0)=0d0
7944 caux = shat(0,inds(1),inds(2),k)
7946 cexpgpf(1,inds0(1),inds0(2),0) = caux/(2*(nl+1))
7947 maxcexpgpf(1,r,0) = maxcexpgpf(1,r,0) + abs(cexpgpf(1,inds0(1),inds0(2),0) )
7949 if (r+1.le.rmax)
then
7950 c(1,inds0(1),inds0(2)) = cexpgpf(1,inds0(1),inds0(2),0)
7956 maxcexpgpf(0,r,0)=0d0
7962 caux = 2*(4+r+r)*cexpgpf(1,n1,n2,0) - 4*cuvexpgpf(1,n1,n2) &
7965 cexpgpf(0,n1,n2,0) = caux/(2d0*m02)
7967 maxcexpgpf(0,r,0) = maxcexpgpf(0,r,0) + abs(cexpgpf(0,n1,n2,0))
7969 c(0,n1,n2) = cexpgpf(0,n1,n2,0)
7976 cerr(r) =
fac_gpf*maxcexpgpf(0,r,0)
7983 c00_err(r+2) = b_err /2d0
7984 cij_err(r) = max(b_err,2*(r+2)*c00_err(r+2))/abs(m02)
7986 c00_err2(r+2) = b_err /2d0
7987 cij_err2(r) = max(b_err,2*(r+2)*c00_err2(r+2))/abs(m02)
7991 write(*,*)
'CalcCgpf leading terms r =',r
7992 write(*,*)
'CalcCgpf Cij_err =',r,cij_err(0:r)
7993 write(*,*)
'CalcCgpf Cexp0(1,0,0)=',r,c(1,0,0),cexpgpf(1,0,0,0)
7994 write(*,*)
'CalcCgpf Cexp0(0,0,0)=',r,c(0,0,0),cexpgpf(0,0,0,0)
8005 gloop:
do g=1,min(gtrunc,r/2)
8009 maxcexpgpf(1,rg,g) = 0d0
8017 caux = -
f(k)*cexpgpf(0,inds(1),inds(2),g-1)
8020 caux = caux -
z(k,l)*cexpgpf(0,inds(1),inds(2),g-1)
8023 inds(lt) = inds(lt)+1
8024 caux = caux -
z(k,lt)*cexpgpf(0,inds(1),inds(2),g-1)
8026 cexpgpf(1,inds0(1),inds0(2),g) = caux/(2*(nl+1))
8028 maxcexpgpf(1,rg,g) = maxcexpgpf(1,rg,g) + abs(cexpgpf(1,inds0(1),inds0(2),g) )
8030 if (g.eq.1.and.abs(cexpgpf(1,inds0(1),inds0(2),g)).gt. &
8031 truncfacexp*max(1d0,maxcexpgpf(1,rg,g-1)) .or. &
8032 g.ge.2.and.abs(cexpgpf(1,inds0(1),inds0(2),g)).gt. &
8033 truncfacexp*maxcexpgpf(1,rg,g-1))
then
8035 write(*,*)
'CalcCgpf exit gloop',n1,n2,g,abs(cexpgpf(1,inds0(1),inds0(2),g)),maxcexpgpf(1,rg,g-1)
8036 write(*,*)
'CalcCgpf exit gloop',g,rg,inds0(1),inds0(2)
8048 if (rg+1.le.rmax)
then
8053 c(1,inds0(1),inds0(2)) = c(1,inds0(1),inds0(2)) &
8054 + cexpgpf(1,inds0(1),inds0(2),g)
8060 maxcexpgpf(0,rg,g) = 0d0
8066 caux = 2*(4+rg+rg)*cexpgpf(1,n1,n2,g)
8072 caux = caux +
z(i,j)*cexpgpf(0,inds(1),inds(2),g-1)
8078 cexpgpf(0,n1,n2,g) = caux/(2*m02)
8080 maxcexpgpf(0,rg,g) = maxcexpgpf(0,rg,g) + abs(cexpgpf(0,n1,n2,g))
8082 if (g.eq.1.and.abs(cexpgpf(0,n1,n2,g)).gt. &
8083 truncfacexp*max(1d0/
m2scale,maxcexpgpf(0,rg,g-1)).or. &
8084 g.ge.2.and.abs(cexpgpf(0,n1,n2,g)).gt. &
8085 truncfacexp*maxcexpgpf(0,rg,g-1))
then
8088 write(*,*)
'CalcCgpf exit gloop',n1,n2,g,rg
8089 write(*,*)
'CalcCgpf exit gloop',abs(cexpgpf(0,n1,n2,g)),maxcexpgpf(0,rg,g-1),1d0/
m2scale
8090 write(*,*)
'CalcCgpf exit gloop',truncfacexp
8105 write(*,*)
'CalcCgpf expansion terms r =',r,g,rg
8106 write(*,*)
'CalcCgpf Cij_err =',r,cij_err(0:r)
8107 write(*,*)
'CalcCgpf Cexp(1,0,0)=',r,g,cexpgpf(1,0,0,g)
8108 write(*,*)
'CalcCgpf Cexp(0,0,0)=',r,g,cexpgpf(0,0,0,g)
8113 c00_err(rg+2) =max(c00_err(rg+2), &
8114 fmax/2d0*cij_err(rg+1), &
8115 maxz/2d0*cij_err(rg+2))
8119 write(*,*)
'CalcCgpf test2',rg,i,j,cij_err(rg)
8120 write(*,*)
'CalcCgpf test2',rg,cij_err(rg+1),c00_err(rg+1)
8123 cij_err(rg)= max( cij_err(rg), &
8124 2*(rg+2)/abs(m02)*c00_err(rg+2), &
8125 maxz/(2*abs(m02))*cij_err(rg+2))
8128 c00_err2(rg+2) =max(c00_err2(rg+2), &
8129 fmax/2d0*cij_err(rg+1), &
8130 maxz/2d0*cij_err(rg+2))
8133 cij_err2(rg)= max( cij_err2(rg), &
8134 2*(rg+2)/abs(m02)*c00_err2(rg+2), &
8135 maxz/(2*abs(m02))*cij_err2(rg+2))
8139 if (rg+1.le.rmax)
then
8144 c(1,inds0(1),inds0(2)) = c(1,inds0(1),inds0(2)) &
8145 + cexpgpf(1,inds0(1),inds0(2),g)
8150 if ((rg.le.rmax))
then
8154 c(0,n1,n2) = c(0,n1,n2) + cexpgpf(0,n1,n2,g)
8157 write(*,*)
'CalcCgpf test1',rg,n1,n2,cerr(rg)
8158 write(*,*)
'CalcCgpf test1',cexpgpf(0,n1,n2,g)
8159 write(*,*)
'CalcCgpf test1',cexpgpf(0,n1,n2,g-1)
8162 if(abs(cexpgpf(0,n1,n2,g-1)).ne.0d0)
then
8163 cerr(rg)=max(cerr(rg),abs(cexpgpf(0,n1,n2,g))*min(1d0,abs(cexpgpf(0,n1,n2,g))/abs(cexpgpf(0,n1,n2,g-1))))
8165 cerr(rg)=max(cerr(rg),abs(cexpgpf(0,n1,n2,g)))
8169 write(*,*)
'CalcCgpf test1',cerr(rg)
8173 write(*,*)
'CalcCgpf expansion terms r =',r,g,rg
8174 write(*,*)
'CalcCgpf Cij_err =',r,cij_err(0:min(r,rmax))
8175 write(*,*)
'CalcCgpf Cerr =',r,cerr(0:min(r,rmax))
8176 write(*,*)
'CalcCgpf C(1,0,0)=',r,c(1,0,0)
8177 write(*,*)
'CalcCgpf C(0,0,0)=',r,c(0,0,0)
8183 if(cij_err(rg).gt.cerr(rg))
then
8184 gtrunc = min(g,gtrunc)
8188 write(*,*)
'CalcCgpf adjust gtrunc',r,g,gtrunc
8200 write(*,*)
'CalcCgpf Cerr r =',r
8201 write(*,*)
'CalcCgpf Cerr r =',r,cerr(0:min(r,rmax))
8202 write(*,*)
'CalcCgpf Cacc r =',r,cerr/abs(c(0,0,0))
8203 write(*,*)
'CalcCgpf Cij_err =',r,cij_err
8204 write(*,*)
'CalcCgpf C(1,0,0)=',r,c(1,0,0)
8205 write(*,*)
'CalcCgpf C(0,0,0)=',r,c(0,0,0)
8206 if(rmax.ge.1.and.r.ge.1)
then
8207 write(*,*)
'CalcCgpf C(0,1,0)=',r,c(0,1,0)
8208 if(rmax.ge.2.and.r.ge.2)
then
8209 write(*,*)
'CalcCgpf C(0,1,1)=',r,c(0,1,1)
8210 if(rmax.ge.3.and.r.ge.3)
then
8211 write(*,*)
'CalcCgpf C(0,1,2)=',r,c(0,1,2)
8213 write(*,*)
'CalcCgpf C(0,0,4)=',r,c(0,0,4)
8220 cerr2 = max(cerr,cij_err2(0:rmax))
8221 cerr = max(cerr,cij_err(0:rmax))
8224 write(*,*)
'CalcCgpf Cerr =',r,cerr,maxval(cerr)
8230 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0)
then
8242 if(maxval(cerr-acc_req_cr*abs(c(0,0,0))).le.0d0.and.r.ge.rmax)
then
8246 write(*,*)
'CalcCgpf exit rloop',r,cerr,maxval(cerr)
8267 caux = shat(n0-1,inds(1),inds(2),k) &
8268 -
f(k)*c(n0-1,inds(1),inds(2)) &
8269 -
z(k,1)*c(n0-1,inds(1)+1,inds(2)) &
8270 -
z(k,2)*c(n0-1,inds(1),inds(2)+1)
8272 c(n0,inds0(1),inds0(2)) = caux/(2*(nl+1))
8283 c(n0,n1,n2) = (b_0(n0-1,n1,n2) + 2*
mm02*c(n0-1,n1,n2) + 4*cuv(n0,n1,n2) &
8284 +
f(1)*c(n0-1,n1+1,n2) +
f(2)*c(n0-1,n1,n2+1)) / (2*r)
8290 write(*,*)
'CalcCgpf final err',cerr
8291 write(*,*)
'CalcCgpf final acc',cerr/abs(c(0,0,0))
8296 write(*,*)
'CalcCgpf rmax',rmax
8300 write(*,*)
'CalcCgpf out',r,n0,n1,r-2*n0-n1,c(n0,n1,r-2*n0-n1)
8323 subroutine copycimp3(C,C_alt,Cerr,Cerr_alt,Cerr1,Cerr1_alt,Cerr2,Cerr2_alt,Crmethod,Crmethod_alt,rmax,r_alt)
8325 integer,
intent(in) :: rmax,r_alt
8326 double complex,
intent(inout) :: C(0:rmax,0:rmax,0:rmax)
8327 double precision,
intent(inout) :: Cerr(0:rmax),Cerr1(0:rmax),Cerr2(0:rmax)
8328 integer,
intent(inout) :: Crmethod(0:rmax)
8329 double complex,
intent(in) :: C_alt(0:r_alt,0:r_alt,0:r_alt)
8330 double precision,
intent(in) :: Cerr_alt(0:r_alt),Cerr2_alt(0:r_alt),Cerr1_alt(0:r_alt)
8331 integer,
intent(in) :: Crmethod_alt(0:rmax)
8339 if (cerr_alt(r).lt.cerr(r))
then
8340 crmethod(r)=crmethod_alt(r)
8342 cerr1(r)=cerr1_alt(r)
8343 cerr2(r)=cerr2_alt(r)
8346 c(n0,n1,r-n0-n1) = c_alt(n0,n1,r-n0-n1)