Krylov iterations:

Centralized check convergence routines, with some fixes for RGMRES.
17 years ago · 9a58493988
parent 649e356750
commit 9a58493988
10 changed files with 104 additions and 160 deletions
--- 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
--- 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
--- 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
--- 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
--- 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)
--- 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)
--- 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'
+    
+  end subroutine log_header

-    if (me == 0) then 
+
+  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).and. (me == 0))&
+         &  call log_header(methdname) 

-    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
    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_)

@ -696,20 +713,12 @@ contains
           & (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
@ -782,20 +786,12 @@ contains
           & (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,51 +825,17 @@ 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_) 

+    call log_end(methdname,me,it,errnum,errden,eps,err,iter)

-    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
-
-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

--- 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
--- 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)
--- a/krylov/psb_zrgmres.f90
+++ b/krylov/psb_zrgmres.f90
@ -254,9 +254,10 @@ 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
@ -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)