From 9a58493988754c41c40c43f96772b3402bac0c13 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 18 Jan 2008 15:34:15 +0000 Subject: [PATCH] Krylov iterations: Centralized check convergence routines, with some fixes for RGMRES. --- krylov/psb_dbicg.f90 | 2 +- krylov/psb_dcg.f90 | 2 +- krylov/psb_dcgs.f90 | 3 +- krylov/psb_dcgstab.F90 | 2 +- krylov/psb_dcgstabl.f90 | 10 +- krylov/psb_drgmres.f90 | 13 +-- krylov/psb_krylov_mod.f90 | 204 ++++++++++++++++---------------------- krylov/psb_zcgs.f90 | 2 +- krylov/psb_zcgstab.f90 | 10 +- krylov/psb_zrgmres.f90 | 16 +-- 10 files changed, 104 insertions(+), 160 deletions(-) diff --git a/krylov/psb_dbicg.f90 b/krylov/psb_dbicg.f90 index 0e5bcbaf..b4d1d509 100644 --- a/krylov/psb_dbicg.f90 +++ b/krylov/psb_dbicg.f90 @@ -215,7 +215,7 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) itx = 0 - call psb_init_conv(methdname,istop_,itrace_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) if (info /= 0) Then call psb_errpush(4011,name) goto 9999 diff --git a/krylov/psb_dcg.f90 b/krylov/psb_dcg.f90 index 7be51987..3746e4a3 100644 --- a/krylov/psb_dcg.f90 +++ b/krylov/psb_dcg.f90 @@ -202,7 +202,7 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) rho = dzero - call psb_init_conv(methdname,istop_,itrace_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) if (info /= 0) Then call psb_errpush(4011,name) goto 9999 diff --git a/krylov/psb_dcgs.f90 b/krylov/psb_dcgs.f90 index 52afedbd..9980aa19 100644 --- a/krylov/psb_dcgs.f90 +++ b/krylov/psb_dcgs.f90 @@ -192,8 +192,7 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) itx = 0 - - call psb_init_conv(methdname,istop_,itrace_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) if (info /= 0) Then call psb_errpush(4011,name) goto 9999 diff --git a/krylov/psb_dcgstab.F90 b/krylov/psb_dcgstab.F90 index 94854215..923842fc 100644 --- a/krylov/psb_dcgstab.F90 +++ b/krylov/psb_dcgstab.F90 @@ -223,7 +223,7 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) call psb_set_coher(ictxt,isvch) itx = 0 - call psb_init_conv(methdname,istop_,itrace_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) if (info /= 0) Then call psb_errpush(4011,name) goto 9999 diff --git a/krylov/psb_dcgstabl.f90 b/krylov/psb_dcgstabl.f90 index 5066ffab..f912abde 100644 --- a/krylov/psb_dcgstabl.f90 +++ b/krylov/psb_dcgstabl.f90 @@ -124,7 +124,7 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is Real(Kind(1.d0)), Pointer :: ww(:), q(:), r(:), rt0(:), p(:), v(:), & & s(:), t(:), z(:), f(:), gamma(:), gamma1(:), gamma2(:), taum(:,:), sigma(:) - Integer :: litmax, naux, mglob, it, itrace_,& + Integer :: itmax_, naux, mglob, it, itrace_,& & np,me, n_row, n_col, nl, err_act Logical, Parameter :: exchange=.True., noexchange=.False. Integer, Parameter :: irmax = 8 @@ -160,9 +160,9 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is endif if (present(itmax)) then - litmax = itmax + itmax_ = itmax else - litmax = 1000 + itmax_ = 1000 endif if (present(itrace)) then @@ -234,7 +234,7 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is call psb_set_coher(ictxt,isvch) - call psb_init_conv(methdname,istop_,itrace_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) if (info /= 0) Then call psb_errpush(4011,name) goto 9999 @@ -247,7 +247,7 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is !!$ if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),' restart: ',itx,it - if (itx >= litmax) exit restart + if (itx >= itmax_) exit restart it = 0 call psb_geaxpby(done,b,dzero,r,desc_a,info) diff --git a/krylov/psb_drgmres.f90 b/krylov/psb_drgmres.f90 index f881f56e..7e3ae5ed 100644 --- a/krylov/psb_drgmres.f90 +++ b/krylov/psb_drgmres.f90 @@ -254,6 +254,7 @@ subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist call psb_errpush(info,name) goto 9999 end if + if ((itrace_ > 0).and.(me==0)) call log_header(methdname) itx = 0 restart: do @@ -452,17 +453,7 @@ subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist if (itrace_ > 0) & & call log_conv(methdname,me,itx,1,errnum,errden,eps) - if (present(err)) then - if (errden /= dzero) then - err = errnum/errden - else - err = errnum - end if - end if - - if (present(iter)) iter = itx - if (errnum > eps*errden) & - & call end_log(methdname,me,itx,errnum,errden,eps) + call log_end(methdname,me,itx,errnum,errden,eps,err=err,iter=iter) deallocate(aux,h,c,s,rs,rst, stat=info) if (info == 0) call psb_gefree(v,desc_a,info) diff --git a/krylov/psb_krylov_mod.f90 b/krylov/psb_krylov_mod.f90 index 275ebfa3..e8f64f98 100644 --- a/krylov/psb_krylov_mod.f90 +++ b/krylov/psb_krylov_mod.f90 @@ -209,7 +209,7 @@ Module psb_krylov_mod integer, parameter, private :: bni_=1, rni_=2, ani_=3, xni_=4, bn2_=5, xn2_=6 integer, parameter, private :: errnum_=7, errden_=8, eps_=9, rn2_=10 - integer, parameter, private :: stopc_=1, trace_=2 + integer, parameter, private :: stopc_=1, trace_=2, itmax_=3 integer, parameter, private :: ivsz_=16 type psb_itconv_type private @@ -448,51 +448,87 @@ contains end subroutine psb_zkrylov - subroutine log_conv(methdname,me,itx,itrace,errnum,errden,eps) + subroutine log_header(methdname) use psb_base_mod - character(len=*), intent(in) :: methdname - integer, intent(in) :: me, itx, itrace - real(kind(1.d0)), intent(in) :: errnum, errden, eps - character(len=*), parameter :: fmt='(a,i4,3(2x,es10.4))' + implicit none + character(len=*), intent(in) :: methdname + character(len=*), parameter :: fmt='(a18,1x,a4,3(2x,a10))' + integer, parameter :: outlen=18 + character(len=len(methdname)) :: mname + character(len=outlen) :: outname + + mname = adjustl(trim(methdname)) + write(outname,'(a)') mname(1:min(len_trim(mname),outlen-1))//':' + write(*,fmt) adjustl(outname),'Iter','Conv. Ind.','Epsilon' - if (me == 0) then + end subroutine log_header + + + subroutine log_conv(methdname,me,itx,itrace,errnum,errden,eps) + use psb_base_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: me, itx, itrace + real(kind(1.d0)), intent(in) :: errnum, errden, eps + character(len=*), parameter :: fmt='(a18,1x,i4,3(2x,es10.4))' + integer, parameter :: outlen=18 + character(len=len(methdname)) :: mname + character(len=outlen) :: outname + + if ((mod(itx,itrace) == 0).and.(me==0)) then + mname = adjustl(trim(methdname)) + write(outname,'(a)') mname(1:min(len_trim(mname),outlen-1))//':' if (errden > dzero ) then - write(*,fmt) methdname//': ',itx,errnum/errden,eps + write(*,fmt) adjustl(outname),itx,errnum/errden,eps else - write(*,fmt) methdname//': ',itx,errnum,eps + write(*,fmt) adjustl(outname),itx,errnum,eps end if - Endif + endif end subroutine log_conv - subroutine end_log(methdname,me,itx,errnum,errden,eps) + subroutine log_end(methdname,me,it,errnum,errden,eps,err,iter) use psb_base_mod + implicit none character(len=*), intent(in) :: methdname - integer, intent(in) :: me, itx + integer, intent(in) :: me, it real(kind(1.d0)), intent(in) :: errnum, errden, eps + real(kind(1.d0)), optional, intent(out) :: err + integer, optional, intent(out) :: iter + character(len=*), parameter :: fmt='(a,2x,es10.4,1x,a,1x,i4,1x,a)' character(len=*), parameter :: fmt1='(a,3(2x,es10.4))' - if (me==0) then - write(*,fmt) methdname//' failed to converge to ',eps,& - & ' in ',ITX,' iterations. ' - if (errden > dzero) then - write(*,fmt1) 'Last iteration convergence indicator: ',& - & errnum/errden - else - write(*,fmt1) 'Last iteration convergence indicator: ',& - & errnum/errden + if (errden == dzero) then + if (errnum > eps) then + if (me==0) then + write(*,fmt) trim(methdname)//' failed to converge to ',eps,& + & ' in ',it,' iterations. ' + write(*,fmt1) 'Last iteration convergence indicator: ',& + & errnum + end if + end if + if (present(err)) err=errnum + else + if (errnum/errden > eps) then + if (me==0) then + write(*,fmt) trim(methdname)//' failed to converge to ',eps,& + & ' in ',it,' iterations. ' + write(*,fmt1) 'Last iteration convergence indicator: ',& + & errnum/errden + end if endif + if (present(err)) err=errnum/errden end if - end subroutine end_log - + if (present(iter)) iter = it + end subroutine log_end - subroutine psb_d_init_conv(methdname,stopc,trace,a,b,eps,desc_a,stopdat,info) + subroutine psb_d_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) use psb_base_mod implicit none character(len=*), intent(in) :: methdname - integer, intent(in) :: stopc, trace + integer, intent(in) :: stopc, trace,itmax type(psb_dspmat_type), intent(in) :: a real(kind(1.d0)), intent(in) :: b(:), eps type(psb_desc_type), intent(in) :: desc_a @@ -500,10 +536,6 @@ contains integer, intent(out) :: info integer :: ictxt, me, np, err_act - character(len=*), parameter :: fmt='(a18,1x,A4,3(2x,A10))' - integer, parameter :: outlen=18 - character(len=len(methdname)) :: mname - character(len=outlen) :: outname character(len=20) :: name info = 0 @@ -520,6 +552,7 @@ contains stopdat%controls(stopc_) = stopc stopdat%controls(trace_) = trace + stopdat%controls(itmax_) = itmax select case(stopdat%controls(stopc_)) case (1) @@ -543,13 +576,9 @@ contains stopdat%values(errnum_) = dzero stopdat%values(errden_) = done - if (stopdat%controls(trace_) > 0) then - if (me == 0) then - mname = adjustl(trim(methdname)) - write(outname,'(a)') mname(1:min(len_trim(mname),outlen-1))//':' - write(*,fmt) adjustl(outname),'Iter','Conv. Ind.','Epsilon' - endif - end if + if ((stopdat%controls(trace_) > 0).and. (me == 0))& + & call log_header(methdname) + call psb_erractionrestore(err_act) return @@ -562,11 +591,11 @@ contains end subroutine psb_d_init_conv - subroutine psb_z_init_conv(methdname,stopc,trace,a,b,eps,desc_a,stopdat,info) + subroutine psb_z_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) use psb_base_mod implicit none character(len=*), intent(in) :: methdname - integer, intent(in) :: stopc, trace + integer, intent(in) :: stopc, trace, itmax type(psb_zspmat_type), intent(in) :: a complex(kind(1.d0)), intent(in) :: b(:) real(kind(1.d0)), intent(in) :: eps @@ -575,10 +604,6 @@ contains integer, intent(out) :: info integer :: ictxt, me, np, err_act - character(len=*), parameter :: fmt='(a18,1x,A4,3(2x,A10))' - integer, parameter :: outlen=18 - character(len=len(methdname)) :: mname - character(len=outlen) :: outname character(len=20) :: name info = 0 @@ -595,6 +620,7 @@ contains stopdat%controls(stopc_) = stopc stopdat%controls(trace_) = trace + stopdat%controls(itmax_) = itmax select case(stopdat%controls(stopc_)) case (1) @@ -618,14 +644,9 @@ contains stopdat%values(errnum_) = dzero stopdat%values(errden_) = done - - if (stopdat%controls(trace_) > 0) then - if (me == 0) then - mname = adjustl(trim(methdname)) - write(outname,'(a)') mname(1:min(len_trim(mname),outlen-1))//':' - write(*,fmt) adjustl(outname),'Iter','Conv. Ind.','Epsilon' - endif - end if + if ((stopdat%controls(trace_) > 0).and. (me == 0))& + & call log_header(methdname) + call psb_erractionrestore(err_act) return @@ -651,11 +672,6 @@ contains integer, intent(out) :: info integer :: ictxt, me, np, err_act - real(kind(1.d0)) :: errnum, errden - character(len=*), parameter :: fmt='(a18,1x,i4,3(2x,es10.4))' - integer, parameter :: outlen=18 - character(len=len(methdname)) :: mname - character(len=outlen) :: outname character(len=20) :: name info = 0 @@ -664,6 +680,7 @@ contains ictxt = psb_cd_get_context(desc_a) call psb_info(ictxt,me,np) + psb_d_check_conv = .false. select case(stopdat%controls(stopc_)) @@ -674,7 +691,7 @@ contains stopdat%values(errden_) = & & (stopdat%values(ani_)*stopdat%values(xni_)+stopdat%values(bni_)) case(2) - stopdat%values(rn2_) = psb_genrm2(r,desc_a,info) + stopdat%values(rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(errnum_) = stopdat%values(rn2_) stopdat%values(errden_) = stopdat%values(bn2_) @@ -695,21 +712,13 @@ contains psb_d_check_conv = & & (stopdat%values(errnum_) <= stopdat%values(eps_)*stopdat%values(errden_)) end if + + psb_d_check_conv = (psb_d_check_conv.or.(stopdat%controls(itmax_) <= it)) if (((stopdat%controls(trace_) > 0).and.(mod(it,stopdat%controls(trace_))==0))& & .or.psb_d_check_conv) then - if (me == 0) then - mname = adjustl(trim(methdname)) - write(outname,'(a)') mname(1:min(len_trim(mname),outlen-1))//':' - if (stopdat%values(errden_) > dzero ) then - write(*,fmt) adjustl(outname),it,& - & stopdat%values(errnum_)/stopdat%values(errden_),stopdat%values(eps_) - else - write(*,fmt) adjustl(outname),it,& - & stopdat%values(errnum_),stopdat%values(eps_) - end if - Endif - + call log_conv(methdname,me,it,1,stopdat%values(errnum_),& + & stopdat%values(errden_),stopdat%values(eps_)) end if call psb_erractionrestore(err_act) @@ -737,11 +746,6 @@ contains integer, intent(out) :: info integer :: ictxt, me, np, err_act - real(kind(1.d0)) :: errnum, errden - character(len=*), parameter :: fmt='(a18,1x,i4,3(2x,es10.4))' - integer, parameter :: outlen=18 - character(len=len(methdname)) :: mname - character(len=outlen) :: outname character(len=20) :: name info = 0 @@ -781,21 +785,13 @@ contains psb_z_check_conv = & & (stopdat%values(errnum_) <= stopdat%values(eps_)*stopdat%values(errden_)) end if + + psb_z_check_conv = (psb_z_check_conv.or.(stopdat%controls(itmax_) <= it)) if (((stopdat%controls(trace_) > 0).and.(mod(it,stopdat%controls(trace_))==0))& & .or.psb_z_check_conv) then - if (me == 0) then - mname = adjustl(trim(methdname)) - write(outname,'(a)') mname(1:min(len_trim(mname),outlen-1))//':' - if (stopdat%values(errden_) > dzero ) then - write(*,fmt) adjustl(outname),it,& - & stopdat%values(errnum_)/stopdat%values(errden_),stopdat%values(eps_) - else - write(*,fmt) adjustl(outname),it,& - & stopdat%values(errnum_),stopdat%values(eps_) - end if - Endif - + call log_conv(methdname,me,it,1,stopdat%values(errnum_),& + & stopdat%values(errden_),stopdat%values(eps_)) end if call psb_erractionrestore(err_act) @@ -829,52 +825,18 @@ contains info = 0 name = 'psb_end_conv' - call psb_erractionsave(err_act) ictxt = psb_cd_get_context(desc_a) call psb_info(ictxt,me,np) - if (present(iter)) iter = it errnum = stopdat%values(errnum_) errden = stopdat%values(errden_) eps = stopdat%values(eps_) - - - if (errden == dzero) then - if (errnum > eps) then - - if (me==0) then - write(*,fmt) trim(methdname)//' failed to converge to ',eps,& - & ' in ',it,' iterations. ' - write(*,fmt1) 'Last iteration convergence indicator: ',& - & errnum - end if - end if - if (present(err)) err=errnum - else - if (errnum > eps) then - - if (me==0) then - write(*,fmt) trim(methdname)//' failed to converge to ',eps,& - & ' in ',it,' iterations. ' - write(*,fmt1) 'Last iteration convergence indicator: ',& - & errnum/errden - end if - endif - if (present(err)) err=errnum/errden - end if - call psb_erractionrestore(err_act) - return + call log_end(methdname,me,it,errnum,errden,eps,err,iter) -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - + end subroutine psb_end_conv diff --git a/krylov/psb_zcgs.f90 b/krylov/psb_zcgs.f90 index e46a97a2..8c5b4850 100644 --- a/krylov/psb_zcgs.f90 +++ b/krylov/psb_zcgs.f90 @@ -191,7 +191,7 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) itx = 0 - call psb_init_conv(methdname,istop_,itrace_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) if (info /= 0) Then call psb_errpush(4011,name) goto 9999 diff --git a/krylov/psb_zcgstab.f90 b/krylov/psb_zcgstab.f90 index b3fe6f7c..098ef7d3 100644 --- a/krylov/psb_zcgstab.f90 +++ b/krylov/psb_zcgstab.f90 @@ -112,7 +112,7 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) Complex(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:) Complex(Kind(1.d0)), Pointer :: q(:),& & r(:), p(:), v(:), s(:), t(:), z(:), f(:) - Integer :: litmax, naux, mglob, it,itrace_,& + Integer :: itmax_, naux, mglob, it,itrace_,& & np,me, n_row, n_col integer :: debug_level, debug_unit Integer :: itx, isvch, ictxt, err_act @@ -177,9 +177,9 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) Z => WWRK(:,8) if (present(itmax)) then - litmax = itmax + itmax_ = itmax else - litmax = 1000 + itmax_ = 1000 endif if (present(itrace)) then @@ -194,7 +194,7 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) itx = 0 - call psb_init_conv(methdname,istop_,itrace_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) if (info /= 0) Then call psb_errpush(4011,name) goto 9999 @@ -205,7 +205,7 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) !!$ !!$ r0 = b-Ax0 !!$ - if (itx >= litmax) exit restart + if (itx >= itmax_) exit restart it = 0 call psb_geaxpby(zone,b,zzero,r,desc_a,info) if (info == 0) call psb_spmm(-zone,a,x,zone,r,desc_a,info,work=aux) diff --git a/krylov/psb_zrgmres.f90 b/krylov/psb_zrgmres.f90 index 775e03ba..855a5349 100644 --- a/krylov/psb_zrgmres.f90 +++ b/krylov/psb_zrgmres.f90 @@ -254,10 +254,11 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist call psb_errpush(info,name) goto 9999 end if + if ((itrace_ > 0).and.(me==0)) call log_header(methdname) itx = 0 - restart: Do - + restart: do + ! compute r0 = b-ax0 ! check convergence ! compute v1 = r0/||r0||_2 @@ -435,16 +436,7 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist if (itrace_ > 0) & & call log_conv(methdname,me,itx,1,errnum,errden,eps) - if (present(err)) then - if (errden /= dzero) then - err = errnum/errden - else - err = errnum - end if - end if - if (present(iter)) iter = itx - if (errnum > eps*errden) & - & call end_log(methdname,me,itx,errnum,errden,eps) + call log_end(methdname,me,itx,errnum,errden,eps,err=err,iter=iter) deallocate(aux,h,c,s,rs,rst, stat=info) if (info == 0) call psb_gefree(v,desc_a,info)