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)
422 m2max = max(abs(mm02),abs(mm12),abs(
mm22))
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
514 xadj(2,2) = 2d0*mm02*
z(1,1) -
f(1)*
f(1)
515 xadj(1,1) = 2d0*mm02*
z(2,2) -
f(2)*
f(2)
516 xadj(2,1) = 2d0*mm02*
z(2,1) -
f(1)*
f(2)
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
616 ccount(0) = ccount(0)+1
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)
785 ccount(1) = ccount(1)+1
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)
811 ccount(2) = ccount(2)+1
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
856 ccount(ccalc+ccountoffset0) = ccount(ccalc+ccountoffset0)+1
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
957 v_gp = max(1d0,abs(mm02)/
fmax)
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
1253 ccount(3) = ccount(3)+1
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
1283 ccount(3) = ccount(3)+1
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
1313 ccount(4) = ccount(4)+1
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
1342 ccount(5) = ccount(5)+1
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
1371 ccount(6) = ccount(6)+1
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
1401 ccount(4) = ccount(4)+1
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
1446 ccount(ccalc+ccountoffset0) = ccount(ccalc+ccountoffset0)+1
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
1680 ccount(9) = ccount(9)+1
1683 crmethod_alt(0:rmax)=1
1684 if (cerr_alt(rmax).lt.cerr(rmax))
then
1685 ccount(8) = ccount(8)+1
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
1718 ccount(ccalc+ccountoffset0) = ccount(ccalc+ccountoffset0)+1
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
1964 ccount(11) = ccount(11)+1
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
2019 ccount(12) = ccount(12)+1
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
2084 ccount(13) = ccount(13)+1
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
2140 ccount(14) = ccount(14)+1
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
2196 ccount(15) = ccount(15)+1
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
2252 ccount(16) = ccount(16)+1
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
2308 ccount(17) = ccount(17)+1
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
2373 ccount(ccalc+ccountoffset0) = ccount(ccalc+ccountoffset0)+1
2375 if (maxval(acc_c(0:rdef)-sqrt(reqacc_coli)).gt.0)
then
2376 ccount(ccalc+ccountoffset3) = ccount(ccalc+ccountoffset3)+1
2379 if (maxval(acc_c(0:rdef)-reqacc_coli).gt.0)
then
2380 ccount(ccalc+ccountoffset1) = ccount(ccalc+ccountoffset1)+1
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
2390 ccount(ccalc+ccountoffset2) = ccount(ccalc+ccountoffset2)+1
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'