323 integer,
intent(in) :: rmax,id
324 double complex,
intent(in) :: p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
327 double complex,
intent(out) :: D(0:rmax,0:rmax,0:rmax,0:rmax)
328 double complex,
intent(out) :: Duv(0:rmax,0:rmax,0:rmax,0:rmax)
329 double complex :: elimminf2_coli
330 double precision,
intent(out) :: Derr1(0:rmax),Derr2(0:rmax)
331 double precision :: D0est,Dtyp
333 double complex :: D0_coli
336 double complex :: chdet
337 integer :: rmaxC,r,rid,n0,n1,n2,n3,g,gy,gp,gr,gm,gpf,i,iexp
338 integer :: bin,k,nid(0:3)
339 logical :: use_pv,use_pv2,use_g,use_gy,use_gp,use_gr,use_gm,use_gpf
341 integer :: r_alt,Drmethod(0:rmax),DrCalc(0:rmax),DCalc
342 double complex,
allocatable :: C_i(:,:,:,:), Cuv_i(:,:,:,:)
343 double complex :: D_alt(0:rmax,0:rmax,0:rmax,0:rmax)
344 double complex :: Duv_alt(0:rmax,0:rmax,0:rmax,0:rmax)
345 double precision :: Derr(0:rmax),Derr_alt(0:rmax),Derr1_alt(0:rmax),Derr2_alt(0:rmax)
346 integer :: Drmethod_alt(0:rmax)
348 double precision :: err_pv(0:rmax),err_pv2(0:rmax),err_g(0:rmax),err_gy(0:rmax),err_gp(0:rmax)
349 double precision :: err_gr(0:rmax),err_gm(0:rmax),err_gpf(0:rmax)
350 double precision :: h_pv,w_pv,v_pv,z_pv,h_pv2,w_pv2,v_pv2,z_pv2,hw_pv2
351 double precision :: u_g,z_g,err_g_C,err_g_Cr,err_g_exp
352 double precision :: u_gm,z_gm,err_gm_C,err_gm_Cr,err_gm_exp
353 double precision :: v1_gy,b_gy,err_gy_C,err_gy_Cr,err_gy_exp
354 double precision :: v_gp,v1_gp,z_gp,err_gp_C,err_gp_Cr,err_gp_exp
355 double precision :: v1_gpf,b_gpf,err_gpf_C,err_gpf_Cr,err_gpf_exp
356 double precision :: x_gr,y_gr,y1_gr,a_gr,err_gr_C,err_gr_Cr,err_gr_exp
357 double precision :: err_C0,Cerr_i(0:rmax_C,0:3),err_C(0:rmax_C),err_D0,acc_D,errfac(0:3),err_req_D,err_inf,Cerr2_i(0:rmax_C,0:3)
358 double precision :: checkest,norm,Dscale2
359 logical :: lerr_D0,errorwriteflag
361 character(len=*),
parameter :: fmt1 =
"(A7,'dcmplx(',d25.18,' , ',d25.18,' )')"
362 character(len=*),
parameter :: fmt10 =
"(A17,'(',d25.18,' , ',d25.18,' )')"
364 #ifdef CritPointsCOLI
365 integer,
parameter :: MaxCritPointD=50
367 integer,
parameter :: MaxCritPointD=0
369 integer,
save :: CritPointCntD
372 data critpointcntd /0/
375 write(*,*)
'CalcDred in',rmax,id,acc_req_d
376 write(*,*)
'CalcDred acc_req',acc_req_d,reqacc_coli
377 write(*,*)
'CalcDred in',p10,p21,p32,p30,p20,p31,m02,m12,m22,m32
388 rmaxc = max(rmax-1,3)
389 allocate(c_i(0:rmaxc,0:rmaxc,0:rmaxc,0:3))
390 allocate(cuv_i(0:rmaxc,0:rmaxc,0:rmaxc,0:3))
396 if (mod(id/bin,2).eq.0)
then
404 call calcc(c_i(:,:,:,0),cuv_i(:,:,:,0),p21,p32,p31,m12,m22,m32,rmaxc,nid(0),cerr_i(0:rmaxc,0),cerr2_i(0:rmaxc,0))
405 call calcc(c_i(:,:,:,1),cuv_i(:,:,:,1),p20,p32,p30,m02,m22,m32,rmaxc,nid(1),cerr_i(0:rmaxc,1),cerr2_i(0:rmaxc,1))
406 call calcc(c_i(:,:,:,2),cuv_i(:,:,:,2),p10,p31,p30,m02,m12,m32,rmaxc,nid(2),cerr_i(0:rmaxc,2),cerr2_i(0:rmaxc,2))
407 call calcc(c_i(:,:,:,3),cuv_i(:,:,:,3),p10,p21,p20,m02,m12,m22,rmaxc,nid(3),cerr_i(0:rmaxc,3),cerr2_i(0:rmaxc,3))
410 write(*,*)
'CalcDred Cerr 0 =',cerr_i(0:rmaxc,0)
411 write(*,*)
'CalcDred Cerr 1 =',cerr_i(0:rmaxc,1)
412 write(*,*)
'CalcDred Cerr 2 =',cerr_i(0:rmaxc,2)
413 write(*,*)
'CalcDred Cerr 3 =',cerr_i(0:rmaxc,3)
414 if (abs(c_i(0,0,0,0)).ne.0d0) &
415 write(*,*)
'CalcDred Cacc 0 =',cerr_i(0:rmaxc,0)/abs(c_i(0,0,0,0))
416 if (abs(c_i(0,0,0,1)).ne.0d0) &
417 write(*,*)
'CalcDred Cacc 1 =',cerr_i(0:rmaxc,1)/abs(c_i(0,0,0,1))
418 if (abs(c_i(0,0,0,2)).ne.0d0) &
419 write(*,*)
'CalcDred Cacc 2 =',cerr_i(0:rmaxc,2)/abs(c_i(0,0,0,2))
420 if (abs(c_i(0,0,0,3)).ne.0d0) &
421 write(*,*)
'CalcDred Cacc 3 =',cerr_i(0:rmaxc,3)/abs(c_i(0,0,0,3))
435 cerr_i(r,i)=cerr_i(r-1,i)*errfac(i)
440 err_c(r)=maxval(cerr_i(r,0:3))
446 write(*,*)
'CalcDred err_C0:',err_c0
447 write(*,*)
'CalcDred Cerr 0 =',cerr_i(0:rmax_c,0)
448 write(*,*)
'CalcDred Cerr 1 =',cerr_i(0:rmax_c,1)
449 write(*,*)
'CalcDred Cerr 2 =',cerr_i(0:rmax_c,2)
450 write(*,*)
'CalcDred Cerr 3 =',cerr_i(0:rmax_c,3)
451 write(*,*)
'CalcDred Cerr =', err_c(0:rmax_c)
464 mm02 = elimminf2_coli(m02)
465 mm12 = elimminf2_coli(m12)
466 mm22 = elimminf2_coli(m22)
467 mm32 = elimminf2_coli(m32)
468 q10 = elimminf2_coli(p10)
469 q21 = elimminf2_coli(p21)
470 q32 = elimminf2_coli(p32)
471 q30 = elimminf2_coli(p30)
472 q31 = elimminf2_coli(p31)
473 q20 = elimminf2_coli(p20)
496 maxz = maxval(abs(
z))
499 if (
detz.ne.0d0)
then
513 zadj(1,1) = (
z(3,3)*
z(2,2)-
z(2,3)*
z(2,3))
514 zadj(2,1) = (
z(1,3)*
z(2,3)-
z(3,3)*
z(1,2))
515 zadj(3,1) = (
z(1,2)*
z(2,3)-
z(2,2)*
z(1,3))
517 zadj(2,2) = (
z(1,1)*
z(3,3)-
z(1,3)*
z(1,3))
518 zadj(3,2) = (
z(1,2)*
z(1,3)-
z(1,1)*
z(2,3))
521 zadj(3,3) = (
z(1,1)*
z(2,2)-
z(1,2)*
z(1,2))
558 fmax = max(abs(f(1)),abs(f(2)),abs(f(3)))
561 mx(1,0) =
q10 - mm12 + mm02
567 mx(1:3,1:3) =
z(1:3,1:3)
571 if (
detx.ne.0d0.and.
maxz.ne.0d0)
then
583 zadj2ff(1,1) = -f(2)*f(2)*
z(3,3) - f(3)*f(3)*
z(2,2) &
585 zadj2ff(2,1) = f(2)*f(1)*
z(3,3) - f(3)*f(1)*
z(2,3) &
586 - f(2)*f(3)*
z(1,3) + f(3)*f(3)*
z(1,2)
587 zadj2ff(3,1) = -f(2)*f(1)*
z(3,2) + f(3)*f(1)*
z(2,2) &
588 + f(2)*f(2)*
z(3,1) - f(3)*f(2)*
z(2,1)
590 zadj2ff(2,2) = -f(1)*f(1)*
z(3,3) - f(3)*f(3)*
z(1,1) &
592 zadj2ff(3,2) = f(1)*f(1)*
z(2,3) - f(1)*f(2)*
z(1,3) &
593 - f(3)*f(1)*
z(2,1) + f(3)*f(2)*
z(1,1)
596 zadj2ff(3,3) = -f(1)*f(1)*
z(2,2) - f(2)*f(2)*
z(1,1) &
641 zadj2f(1,2,1) =
z(3,2)*f(3) -
z(3,3)*f(2)
642 zadj2f(1,3,1) = -
z(2,2)*f(3) +
z(2,3)*f(2)
643 zadj2f(2,3,1) =
z(1,2)*f(3) -
z(1,3)*f(2)
644 zadj2f(1,2,2) = -
z(3,1)*f(3) +
z(3,3)*f(1)
645 zadj2f(1,3,2) =
z(2,1)*f(3) -
z(2,3)*f(1)
646 zadj2f(2,3,2) = -
z(1,1)*f(3) +
z(1,3)*f(1)
647 zadj2f(1,2,3) =
z(3,1)*f(2) -
z(3,2)*f(1)
648 zadj2f(1,3,3) = -
z(2,1)*f(2) +
z(2,2)*f(1)
649 zadj2f(2,3,3) =
z(1,1)*f(2) -
z(1,2)*f(1)
650 zadj2f(2,1,1) = -zadj2f(1,2,1)
651 zadj2f(3,1,1) = -zadj2f(1,3,1)
652 zadj2f(3,2,1) = -zadj2f(2,3,1)
653 zadj2f(2,1,2) = -zadj2f(1,2,2)
654 zadj2f(3,1,2) = -zadj2f(1,3,2)
655 zadj2f(3,2,2) = -zadj2f(2,3,2)
656 zadj2f(2,1,3) = -zadj2f(1,2,3)
657 zadj2f(3,1,3) = -zadj2f(1,3,3)
658 zadj2f(3,2,3) = -zadj2f(2,3,3)
660 maxzadj2f=maxval(abs(zadj2f))
708 dscale2 = max(abs(p10*p32),abs(p21*p30),abs(p20*p31),abs(m02*m02), \
709 abs(m12*m12),abs(m22*m22),abs(m32*m32))
711 d0est = max(abs(d0_coli(p10,p21,p32,p30,p20,p31,m02,m12,m22,m32)), \
716 if(dscale2.ne.0d0)
then
732 err_inf = acc_inf*d0est
736 if (
adetx.ne.0d0)
write(*,*)
'D0est',1d0/sqrt(
adetx)
737 if (
m2max.ne.0d0)
write(*,*)
'D0est',1d0/
m2max**2
738 if (
maxz.ne.0d0)
write(*,*)
'D0est',1d0/
maxz**2
739 write(*,*)
'D0est',d0est
749 dcount(0) = dcount(0)+1
752 if (
adetx.ne.0d0)
then
753 err_d0 = acc_def_d0*max( d0est, 1d0/sqrt(
adetx) )
755 err_d0 = acc_def_d0*d0est
760 err_req_d = acc_req_d * d0est
792 if (mod(rmax,2).eq.1)
then
793 err_pv(rmax) = max( w_pv**((rmax-1)/2) * v_pv * err_d0, &
794 w_pv**((rmax-1)/2) * z_pv * err_c0, z_pv * err_c(rmax-1) )
796 write(*,*)
'CalcDred w_pv: ',w_pv,v_pv,z_pv,err_d0,err_c0,rmax
798 write(*,*)
'CalcDred err_pv con: ',err_pv(rmax), w_pv**((rmax-1)/2) * v_pv * err_d0, &
799 w_pv**((rmax-1)/2) * z_pv * err_c0, z_pv * err_c(rmax-1)
802 err_pv(rmax) = max( w_pv**(rmax/2) * err_d0, &
803 w_pv**(rmax/2-1) * v_pv * z_pv * err_c0, z_pv * err_c(rmax-1) )
805 write(*,*)
'CalcDred w_pv: ',w_pv,v_pv,z_pv
806 write(*,*)
'CalcDred err_pv con: ',err_pv(rmax), w_pv**((rmax)/2) * err_d0, &
807 w_pv**(rmax/2-1) * v_pv * z_pv * err_c0, z_pv * err_c(rmax-1)
808 write(*,*)
'CalcDred err_pv con: ',err_pv(rmax), w_pv**((rmax)/2),err_d0, &
809 w_pv**(rmax/2-1) * v_pv,err_c0, z_pv,err_c(rmax-1)
847 write(*,*)
'CalcDred w_pv2: ',w_pv2,v_pv2,z_pv2
848 write(*,*)
'CalcDred w_pv2: ',mod(rmax,2).eq.1
851 if (mod(rmax,2).eq.1)
then
859 err_pv2(rmax) = max( err_d0 * max(hw_pv2**rmax, &
860 hw_pv2*v_pv2**((rmax-1)/2) ), &
861 err_c0 * z_pv2*max(w_pv2*hw_pv2**(rmax), &
862 max(1d0,w_pv2)*hw_pv2*v_pv2**((rmax-1)/2), &
863 v_pv2**((rmax+1)/2)), &
864 err_c(rmax-1) * z_pv2 * max(hw_pv2,hw_pv2*w_pv2,v_pv2) )
867 write(*,*)
'CalcDred err_pv2: ', &
868 err_pv2(rmax) , err_d0,err_d0*w_pv2**rmax,err_d0*v_pv2**((rmax-1)/2), &
869 err_d0*w_pv2*v_pv2**((rmax-1)/2)
870 write(*,*)
'CalcDred err_pv2: ' &
871 ,err_pv2(rmax) , err_d0 * max(1d0,w_pv2**rmax,v_pv2**((rmax-1)/2), &
872 w_pv2*v_pv2**((rmax-1)/2) ) &
873 , err_c0 * z_pv2*max(w_pv2**(rmax+1), &
874 w_pv2*v_pv2**((rmax-1)/2), &
875 v_pv2**((rmax+1)/2)) &
876 , err_c(rmax-1) * max(z_pv2*w_pv2, &
877 z_pv2*w_pv2**2,z_pv2*v_pv2)
878 write(*,*)
'CalcDred err_pv2: ', &
879 err_c0 * z_pv2*w_pv2**(rmax+1), &
880 err_c0 * z_pv2* w_pv2*v_pv2**((rmax-1)/2), &
881 err_c0 * z_pv2* v_pv2**((rmax+1)/2) , &
891 err_pv2(rmax) = max( err_d0 * max(hw_pv2**rmax,v_pv2**(rmax/2)), &
892 err_c0 * z_pv2 * max(w_pv2*hw_pv2**(rmax), &
893 v_pv2**(rmax/2),w_pv2*v_pv2**(rmax/2)), &
894 err_c(rmax-1) * z_pv2 * max(hw_pv2,hw_pv2*w_pv2, v_pv2) )
913 err_pv(rmax) = err_pv(rmax)/impest_d
914 err_pv2(rmax) = err_pv2(rmax)/impest_d
920 write(*,*)
'CalcDred: err_pv',err_pv(rmax),err_pv2(rmax),err_req_d
921 write(*,*)
'CalcDred: acc_pv',err_pv(rmax)/d0est,err_pv2(rmax)/d0est,acc_req_d
937 if(use_pv.or.use_pv2)
then
939 if (min(err_pv(rmax),err_pv2(rmax)).le.err_req_d)
then
941 if (err_pv(rmax).le.err_pv2(rmax))
then
944 write(*,*)
'CalcDred: call Dpv 1 ',rmax,id,err_pv(rmax)
948 call calcdpv1(d,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,derr1,derr2)
955 dcount(1) = dcount(1)+1
956 drcalc(0:rmax) = drcalc(0:rmax)+1
962 checkest=derr(rmax)/err_pv(rmax)
963 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
964 write(*,*)
'CalcDred: estimate err_pv imprecise',err_pv(rmax),derr(rmax)
973 write(*,*)
'CalcDred: call Dpv2 1',rmax,id,err_pv2(rmax)
977 call calcdpv2(d,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,id,derr1,derr2)
983 dcount(2) = dcount(2)+1
984 drcalc(0:rmax)=drcalc(0:rmax)+2
989 checkest=derr(rmax)/err_pv2(rmax)
990 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
991 write(*,*)
'CalcDred: estimate err_pv2 imprecise',err_pv2(rmax),derr(rmax)
1001 err_d0 = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
1002 err_req_d = acc_req_d * abs(d(0,0,0,0))
1007 dtyp = max(abs(d(0,0,0,0)), &
1008 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
1010 dtyp = abs(d(0,0,0,0))
1012 if(dtyp.eq.0d0) dtyp = d0est
1013 err_req_d = acc_req_d * dtyp
1017 write(*,*)
'CalcDred Derr1 after PV = ',derr1
1019 write(*,*)
'CalcDred Dacc1 after PV = ',derr1/dtyp
1020 write(*,*)
'CalcDred err1_D',derr1(rmax)
1021 write(*,*)
'CalcDred Derr2 after PV = ',derr2
1023 write(*,*)
'CalcDred Dacc2 after PV = ',derr2/dtyp
1024 write(*,*)
'CalcDred err2_D',derr2(rmax)
1028 if (derr1(rmax).lt.err_req_d)
then
1029 dcount(dcalc+dcountoffset0) = dcount(dcalc+dcountoffset0)+1
1063 if (
fac_g.ge.1)
then
1066 err_g_c = err_c(rmax)
1073 err_g_cr = max( err_c(rmax), err_c0 * u_g**rmax ) * z_g
1075 err_g_exp = u_g**(rmax-1) * dtyp
1082 err_g_c = err_c(rmax)
1089 write(*,*)
'CalcDred: after Gram pars',use_g,
fac_g,x_g,u_g,z_g,err_g_cr,err_c(rmax),err_c0,err_g_exp
1091 write(*,*)
'CalcDred: after Gram pars',err_c(rmax), err_c0 * u_g**rmax
1109 err_gm_exp = err_inf
1110 err_gm_c = err_c(rmax)
1115 err_gm_cr = max( err_c(rmax), err_c0 * u_gm**rmax ) * z_gm
1116 err_gm_c = err_gm_cr
1117 err_gm_exp = u_gm**(rmax-1) * dtyp
1123 err_gm_exp = err_inf
1124 err_gm_c = err_c(rmax)
1128 write(*,*)
'CalcDred: after mod Gram pars',use_gm,
fac_gm,
x_gm,u_gm,z_gm,err_gm_cr,err_c(rmax),err_c0,err_gm_exp
1130 write(*,*)
'CalcDred: after mod Gram pars',err_c(rmax), err_c0 * u_gm**rmax
1138 err_gm_exp = err_inf
1139 err_gm_c = err_c(rmax)
1148 v1_gy = max(1d0,v_gy)
1149 fac_gy = max(x_gy,y_gy)*v1_gy
1154 err_gy_exp = err_inf
1155 err_gy_c = err_c(rmax+1)
1162 err_gy_cr = max( err_c(rmax) * v1_gy, err_c(rmax+1) )
1163 err_gy_c = err_gy_cr * b_gy
1164 err_gy_exp = 1d0 * dtyp
1170 err_gy_exp = err_inf
1171 err_gy_c = err_c(rmax+1)
1179 write(*,*)
'CalcDred: after GramCay pars',use_gy,
fac_gy,x_gy,y_gy,v_gy,b_gy,err_gy_cr,err_gy_exp
1189 v_gp = abs(mm02/
fmax)
1190 v1_gp = max(1d0,v_gp)
1196 err_gp_exp = err_inf
1197 err_gp_c = err_c(rmax)
1204 err_gp_cr = max(err_c0 * v_gp**rmax , err_c(rmax)) * z_gp
1205 err_gp_c = err_gp_cr
1206 err_gp_exp = v1_gp**(rmax-1) * dtyp
1212 err_gp_exp = err_inf
1213 err_gp_c = err_c(rmax)
1221 write(*,*)
'CalcDred: after Mom pars',use_gp,
fac_gp,w_gp,v_gp,z_gp,err_gp_cr,err_gp_exp
1229 y1_gr = max(1d0,y_gr)
1234 if (
fac_gr.ge.1.or.2*rmax.gt.rmax_c)
then
1236 err_gr_exp = err_inf
1237 err_gr_c = err_c(rmax)
1241 err_gr_cr = err_c(rmax)
1242 err_gr_c = err_gr_cr * a_gr
1243 err_gr_exp = y1_gr * dtyp
1249 err_gr_exp = err_inf
1250 err_gr_c = err_c(rmax)
1258 write(*,*)
'CalcDred: after revGram pars',use_gr,
fac_gr,x_gr,y_gr,y1_gr,a_gr,err_gr_cr,err_c(rmax),err_c0,err_gr_exp
1259 write(*,*)
'CalcDred: after revGram pars',err_gr_c,dtyp
1265 if (abs(m02).gt.
m2scale*dprec_cll)
then
1266 x_gpf =
fmax/abs(m02)
1267 y_gpf =
maxz/abs(m02)
1269 v1_gpf = max(1d0,v_gpf)
1275 err_gpf_exp = err_inf
1276 err_gpf_c = err_c(rmax+1)
1281 b_gpf = 1d0/abs(m02)
1282 err_gpf_cr = max( err_c(rmax), err_c(rmax+1) )
1283 err_gpf_c = err_gpf_cr * b_gpf
1284 err_gpf_exp = 1d0 * dtyp
1290 err_gpf_exp = err_inf
1291 err_gpf_c = err_c(rmax+1)
1298 write(*,*)
'CalcDred: after pf pars',use_gpf,
fac_gpf,x_gpf,y_gpf,v_gpf,b_gpf,err_gpf_cr,err_gpf_exp,err_gpf
1305 if(use_pv.or.use_pv2.or.use_g.or.use_gy.or.use_gp.or.use_gr.or.use_gm.or.use_gpf.eqv..false.)
then
1306 call seterrflag_coli(-6)
1307 call errout_coli(
'CalcDred',
' no reduction method works', &
1310 if (errorwriteflag)
then
1311 write(nerrout_coli,fmt10)
' CalcDred: p10 = ',p10
1312 write(nerrout_coli,fmt10)
' CalcDred: p21 = ',p21
1313 write(nerrout_coli,fmt10)
' CalcDred: p32 = ',p32
1314 write(nerrout_coli,fmt10)
' CalcDred: p30 = ',p30
1315 write(nerrout_coli,fmt10)
' CalcDred: p20 = ',p20
1316 write(nerrout_coli,fmt10)
' CalcDred: p31 = ',p31
1317 write(nerrout_coli,fmt10)
' CalcDred: m02 = ',m02
1318 write(nerrout_coli,fmt10)
' CalcDred: m12 = ',m12
1319 write(nerrout_coli,fmt10)
' CalcDred: m22 = ',m22
1320 write(nerrout_coli,fmt10)
' CalcDred: m32 = ',m32
1328 write(*,*)
'CalcDred: exit'
1351 if (err_g_exp.gt.err_g_c)
then
1353 err_g_exp = err_g_exp*
fac_g
1354 err_g_c = max(err_g_cr,err_c(rmax+g)*z_g*x_g**g)
1355 err_g(rmax) = max(err_g_exp,err_g_c)
1356 if(err_g(rmax).lt.err_req_d)
then
1359 g = min(max(g+2,3*g/2),rmax_d-rmax)
1373 if (err_gm_exp.gt.err_gm_c)
then
1375 err_gm_exp = err_gm_exp*
fac_gm
1376 err_gm_c = max(err_gm_cr,err_c(rmax+gm)*z_gm*
x_gm**gm)
1377 err_gm(rmax) = max(err_gm_exp,err_gm_c)
1378 if(err_gm(rmax).lt.err_req_d)
then
1381 gm = min(max(gm+2,3*gm/2),rmax_d-rmax)
1395 if (mod(i,2).eq.1)
then
1398 write(*,*)
'CalcDred: it gy',use_gy,err_gy_exp,err_gy_c,err_gy(rmax),err_req_d
1402 if (err_gy_exp.gt.err_gy_c.and.err_gy(rmax).gt.err_req_d)
then
1404 err_gy_exp = err_gy_exp*
fac_gy
1405 err_gy_c = b_gy*max(err_gy_cr, &
1406 max(err_c(rmax+2*gy)*v1_gy,err_c(rmax+2*gy+1))*y_gy**gy, &
1407 max(err_c(rmax+gy)*v1_gy,err_c(rmax+gy+1))*(max(x_gy,v_gy*y_gy))**gy)
1408 err_gy(rmax) = max(err_gy_exp,err_gy_c)
1411 write(*,*)
'CalcDred i gy',i,gy,err_gy_exp,err_gy_c,err_gy(rmax)
1412 write(*,*)
'CalcDred i ',err_gy_cr, &
1413 max(err_c(rmax+2*gy)*v1_gy,err_c(rmax+2*gy+1))*y_gy**gy, &
1414 max(err_c(rmax+gy)*v1_gy,err_c(rmax+gy+1))*(max(x_gy,v_gy*y_gy))**gy
1415 write(*,*)
'CalcDred i ', b_gy*err_c(rmax+2*gy)*v1_gy*y_gy**gy, &
1416 b_gy*err_c(rmax+2*gy+1)*y_gy**gy
1417 write(*,*)
'CalcDred i ', &
1418 b_gy,err_c(rmax+2*gy+1),y_gy**gy
1419 write(*,*)
'CalcDred i ', b_gy*x_gy**gy*err_c(rmax+gy)*v1_gy, &
1420 b_gy*err_c(rmax+gy+1)*x_gy**gy
1421 write(*,*)
'CalcDred i ', &
1422 b_gy,err_c(rmax+gy+1),x_gy**gy,x_gy,gy
1425 if(err_gy(rmax).lt.err_req_d)
then
1428 gy = min(max(gy+2,2*gy),(rmax_d-rmax)/2)
1443 if (err_gp_exp.gt.err_gp_c.and.err_gp(rmax).gt.err_req_d)
then
1445 err_gp_exp = err_gp_exp*
fac_gp
1446 err_gp_c = max(err_c(rmax+gp)*z_gp*w_gp**gp,err_gp_cr)
1447 err_gp(rmax) = max(err_gp_exp,err_gp_c)
1448 if(err_gp(rmax).lt.err_req_d)
then
1451 gp = min(max(gp+2,3*gp/2),rmax_d-rmax)
1459 if (mod(i,2).eq.1.and.i.le.rmax_c-2*rmax)
then
1467 if (err_gr_exp.gt.err_gr_c.and.err_gr(rmax).gt.err_req_d)
then
1469 err_gr_exp = err_gr_exp*
fac_gr
1470 err_gr_c = a_gr*max(err_gr_cr, &
1471 max(err_c(rmax+gr),err_c(rmax+gr+1)*y_gr)*
fac_gr**gr)
1472 err_gr(rmax) = max(err_gr_exp,err_gr_c)
1476 if(err_gr(rmax).lt.err_req_d)
then
1481 gr = min(max(gr+2,2*gr),rmax_d-rmax,max(0,(rmax_c-2*rmax)/2))
1491 if (mod(i,2).eq.1)
then
1494 write(*,*)
'CalcDred: it gy',use_gy,err_gy_exp,err_gy_c,err_gy(rmax),err_req_d
1498 if (err_gpf_exp.gt.err_gpf_c.and.err_gpf(rmax).gt.err_req_d)
then
1500 err_gpf_exp = err_gpf_exp*
fac_gpf
1501 err_gpf_c = b_gpf*max(err_gpf_cr, &
1502 max(err_c(rmax+2*gpf)*v1_gpf,err_c(rmax+2*gpf+1))*y_gpf**gpf, &
1503 max(err_c(rmax+gpf)*v1_gpf,err_c(rmax+gpf+1))*(max(x_gpf,v_gpf*y_gpf))**gpf)
1504 err_gpf(rmax) = max(err_gpf_exp,err_gpf_c)
1507 write(*,*)
'CalcDred i gpf',i,gpf,err_gpf_exp,err_gpf_c,err_gpf(rmax)
1508 write(*,*)
'CalcDred i ',err_gpf_cr, &
1509 max(err_c(rmax+2*gpf)*v1_gpf,err_c(rmax+2*gpf+1))*y_gpf**gpf, &
1510 max(err_c(rmax+gpf)*v1_gpf,err_c(rmax+gpf+1))*(max(x_gpf,v_gpf*y_gpf))**gpf
1511 write(*,*)
'CalcDred i ', b_gpf*err_c(rmax+2*gpf)*v1_gpf*y_gpf**gpf, &
1512 b_gpf*err_c(rmax+2*gpf+1)*y_gpf**gpf
1513 write(*,*)
'CalcDred i ', &
1514 b_gpf,err_c(rmax+2*gpf+1),y_gpf**gpf
1515 write(*,*)
'CalcDred i ', b_gpf*x_gpf**gpf*err_c(rmax+gpf)*v1_gpf, &
1516 b_gpf*err_c(rmax+gpf+1)*x_gpf**gpf
1517 write(*,*)
'CalcDred i ', &
1518 b_gpf,err_c(rmax+gpf+1),x_gpf**gpf,x_gpf,gpf
1521 if(err_gpf(rmax).lt.err_req_d)
then
1524 gpf = min(max(gpf+2,2*gpf),(rmax_d-rmax)/2)
1541 err_g(rmax) = err_g(rmax)/impest_d
1542 err_gy(rmax) = err_gy(rmax)/impest_d
1543 err_gp(rmax) = err_gp(rmax)/impest_d
1544 err_gr(rmax) = err_gr(rmax)/impest_d
1545 err_gm(rmax) = err_gm(rmax)/impest_d
1546 err_gpf(rmax) = err_gpf(rmax)/impest_d
1549 write(*,*)
'iexp=',iexp
1551 write(*,*)
'errexp=',err_g_exp,err_gy_exp,err_gp_exp,err_gr_exp,err_gm_exp,err_gpf_exp,err_req_d
1552 write(*,*)
'errexptot=',i,g,err_g(rmax),gy,err_gy(rmax),gp,err_gp(rmax), &
1554 gr,err_gr(rmax),gm,err_gm(rmax),gpf,err_gpf(rmax)
1555 write(*,*)
'accexptot=',i,g,err_g(rmax)/dtyp,gy,err_gy(rmax)/dtyp, &
1556 gp,err_gp(rmax)/dtyp,gr,err_gr(rmax)/dtyp,gm,err_gm(rmax)/dtyp, &
1557 gpf,err_gpf(rmax)/dtyp
1577 call calcdg(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,g,g,id,derr1_alt,derr2_alt)
1579 derr_alt = derr2_alt
1581 derr_alt = derr1_alt
1583 dcount(3) = dcount(3)+1
1584 drcalc(0:rmax)=drcalc(0:rmax)+4
1586 drmethod_alt(0:rmax)=4
1589 checkest=derr_alt(rmax)/err_g(rmax)
1590 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1591 write(*,*)
'CalcDred: estimate err_g imprecise ',err_g(rmax),derr_alt(rmax)
1597 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1601 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1603 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1604 write(*,*)
'CalcDred after exp method=',drmethod
1608 call calcdgy(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,gy,gy,id,derr1_alt,derr2_alt)
1610 derr_alt = derr2_alt
1612 derr_alt = derr1_alt
1614 dcount(4) = dcount(4)+1
1615 drcalc(0:rmax)=drcalc(0:rmax)+8
1617 drmethod_alt(0:rmax)=8
1620 checkest=derr_alt(rmax)/err_gy(rmax)
1621 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1622 write(*,*)
'CalcDred: estimate err_gy imprecise',err_gy(rmax),derr_alt(rmax),checkest
1627 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1630 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1631 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1632 write(*,*)
'CalcDred after exp method=',drmethod
1640 call calcdgp(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,gp,gp,id,derr1_alt,derr2_alt)
1642 derr_alt = derr2_alt
1644 derr_alt = derr1_alt
1646 dcount(5) = dcount(5)+1
1647 drcalc(0:rmax)=drcalc(0:rmax)+16
1649 drmethod_alt(0:rmax)=16
1652 checkest=derr_alt(rmax)/err_gp(rmax)
1653 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1654 write(*,*)
'CalcDred: estimate err_gp imprecise',err_gp(rmax),derr_alt(rmax)
1659 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1662 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1663 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1664 write(*,*)
'CalcDred after exp method=',drmethod
1668 call calcdgr(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,gr,gr,id,derr1_alt,derr2_alt)
1670 derr_alt = derr2_alt
1672 derr_alt = derr1_alt
1674 dcount(6) = dcount(6)+1
1675 drcalc(0:rmax)=drcalc(0:rmax)+32
1677 drmethod_alt(0:rmax)=32
1680 checkest=derr_alt(rmax)/err_gr(rmax)
1681 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1682 write(*,*)
'CalcDred: estimate err_gr imprecise',err_gr(rmax),derr_alt(rmax)
1687 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1690 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1691 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1692 write(*,*)
'CalcDred after exp method=',drmethod
1697 call calcdgm(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,gm,gm,id,derr1_alt,derr2_alt)
1699 derr_alt = derr2_alt
1701 derr_alt = derr1_alt
1703 dcount(7) = dcount(7)+1
1704 drcalc(0:rmax)=drcalc(0:rmax)+64
1706 drmethod_alt(0:rmax)=64
1709 checkest=derr_alt(rmax)/err_gm(rmax)
1710 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1711 write(*,*)
'CalcDred: estimate err_gm imprecise',err_gm(rmax),derr_alt(rmax)
1716 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1719 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1720 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1721 write(*,*)
'CalcDred after exp method=',drmethod
1726 call calcdgpf(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,rmax,gpf,gpf,id,derr1_alt,derr2_alt)
1728 derr_alt = derr2_alt
1730 derr_alt = derr1_alt
1732 dcount(7) = dcount(7)+1
1733 drcalc(0:rmax)=drcalc(0:rmax)+64
1735 drmethod_alt(0:rmax)=64
1738 checkest=derr_alt(rmax)/err_gpf(rmax)
1739 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
1740 write(*,*)
'CalcDred: estimate err_gpf imprecise',err_gpf(rmax),derr_alt(rmax),checkest
1745 call copydimp3(d,d_alt,derr,derr_alt,derr1,derr1_alt,derr2,derr2_alt,drmethod,drmethod_alt,rmax,rmax)
1748 write(*,*)
'CalcDred after exp Derr=',derr,err_req_d
1749 write(*,*)
'CalcDred after exp Dacc=',derr/dtyp
1750 write(*,*)
'CalcDred after exp method=',drmethod
1760 if(.not.lerr_d0.and.iexp.ne.0)
then
1762 err_d0 = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
1771 dtyp = max(abs(d(0,0,0,0)), &
1772 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
1774 dtyp = abs(d(0,0,0,0))
1776 err_req_d = acc_req_d * dtyp
1779 write(*,*)
'CalcDred ',rmax,derr1(rmax),err_req_d
1782 if (derr1(rmax).le.err_req_d)
then
1783 dcount(dcalc+dcountoffset0) = dcount(dcalc+dcountoffset0)+1
1793 write(*,*)
'CalcDred no optimal method'
1794 write(*,*)
'err_req_D=',err_req_d
1795 write(*,*)
'err_est=',err_pv(rmax),err_pv2(rmax),err_g(rmax),err_gy(rmax),err_gp(rmax),err_gr(rmax),err_gpf(rmax)
1802 if(use_pv.and.mod(drcalc(r),2).ne.1)
then
1806 if (mod(r,2).eq.1)
then
1807 err_pv(r) = max( w_pv**((r-1)/2) * v_pv * err_d0, &
1808 w_pv**((r-1)/2) * z_pv * err_c0, z_pv * err_c(r-1) )
1815 else if (r.ne.0)
then
1816 err_pv(r) = max( w_pv**(r/2) * err_d0, &
1817 w_pv**(r/2-1) * v_pv * z_pv * err_c0, z_pv * err_c(r-1) )
1827 err_pv(r) = err_pv(r)/impest_d
1830 if (use_pv2.and.mod(drcalc(r),4)-mod(drcalc(r),2).ne.2)
then
1836 if (mod(r,2).eq.1)
then
1844 err_pv2(r) = max( err_d0 * max(hw_pv2**r, &
1845 hw_pv2*v_pv2**((r-1)/2) ), &
1846 err_c0 * z_pv2* max(w_pv2*w_pv2**(r), &
1847 hw_pv2*v_pv2**((r-1)/2), &
1848 w_pv2*hw_pv2*v_pv2**((r-1)/2), &
1849 v_pv2**((r+1)/2)), &
1850 err_c(r-1) * z_pv2 * max(hw_pv2,w_pv2*hw_pv2**2,v_pv2) )
1863 else if (r.ne.0)
then
1869 err_pv2(r) = max( err_d0 * max(hw_pv2**r,v_pv2**(r/2)), &
1870 err_c0 * z_pv2 * max(w_pv2*hw_pv2**(r), &
1871 v_pv2**(r/2),w_pv2*v_pv2**(r/2)), &
1872 err_c(r-1) * z_pv2 * max(hw_pv2, w_pv2*hw_pv2**2, v_pv2) )
1887 err_pv2(r) = err_pv2(r)/impest_d
1890 if (mod(drcalc(r),8)-mod(drcalc(r),4).ne.4.and.use_g)
then
1892 err_g_cr = max( err_c(r), err_c0 * u_g**r ) * z_g
1894 err_g_exp = u_g**(r-1) * dtyp
1899 err_g_exp = err_g_exp*
fac_g
1900 err_g_c = max(err_g_cr,err_c(r+g)*z_g*x_g**g)
1901 err_g(r) = max(err_g_exp,err_g_c)
1902 if (err_g_exp.lt.err_g_c.or.err_g(r).lt.err_req_d)
exit
1905 g = min(max(g+2,2*g),rmax_d-r)
1907 err_g(r) = err_g(r)/impest_d
1910 if (mod(drcalc(r),16)-mod(drcalc(r),8).ne.8.and.use_gy)
then
1912 err_gy_cr = max( err_c(r) * v1_gy, err_c(r+1) )
1913 err_gy_c = err_gy_cr * b_gy
1914 err_gy_exp = 1d0 * dtyp
1919 if (mod(i,2).eq.1)
then
1921 err_gy_exp = err_gy_exp*
fac_gy
1922 err_gy_c = b_gy*max(err_gy_cr, &
1923 max(err_c(r+2*gy)*v1_gy,err_c(r+2*gy+1))*y_gy**gy, &
1924 max(err_c(r+gy)*v1_gy,err_c(r+gy+1))*(max(x_gy,v_gy*y_gy))**gy)
1925 err_gy(r) = max(err_gy_exp,err_gy_c)
1926 if (err_gy_exp.lt.err_gy_c.or.err_gy(r).lt.err_req_d)
exit
1930 gy = min(max(gy+2,2*gy),(rmax_d-r)/2)
1932 err_gy(r) = err_gy(r)/impest_d
1935 if (mod(drcalc(r),32)-mod(drcalc(r),16).ne.16.and.use_gp)
then
1937 err_gp_cr = max(err_c0*v_gp**r,err_c(r))*z_gp
1938 err_gp_exp = v1_gp**(r-1) * dtyp
1943 err_gp_exp = err_gp_exp*
fac_gp
1944 err_gp_c = max(err_c(r+gp)*z_gp*w_gp**gp,err_gp_cr)
1945 err_gp(r) = max(err_gp_exp,err_gp_c)
1946 if (err_gp_exp.lt.err_gp_c.or.err_gp(r).lt.err_req_d)
exit
1949 gp = min(max(gp+2,3*gp/2),rmax_d-r)
1951 err_gp(r) = err_gp(r)/impest_d
1954 if (mod(drcalc(r),64)-mod(drcalc(r),32).ne.32.and.use_gr)
then
1956 err_gr_cr = err_c(r)
1957 err_gr_c = err_gr_cr * a_gr
1958 err_gr_exp = y1_gr * dtyp
1962 do i=0,min(rmax_d-r,rmax_c-2*r)
1963 if (mod(i,2).eq.1)
then
1965 err_gr_exp = err_gr_exp*
fac_gr
1966 err_gr_c = a_gr*max(err_gr_cr, &
1967 max(err_c(r+gr),err_c(r+gr+1)*y_gr)*
fac_gr**gr)
1968 err_gr(r) = max(err_gr_exp,err_gr_c)
1971 write(*,*)
'CalcDgr err_gr',i,gr,err_gr_exp,err_gr_c,err_gr(r),err_req_d
1974 if (err_gr_exp.lt.err_gr_c.or.err_gr(r).lt.err_req_d)
exit
1980 gr = min(max(gr+2,2*gr),rmax_d-r,max(0,(rmax_c-2*r)/2))
1982 err_gr(r) = err_gr(r)/impest_d
1986 if (mod(drcalc(r),128)-mod(drcalc(r),64).ne.64.and.use_gpf)
then
1988 err_gpf_cr = max( err_c(r) * v1_gpf, err_c(r+1) )
1989 err_gpf_c = err_gpf_cr * b_gpf
1990 err_gpf_exp = 1d0 * dtyp
1995 if (mod(i,2).eq.1)
then
1997 err_gpf_exp = err_gpf_exp*
fac_gpf
1998 err_gpf_c = b_gpf*max(err_gpf_cr, &
1999 max(err_c(r+2*gpf)*v1_gpf,err_c(r+2*gpf+1))*y_gpf**gpf, &
2000 max(err_c(r+gpf)*v1_gpf,err_c(r+gpf+1))*(max(x_gpf,v_gpf*y_gpf))**gpf)
2001 err_gpf(r) = max(err_gpf_exp,err_gpf_c)
2002 if (err_gpf_exp.lt.err_gpf_c.or.err_gpf(r).lt.err_req_d)
exit
2006 gpf = min(max(gpf+2,2*gpf),(rmax_d-r)/2)
2008 err_gpf(r) = err_gpf(r)/impest_d
2014 if (mod(drcalc(r),128)-mod(drcalc(r),64).ne.64.and.use_gm)
then
2016 err_gm_cr = max( err_c(r), err_c0 * u_gm**r ) * z_gm
2017 err_gm_c = err_gm_cr
2018 err_gm_exp = u_gm**(r-1) * dtyp
2023 err_gm_exp = err_gm_exp*
fac_gm
2024 err_gm_c = max(err_gm_cr,err_c(r+gm)*z_gm*
x_gm**gm)
2025 err_gm(r) = max(err_gm_exp,err_gm_c)
2026 if (err_gm_exp.lt.err_gm_c.or.err_gm(r).lt.err_req_d)
exit
2029 gm = min(max(gm+2,2*gm),rmax_d-r)
2031 err_gm(r) = err_gm(r)/impest_d
2037 write(*,*)
'CalcDred: bef final loop ord methods',r,g,gy,gp,gr,gm,gpf
2038 write(*,*)
'CalcDred: bef final loop err methods',r,err_pv(r),err_pv2(r), &
2039 err_g(r),err_gy(r),err_gp(r),err_gr(r),err_gm(r),err_gpf(r)
2040 write(*,*)
'CalcDred: bef final loop acc methods',r,err_pv(r)/dtyp,err_pv2(r)/dtyp, &
2041 err_g(r)/dtyp,err_gy(r)/dtyp,err_gp(r)/dtyp, &
2042 err_gr(r)/dtyp,err_gm(r)/dtyp,err_gpf(r)/dtyp
2043 write(*,*)
'CalcDred: bef final loop method',r,drcalc(r),drmethod(r)
2047 if (min(err_pv(r),err_pv2(r)).le.min(err_g(r),err_gy(r),err_gp(r),err_gr(r),err_gpf(r)) &
2048 .and.min(err_pv(r),err_pv2(r)).lt.err_inf)
then
2050 if (use_pv.and.err_pv(r).le.err_pv2(r).and.mod(drcalc(r),2).ne.1)
then
2064 write(*,*)
'CalcDred: call Dpv 2',r,id
2073 call calcdpv1(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,id,derr1_alt,derr2_alt)
2075 call calcdpv1(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2076 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,id,derr1_alt(0:r),derr2_alt(0:r))
2079 derr_alt = derr2_alt
2081 derr_alt = derr1_alt
2083 dcount(11) = dcount(11)+1
2084 drcalc(0:r)=drcalc(0:r)+1
2087 checkest=derr_alt(r)/err_pv(r)
2090 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2091 write(*,*)
'CalcDred: estimate err_pv imprecise',err_pv(r),derr_alt(r)
2096 write(*,*)
'final loop r Dpv D(1,0,1,0)',r,d_alt(1,0,1,0),d(1,0,1,0)
2097 write(*,*)
'final loop r Dpv Derr',derr_alt(2),derr(2)
2099 err_pv(0:r)=derr_alt(0:r)
2101 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2102 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2105 dtyp = max(abs(d(0,0,0,0)), &
2106 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2108 dtyp = abs(d(0,0,0,0))
2110 err_req_d = acc_req_d * dtyp
2113 write(*,*)
'CalcDred: after pv 2nd try Dmethod_alt=',drmethod_alt
2114 write(*,*)
'CalcDred: after pv 2nd try Derr_alt(r)=',derr_alt
2115 write(*,*)
'CalcDred: after pv 2nd try Dacc_alt(r)=',derr_alt/dtyp
2116 write(*,*)
'CalcDred: after pv 2nd try Dmethod=',drmethod
2117 write(*,*)
'CalcDred: after pv 2nd try Derr(r)=',derr
2118 write(*,*)
'CalcDred: after pv 2nd try Dacc(r)=',derr/dtyp
2126 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2128 elseif (use_pv2.and.err_pv2(r).le.err_pv(r).and.mod(drcalc(r),4)-mod(drcalc(r),2).ne.2)
then
2142 write(*,*)
'CalcDred: call Dpv2 2',r,id
2145 call calcdpv2(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,id,derr1_alt,derr2_alt)
2147 call calcdpv2(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2148 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,id,derr1_alt(0:r),derr2_alt(0:r))
2152 derr_alt = derr2_alt
2154 derr_alt = derr1_alt
2156 dcount(12) = dcount(12)+1
2157 drcalc(0:r)=drcalc(0:r)+2
2160 checkest=derr_alt(r)/err_pv2(r)
2163 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2164 write(*,*)
'CalcDred: estimate err_pv2 imprecise',err_pv2(r),derr_alt(r)
2167 err_pv2(0:r)=derr_alt(0:r)
2171 write(*,*)
'final loop r Dpv2 D(1,0,1,0)',r,d_alt(1,0,1,0),d(1,0,1,0)
2172 write(*,*)
'final loop r Dpv2 Derr',derr_alt(2),derr(2)
2175 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2176 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2179 dtyp = max(abs(d(0,0,0,0)), &
2180 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2182 dtyp = abs(d(0,0,0,0))
2184 err_req_d = acc_req_d * dtyp
2186 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2189 write(*,*)
'CalcDred: after pv 2nd try Dmethod=',drmethod
2190 write(*,*)
'CalcDred: after pv 2nd try Derr(r)=',derr
2191 write(*,*)
'CalcDred: after pv 2nd try Dacc(r)=',derr/dtyp
2199 write(*,*)
'CalcDred: explore exps once more'
2202 if (use_g.and.err_g(r).le.min(err_gy(r),err_gp(r),err_gr(r),err_gpf(r)) &
2203 .and.mod(drcalc(r),8)-mod(drcalc(r),4).ne.4)
then
2217 call calcdg(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,g,rmax_d,id,derr1_alt,derr2_alt)
2219 call calcdg(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2220 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,g,rmax_d,id,derr1_alt(0:r),derr2_alt(0:r))
2223 derr_alt = derr2_alt
2225 derr_alt = derr1_alt
2227 dcount(13) = dcount(13)+1
2228 drcalc(0:r)=drcalc(0:r)+4
2231 checkest=derr_alt(r)/err_g(r)
2234 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2235 write(*,*)
'CalcDred: estimate err_g imprecise ',err_g(r),derr_alt(r)
2239 err_g(0:r)=derr_alt(0:r)
2241 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2242 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2245 dtyp = max(abs(d(0,0,0,0)), &
2246 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2248 dtyp = abs(d(0,0,0,0))
2250 err_req_d = acc_req_d * dtyp
2253 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2254 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2255 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2264 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2266 else if (use_gy.and.err_gy(r).le.min(err_g(r),err_gp(r),err_gr(r),err_gpf(r)) &
2267 .and.mod(drcalc(r),16)-mod(drcalc(r),8).ne.8)
then
2281 call calcdgy(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gy,(rmax_d)/2,id,derr1_alt,derr2_alt)
2283 call calcdgy(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2284 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gy,(rmax_d)/2,id,derr1_alt(0:r),derr2_alt(0:r))
2287 derr_alt = derr2_alt
2289 derr_alt = derr1_alt
2291 dcount(14) = dcount(14)+1
2292 drcalc(0:r)=drcalc(0:r)+8
2295 checkest=derr_alt(r)/err_gy(r)
2298 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2299 write(*,*)
'CalcDred: estimate err_gy imprecise',err_gy(r),derr_alt(r)
2303 err_gy(0:r)=derr_alt(0:r)
2305 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2306 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2309 dtyp = max(abs(d(0,0,0,0)), &
2310 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2312 dtyp = abs(d(0,0,0,0))
2314 err_req_d = acc_req_d * dtyp
2317 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2318 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2319 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2328 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2330 elseif (use_gp.and.err_gp(r).le.min(err_g(r),err_gy(r),err_gr(r),err_gpf(r)) &
2331 .and.mod(drcalc(r),32)-mod(drcalc(r),16).ne.16)
then
2345 call calcdgp(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gp,rmax_d,id,derr1_alt,derr2_alt)
2347 call calcdgp(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2348 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gp,rmax_d,id,derr1_alt(0:r),derr2_alt(0:r))
2351 derr_alt = derr2_alt
2353 derr_alt = derr1_alt
2355 dcount(15) = dcount(15)+1
2356 drcalc(0:r)=drcalc(0:r)+16
2358 drmethod_alt(0:r)=16
2359 checkest=derr_alt(r)/err_gp(r)
2362 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2363 write(*,*)
'CalcDred: estimate err_gp imprecise',err_gp(r),derr_alt(r)
2367 err_gp(0:r)=derr_alt(0:r)
2369 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2370 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2373 dtyp = max(abs(d(0,0,0,0)), &
2374 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2376 dtyp = abs(d(0,0,0,0))
2378 err_req_d = acc_req_d * dtyp
2381 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2382 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2383 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2386 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2388 elseif (use_gr.and.err_gr(r).le.min(err_g(r),err_gy(r),err_gp(r),err_gpf(r)) &
2389 .and.mod(drcalc(r),64)-mod(drcalc(r),32).ne.32)
then
2403 call calcdgr(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gr,rmax_d,id,derr1_alt,derr2_alt)
2405 call calcdgr(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2406 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gr,rmax_d,id,derr1_alt(0:r),derr2_alt(0:r))
2409 derr_alt = derr2_alt
2411 derr_alt = derr1_alt
2413 dcount(16) = dcount(16)+1
2414 drcalc(0:r)=drcalc(0:r)+32
2416 drmethod_alt(0:r)=32
2417 checkest=derr_alt(r)/err_gr(r)
2420 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2421 write(*,*)
'CalcDred: estimate err_gr imprecise',err_gr(r),derr_alt(r)
2425 err_gr(0:r)=derr_alt(0:r)
2427 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2428 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2431 dtyp = max(abs(d(0,0,0,0)), &
2432 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2434 dtyp = abs(d(0,0,0,0))
2436 err_req_d = acc_req_d * dtyp
2439 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2440 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2441 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2444 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2453 else if (use_gpf.and.err_gpf(r).le.min(err_g(r),err_gy(r),err_gp(r),err_gr(r)) &
2454 .and.mod(drcalc(r),128)-mod(drcalc(r),64).ne.64)
then
2457 call calcdgpf(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gpf,(rmax_d)/2,id,derr1_alt,derr2_alt)
2459 call calcdgpf(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2460 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,gpf,(rmax_d)/2,id,derr1_alt(0:r),derr2_alt(0:r))
2463 derr_alt = derr2_alt
2465 derr_alt = derr1_alt
2467 dcount(17) = dcount(17)+1
2468 drcalc(0:r)=drcalc(0:r)+64
2470 drmethod_alt(0:r)=64
2471 checkest=derr_alt(r)/err_gpf(r)
2474 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2475 write(*,*)
'CalcDred: estimate err_gpf imprecise',err_gpf(r),derr_alt(r)
2479 err_gpf(0:r)=derr_alt(0:r)
2480 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2481 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2484 dtyp = max(abs(d(0,0,0,0)), &
2485 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2487 dtyp = abs(d(0,0,0,0))
2489 err_req_d = acc_req_d * dtyp
2492 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2493 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2494 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2503 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2506 else if (use_gm.and.err_gm(r).le.min(err_gy(r),err_gp(r),err_gr(r),err_g(r)) &
2507 .and.mod(drcalc(r),128)-mod(drcalc(r),64).ne.64)
then
2521 call calcdgm(d_alt,duv,p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,g,rmax_d,id,derr1_alt,derr2_alt)
2523 call calcdgm(d_alt(0:r,0:r,0:r,0:r),duv_alt(0:r,0:r,0:r,0:r), &
2524 p10,p21,p32,p30,p20,p31,m02,m12,m22,m32,r,g,rmax_d,id,derr1_alt(0:r),derr2_alt(0:r))
2527 derr_alt = derr2_alt
2529 derr_alt = derr1_alt
2531 dcount(17) = dcount(17)+1
2532 drcalc(0:r)=drcalc(0:r)+64
2534 drmethod_alt(0:r)=64
2535 checkest=derr_alt(r)/err_gm(r)
2538 if(checkest.gt.1d2*impest_d.or.checkest.lt.1d-2*impest_d)
then
2539 write(*,*)
'CalcDred: estimate err_g imprecise ',err_gm(r),derr_alt(r)
2543 err_gm(0:r)=derr_alt(0:r)
2545 call copydimp3(d,d_alt(0:r,0:r,0:r,0:r),derr,derr_alt(0:r),derr1,derr1_alt(0:r), &
2546 derr2,derr2_alt(0:r),drmethod,drmethod_alt(0:r),rmax,r)
2549 dtyp = max(abs(d(0,0,0,0)), &
2550 abs(d(0,1,0,0)),abs(d(0,0,1,0)),abs(d(0,0,0,1)))
2552 dtyp = abs(d(0,0,0,0))
2554 err_req_d = acc_req_d * dtyp
2557 write(*,*)
'CalcDred: after exp 2nd try Dmethod=',drmethod
2558 write(*,*)
'CalcDred: after exp 2nd try Derr(r)=',derr
2559 write(*,*)
'CalcDred: after exp 2nd try Dacc(r)=',derr/dtyp
2568 if(checkest.gt.impest_d.and.mode_coli.lt.1)
goto 100
2577 if(.not.lerr_d0)
then
2579 err_d0 = acc_def_d0*max( abs(d(0,0,0,0)), 1d0/sqrt(
adetx) )
2595 write(*,*)
'CalcDred: final loop err methods',r,err_pv(r),err_pv2(r), &
2596 err_g(r),err_gy(r),err_gp(r),err_gr(r),err_gm(r),err_gpf(r)
2597 write(*,*)
'CalcDred: final loop acc methods',r,err_pv(r)/dtyp,err_pv2(r)/dtyp, &
2598 err_g(r)/dtyp,err_gy(r)/dtyp,err_gp(r)/dtyp, &
2599 err_gr(r)/dtyp,err_gm(r)/dtyp,err_gpf(r)/dtyp
2600 write(*,*)
'CalcDred: final loop method',r,drcalc(r),drmethod(r)
2605 norm = abs(d(0,0,0,0))
2611 norm = max(norm,abs(d(0,n1,n2,n3)))
2615 acc_d = derr(rmax)/norm
2617 dcount(dcalc+dcountoffset0) = dcount(dcalc+dcountoffset0)+1
2620 write(*,*)
'CalcDred final err_D=',derr
2621 write(*,*)
'CalcDred final acc_D=',derr/norm,critacc_coli
2622 write(*,*)
'CalcDred final method_D=',drmethod
2625 if (acc_d.gt.sqrt(reqacc_coli))
then
2626 dcount(dcalc+dcountoffset3) = dcount(dcalc+dcountoffset3)+1
2629 if (acc_d.gt.reqacc_coli)
then
2630 dcount(dcalc+dcountoffset1) = dcount(dcalc+dcountoffset1)+1
2633 if (acc_d.gt.critacc_coli)
then
2635 dcount(dcalc+dcountoffset2) = dcount(dcalc+dcountoffset2)+1
2638 write(*,*)
'CritPoint D',critacc_coli,acc_d
2639 write(*,*)
'CritPoint D',critpointcntd,maxcritpointd
2646 #ifdef CritPointsCOLI
2647 critpointcntd = critpointcntd + 1
2649 if (critpointcntd.le.maxcritpointd.and.monitoring)
then
2651 call critpointsout_coli(
'D_coli',acc_d)
2652 write(ncpout_coli,*)
'arguments of CalcDred_coli:'
2653 write(ncpout_coli,*)
'rank = ', rmax
2654 write(ncpout_coli,fmt1)
'p10 = ', p10
2655 write(ncpout_coli,fmt1)
'p21 = ', p21
2656 write(ncpout_coli,fmt1)
'p32 = ', p32
2657 write(ncpout_coli,fmt1)
'p30 = ', p30
2658 write(ncpout_coli,fmt1)
'p20 = ', p20
2659 write(ncpout_coli,fmt1)
'p31 = ', p31
2660 write(ncpout_coli,fmt1)
'm02 = ', m02
2661 write(ncpout_coli,fmt1)
'm12 = ', m12
2662 write(ncpout_coli,fmt1)
'm22 = ', m22
2663 write(ncpout_coli,fmt1)
'm32 = ', m32
2664 if (critpointcntd.eq.maxcritpointd)
then
2665 write(ncpout_coli,*)
2666 write(ncpout_coli,*)
2667 write(ncpout_coli,*)
2668 write(ncpout_coli,*)
'***********************************************************'
2669 write(ncpout_coli,*)
2670 write(ncpout_coli,*)
' Further output of bad D functions will be suppressed '
2677 write(*,*)
'CalcDred exit D(1,0,0,0)',r,d_alt(1,0,0,0),d(1,0,0,0)