From bc980ed5d67f45bc528e289fae49d19414fdc482 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Thu, 10 Jan 2008 17:10:14 +0000 Subject: [PATCH] Modified application of stopping criterion to be errnum <= eps*errden and relevant informational messages. --- krylov/psb_dbicg.f90 | 69 ++++++++++++++++++++++++----------------- krylov/psb_dcg.f90 | 51 ++++++++++++++++++++---------- krylov/psb_dcgs.f90 | 52 ++++++++++++++++++++----------- krylov/psb_dcgstab.F90 | 59 +++++++++++++++++++++++------------ krylov/psb_dcgstabl.f90 | 56 +++++++++++++++++++++------------ krylov/psb_drgmres.f90 | 60 ++++++++++++++++++++++------------- krylov/psb_zcgs.f90 | 52 ++++++++++++++++++++----------- krylov/psb_zcgstab.f90 | 53 ++++++++++++++++++++----------- krylov/psb_zrgmres.f90 | 56 +++++++++++++++++++++------------ 9 files changed, 328 insertions(+), 180 deletions(-) diff --git a/krylov/psb_dbicg.f90 b/krylov/psb_dbicg.f90 index 0f295c29..408e3ff2 100644 --- a/krylov/psb_dbicg.f90 +++ b/krylov/psb_dbicg.f90 @@ -76,16 +76,21 @@ ! performed. ! iter - integer(optional) Output: how many iterations have been ! performed. -! err - real (optional) Output: error estimate on exit -! itrace - integer(optional) Input: print an informational message -! with the error estimate every itrace -! iterations -! istop - integer(optional) Input: stopping criterion, or how -! to estimate the error. -! 1: err = |r|/|b| -! 2: err = |r|/(|a||x|+|b|) -! where r is the (preconditioned, recursive -! estimate of) residual +! performed. +! err - real (optional) Output: error estimate on exit. If the +! denominator of the estimate is exactly +! 0, it is changed into 1. +! itrace - integer(optional) Input: print an informational message +! with the error estimate every itrace +! iterations +! istop - integer(optional) Input: stopping criterion, or how +! to estimate the error. +! 1: err = |r|/|b|; here the iteration is +! stopped when |r| <= eps * |b| +! 2: err = |r|/(|a||x|+|b|); here the iteration is +! stopped when |r| <= eps * (|a||x|+|b|) +! where r is the (preconditioned, recursive +! estimate of) residual. ! ! subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) @@ -109,7 +114,6 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) real(kind(1.d0)), pointer :: ww(:), q(:),& & r(:), p(:), zt(:), pt(:), z(:), rt(:),qt(:) integer :: int_err(5) - real(kind(1.d0)) ::rerr integer ::litmax, naux, mglob, it, itrace_,& & np,me, n_row, n_col, istop_, err_act integer :: debug_level, debug_unit @@ -118,6 +122,7 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) integer :: itx, isvch, ictxt real(kind(1.d0)) :: alpha, beta, rho, rho_old, rni, xni, bni, ani,& & sigma,bn2 + real(kind(1.d0)) :: errnum, errden character(len=20) :: name,ch_err info = 0 @@ -214,6 +219,8 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) else if (istop_ == 2) then bn2 = psb_genrm2(b,desc_a,info) endif + errnum = dzero + errden = done if(info.ne.0) then info=4011 @@ -256,9 +263,11 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) if (istop_ == 1) then xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) else if (istop_ == 2) then - rerr = rni/bn2 + errnum = rni + errden = bn2 endif if(info.ne.0) then @@ -267,12 +276,12 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) goto 9999 end if - if (rerr<=eps) then + if (errnum <= eps*errden) Then exit restart end if If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'bicg: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'bicg: ',itx,errnum,eps*errden end If @@ -327,17 +336,14 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) if (istop_ == 1) then rni = psb_geamax(r,desc_a,info) xni = psb_geamax(x,desc_a,info) + errnum = rni + errden = (ani*xni+bni) else if (istop_ == 2) then rni = psb_genrm2(r,desc_a,info) + errnum = rni + errden = bn2 endif - - if (istop_ == 1) then - xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) - else if (istop_ == 2) then - rerr = rni/bn2 - endif - if (rerr<=eps) then + If (errnum <= eps*errden) Then exit restart end if @@ -345,18 +351,25 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'bicg: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'bicg: ',itx,errnum,eps*errden end If end do iteration end do restart If (itrace_ > 0) then - if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicg: ',itx,rerr + if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicg: ',itx,errnum,eps*errden end If - if (present(err)) err=rerr + if (present(err)) then + if (errden /= dzero) then + err = errnum/errden + else + err = errnum + end if + end if + if (present(iter)) iter = itx - if (rerr>eps) then - write(debug_unit,*) 'bicg failed to converge to ',eps,& + If ((errnum > eps*errden).and.(me==0)) Then + write(debug_unit,*) 'bicg failed to converge to ',eps*errden,& & ' in ',itx,' iterations ' end if diff --git a/krylov/psb_dcg.f90 b/krylov/psb_dcg.f90 index d2b300df..2ac03508 100644 --- a/krylov/psb_dcg.f90 +++ b/krylov/psb_dcg.f90 @@ -77,16 +77,21 @@ ! performed. ! iter - integer(optional) Output: how many iterations have been ! performed. -! err - real (optional) Output: error estimate on exit +! performed. +! err - real (optional) Output: error estimate on exit. If the +! denominator of the estimate is exactly +! 0, it is changed into 1. ! itrace - integer(optional) Input: print an informational message ! with the error estimate every itrace ! iterations ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/|b| -! 2: err = |r|/(|a||x|+|b|) +! 1: err = |r|/|b|; here the iteration is +! stopped when |r| <= eps * |b| +! 2: err = |r|/(|a||x|+|b|); here the iteration is +! stopped when |r| <= eps * (|a||x|+|b|) ! where r is the (preconditioned, recursive -! estimate of) residual +! estimate of) residual. ! ! Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) @@ -108,12 +113,12 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) !!$ Local data real(kind(1.d0)), allocatable, target :: aux(:), wwrk(:,:) real(kind(1.d0)), pointer :: q(:), p(:), r(:), z(:), w(:) - real(kind(1.d0)) ::rerr real(kind(1.d0)) ::alpha, beta, rho, rho_old, rni, xni, bni, ani,bn2,& & sigma integer :: litmax, istop_, naux, mglob, it, itx, itrace_,& & np,me, n_col, isvch, ictxt, n_row,err_act, int_err(5) logical, parameter :: exchange=.true., noexchange=.false. + real(kind(1.d0)) :: errnum, errden character(len=20) :: name info = 0 @@ -216,6 +221,8 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) else if (istop_ == 2) then bn2 = psb_genrm2(b,desc_a,info) endif + errnum = dzero + errden = done if (info.ne.0) then info=4011 call psb_errpush(info,name) @@ -254,33 +261,43 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) call psb_geaxpby(-alpha,q,done,r,desc_a,info) - if (istop_ == 1) Then + if (istop_ == 1) then rni = psb_geamax(r,desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) - Else If (istop_ == 2) Then - + errnum = rni + errden = (ani*xni+bni) + else if (istop_ == 2) then rni = psb_genrm2(r,desc_a,info) - rerr = rni/bn2 - Endif - if (rerr<=eps) exit restart + errnum = rni + errden = bn2 + endif + if (errnum <= eps*errden) then + exit restart + end if if (itx>= litmax) exit restart If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'cg: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'cg: ',itx,errnum,eps*errden end If end do iteration end do restart If (itrace_ > 0) then - if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'cg: ',itx,rerr + if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'cg: ',itx,errnum,eps*errden + end If + + If (Present(err)) then + if (errden /= dzero) then + err = errnum/errden + else + err = errnum + end if end If - if (present(err)) err=rerr if (present(iter)) iter = itx - if (rerr>eps) then - write(0,*) 'CG Failed to converge to ',eps,& + If ((errnum > eps*errden).and.(me==0)) Then + write(0,*) 'CG Failed to converge to ',eps*errden,& & ' in ',litmax,' iterations ' info=itx end if diff --git a/krylov/psb_dcgs.f90 b/krylov/psb_dcgs.f90 index 85b4aa18..20712484 100644 --- a/krylov/psb_dcgs.f90 +++ b/krylov/psb_dcgs.f90 @@ -76,17 +76,21 @@ ! performed. ! iter - integer(optional) Output: how many iterations have been ! performed. -! err - real (optional) Output: error estimate on exit +! performed. +! err - real (optional) Output: error estimate on exit. If the +! denominator of the estimate is exactly +! 0, it is changed into 1. ! itrace - integer(optional) Input: print an informational message ! with the error estimate every itrace ! iterations ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/|b| -! 2: err = |r|/(|a||x|+|b|) +! 1: err = |r|/|b|; here the iteration is +! stopped when |r| <= eps * |b| +! 2: err = |r|/(|a||x|+|b|); here the iteration is +! stopped when |r| <= eps * (|a||x|+|b|) ! where r is the (preconditioned, recursive -! estimate of) residual -! +! estimate of) residual. ! ! Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) @@ -109,7 +113,6 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) Real(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:) Real(Kind(1.d0)), Pointer :: ww(:), q(:),& & r(:), p(:), v(:), s(:), z(:), f(:), rt(:),qt(:),uv(:) - Real(Kind(1.d0)) :: rerr Integer :: litmax, naux, mglob, it, itrace_,int_err(5),& & np,me, n_row, n_col,istop_, err_act Logical, Parameter :: exchange=.True., noexchange=.False. @@ -118,6 +121,7 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) integer :: debug_level, debug_unit Real(Kind(1.d0)) :: alpha, beta, rho, rho_old, rni, xni, bni, ani,bn2,& & sigma + real(kind(1.d0)) :: errnum, errden character(len=20) :: name info = 0 @@ -212,6 +216,8 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) else if (istop_ == 2) then bn2 = psb_genrm2(b,desc_a,info) endif + errnum = dzero + errden = done if(info/=0)then info=4011 call psb_errpush(info,name) @@ -241,22 +247,24 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) if (istop_ == 1) then rni = psb_geamax(r,desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) else if (istop_ == 2) then rni = psb_genrm2(r,desc_a,info) - rerr = rni/bn2 + errnum = rni + errden = bn2 endif if(info/=0)then info=4011 call psb_errpush(info,name) goto 9999 end if - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden end If iteration: Do @@ -318,13 +326,15 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) if (istop_ == 1) then rni = psb_geamax(r,desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) else if (istop_ == 2) then rni = psb_genrm2(r,desc_a,info) - rerr = rni/bn2 + errnum = rni + errden = bn2 endif - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If @@ -332,7 +342,7 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden end If End Do iteration @@ -340,14 +350,20 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden end If - If (Present(err)) err=rerr If (Present(iter)) iter = itx - If (rerr>eps) Then - write(debug_unit,*) 'cgs failed to converge to ',eps,& + If (Present(err)) then + if (errden /= dzero) then + err = errnum/errden + else + err = errnum + end if + end If + If ((errnum > eps*errden).and.(me==0)) Then + write(debug_unit,*) 'cgs failed to converge to ',eps*errden,& & ' in ',itx,' iterations ' End If diff --git a/krylov/psb_dcgstab.F90 b/krylov/psb_dcgstab.F90 index 00b98d51..6f6dda78 100644 --- a/krylov/psb_dcgstab.F90 +++ b/krylov/psb_dcgstab.F90 @@ -76,16 +76,20 @@ ! performed. ! iter - integer(optional) Output: how many iterations have been ! performed. -! err - real (optional) Output: error estimate on exit +! err - real (optional) Output: error estimate on exit. If the +! denominator of the estimate is exactly +! 0, it is changed into 1. ! itrace - integer(optional) Input: print an informational message ! with the error estimate every itrace ! iterations ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/|b| -! 2: err = |r|/(|a||x|+|b|) +! 1: err = |r|/|b|; here the iteration is +! stopped when |r| <= eps * |b| +! 2: err = |r|/(|a||x|+|b|); here the iteration is +! stopped when |r| <= eps * (|a||x|+|b|) ! where r is the (preconditioned, recursive -! estimate of) residual +! estimate of) residual. ! ! ! @@ -108,7 +112,6 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) Real(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:) Real(Kind(1.d0)), Pointer :: q(:),& & r(:), p(:), v(:), s(:), t(:), z(:), f(:) - Real(Kind(1.d0)) :: rerr Integer :: litmax, naux, mglob, it,itrace_,& & np,me, n_row, n_col integer :: debug_level, debug_unit @@ -116,8 +119,9 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) Integer, Parameter :: irmax = 8 Integer :: itx, isvch, ictxt, err_act, int_err(5) Integer :: istop_ - Real(Kind(1.d0)) :: alpha, beta, rho, rho_old, rni, xni, bni, ani,& + Real(Kind(1.d0)) :: alpha, beta, rho, rho_old, rni, xni, bni, ani,& & sigma, omega, tau, rn0, bn2 + real(kind(1.d0)) :: errnum, errden #ifdef MPE_KRYLOV Integer istpb, istpe, ifctb, ifcte, imerr, irank, icomm,immb,imme Integer mpe_log_get_event_number,mpe_Describe_state,mpe_log_event @@ -231,6 +235,8 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) Else If (istop_ == 2) Then bn2 = psb_genrm2(b,desc_a,info) Endif + errnum = dzero + errden = done if (info /= 0) Then info=4011 call psb_errpush(info,name) @@ -282,18 +288,22 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) If (itx == 0) Then rn0 = rni End If - If (rn0 == 0.d0 ) Then + If (rn0 == dzero ) Then If (itrace_ > 0 ) Then If (me == 0) Write(*,*) 'BiCGSTAB: ',itx,rn0 Endif + errnum = dzero + errden = done Exit restart End If If (istop_ == 1) Then xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) Else If (istop_ == 2) Then - rerr = rni/bn2 + errnum = rni + errden = bn2 Endif if (info /= 0) Then info=4011 @@ -301,13 +311,14 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) goto 9999 End If - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If If (itrace_ > 0) then if (((itx==0).or.(mod(itx,itrace_)==0)).and.(me == 0)) & - & write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',& + & itx,errnum,eps*errden Endif iteration: Do @@ -428,13 +439,14 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) If (istop_ == 1) Then rni = psb_geamax(r,desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) - + errnum = rni + errden = (ani*xni+bni) Else If (istop_ == 2) Then rni = psb_genrm2(r,desc_a,info) - rerr = rni/bn2 + errnum = rni + errden = bn2 Endif - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If @@ -443,20 +455,27 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0)) & & write(*,'(a,i4,3(2x,es10.4))') & - & 'bicgstab: ',itx,rerr + & 'bicgstab: ',itx,errnum,eps*errden Endif End Do iteration End Do restart If (itrace_ > 0) then - if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',itx,rerr + if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',& + & itx,errnum,eps*errden Endif - If (Present(err)) err=rerr + If (Present(err)) then + if (errden /= dzero) then + err = errnum/errden + else + err = errnum + end if + end If If (Present(iter)) iter = itx - If (rerr>eps) Then - write(debug_unit,*) 'BI-CGSTAB failed to converge to ',EPS,& + If ((errnum > eps*errden).and.(me==0)) Then + write(debug_unit,*) 'BI-CGSTAB failed to converge to ',eps*errden,& & ' in ',ITX,' iterations. ' End If diff --git a/krylov/psb_dcgstabl.f90 b/krylov/psb_dcgstabl.f90 index 478ff446..c37da9fc 100644 --- a/krylov/psb_dcgstabl.f90 +++ b/krylov/psb_dcgstabl.f90 @@ -83,19 +83,23 @@ ! performed. ! iter - integer(optional) Output: how many iterations have been ! performed. -! err - real (optional) Output: error estimate on exit +! performed. +! err - real (optional) Output: error estimate on exit. If the +! denominator of the estimate is exactly +! 0, it is changed into 1. ! itrace - integer(optional) Input: print an informational message ! with the error estimate every itrace ! iterations -! irst - integer(optional) Input: restart parameter L -! ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/|b| -! 2: err = |r|/(|a||x|+|b|) +! 1: err = |r|/|b|; here the iteration is +! stopped when |r| <= eps * |b| +! 2: err = |r|/(|a||x|+|b|); here the iteration is +! stopped when |r| <= eps * (|a||x|+|b|) ! where r is the (preconditioned, recursive -! estimate of) residual -! +! estimate of) residual. +! irst - integer(optional) Input: restart parameter L +! ! ! Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) @@ -119,7 +123,6 @@ 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(:) - Real(Kind(1.d0)) :: rerr Integer :: litmax, naux, mglob, it, itrace_,& & np,me, n_row, n_col, nl, err_act Logical, Parameter :: exchange=.True., noexchange=.False. @@ -128,6 +131,7 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is integer :: debug_level, debug_unit Real(Kind(1.d0)) :: alpha, beta, rho, rho_old, rni, xni, bni, ani,bn2,& & omega + real(kind(1.d0)) :: errnum, errden character(len=20) :: name info = 0 @@ -249,6 +253,8 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is else if (istop_ == 2) then bn2 = psb_genrm2(b,desc_a,info) endif + errnum = dzero + errden = done if (info.ne.0) Then info=4011 call psb_errpush(info,name) @@ -289,10 +295,12 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is if (istop_ == 1) then rni = psb_geamax(r,desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) else if (istop_ == 2) then rni = psb_genrm2(r,desc_a,info) - rerr = rni/bn2 + errnum = rni + errden = bn2 endif if (info.ne.0) Then info=4011 @@ -300,12 +308,12 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is goto 9999 End If - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'bicgstab(l): ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'bicgstab(l): ',itx,errnum,eps*errden end If iteration: Do @@ -414,32 +422,40 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is if (istop_ == 1) then rni = psb_geamax(rh(:,0),desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) else if (istop_ == 2) then rni = psb_genrm2(rh(:,0),desc_a,info) - rerr = rni/bn2 + errnum = rni + errden = bn2 endif - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If If (itx.Ge.litmax) Exit restart If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'bicgstab(l): ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'bicgstab(l): ',itx,errnum,eps*errden end If End Do iteration End Do restart If (itrace_ > 0) then - if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicgstab(l): ',itx,rerr + if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicgstab(l): ',itx,errnum,eps*errden + end If + If (Present(err)) then + if (errden /= dzero) then + err = errnum/errden + else + err = errnum + end if end If - If (Present(err)) err=rerr If (Present(iter)) iter = itx - If (rerr>eps) Then - write(debug_unit,*) 'bi-cgstabl failed to converge to ',eps,& + If ((errnum > eps*errden).and.(me==0)) Then + write(debug_unit,*) 'bi-cgstabl failed to converge to ',eps*errden,& & ' in ',itx,' iterations ' End If diff --git a/krylov/psb_drgmres.f90 b/krylov/psb_drgmres.f90 index 5e0a4f48..ec033e1e 100644 --- a/krylov/psb_drgmres.f90 +++ b/krylov/psb_drgmres.f90 @@ -87,18 +87,23 @@ ! performed. ! iter - integer(optional) Output: how many iterations have been ! performed. -! err - real (optional) Output: error estimate on exit +! performed. +! err - real (optional) Output: error estimate on exit. If the +! denominator of the estimate is exactly +! 0, it is changed into 1. ! itrace - integer(optional) Input: print an informational message ! with the error estimate every itrace ! iterations -! irst - integer(optional) Input: restart parameter -! ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/|b| -! 2: err = |r|/(|a||x|+|b|) +! 1: err = |r|/|b|; here the iteration is +! stopped when |r| <= eps * |b| +! 2: err = |r|/(|a||x|+|b|); here the iteration is +! stopped when |r| <= eps * (|a||x|+|b|) ! where r is the (preconditioned, recursive -! estimate of) residual +! estimate of) residual. +! irst - integer(optional) Input: restart parameter +! ! Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_base_mod @@ -119,7 +124,7 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist !!$ local data Real(Kind(1.d0)), allocatable, target :: aux(:),w(:),w1(:), v(:,:) Real(Kind(1.d0)), allocatable :: c(:),s(:), h(:,:), rs(:),rst(:),xt(:) - Real(Kind(1.d0)) :: rerr, scal, gm, rti, rti1 + Real(Kind(1.d0)) :: scal, gm, rti, rti1 Integer ::litmax, naux, mglob, it,k, itrace_,& & np,me, n_row, n_col, nl, int_err(5) Logical, Parameter :: exchange=.True., noexchange=.False., use_drot=.true. @@ -128,6 +133,7 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist integer :: debug_level, debug_unit Real(Kind(1.d0)) :: rni, xni, bni, ani,bn2, dt real(kind(1.d0)), external :: dnrm2 + real(kind(1.d0)) :: errnum, errden character(len=20) :: name info = 0 @@ -239,6 +245,8 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist else if (istop_ == 2) then bn2 = psb_genrm2(b,desc_a,info) endif + errnum = dzero + errden = done if (info.ne.0) Then info=4011 call psb_errpush(info,name) @@ -289,10 +297,12 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist if (istop_ == 1) then rni = psb_geamax(v(:,1),desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) else if (istop_ == 2) then rni = psb_genrm2(v(:,1),desc_a,info) - rerr = rni/bn2 + errnum = rni + errden = bn2 endif if (info.ne.0) Then info=4011 @@ -300,13 +310,13 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist goto 9999 End If - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden end If v(:,1) = v(:,1) * scal @@ -379,7 +389,8 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist call psb_spmm(-done,a,xt,done,w1,desc_a,info,work=aux) rni = psb_geamax(w1,desc_a,info) xni = psb_geamax(xt,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) ! else if (istop_ == 2) then @@ -388,12 +399,12 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist ! procedure of the least-squares problem ! rni = abs(rs(i+1)) - rerr = rni/bn2 - + errnum = rni + errden = bn2 endif - if (rerr < eps ) then - + If (errnum <= eps*errden) Then + if (istop_ == 1) then x = xt else if (istop_ == 2) then @@ -418,7 +429,7 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden end If end Do inner @@ -443,13 +454,20 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist End Do restart If (itrace_ > 0) then - if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,rerr + if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden + end If + + If (Present(err)) then + if (errden /= dzero) then + err = errnum/errden + else + err = errnum + end if end If - If (Present(err)) err=rerr If (Present(iter)) iter = itx - If ((rerr>eps).and. (me == 0)) Then - write(debug_unit,*) 'gmresr(l) failed to converge to ',eps,& + If ((errnum > eps*errden).and.(me==0)) Then + write(debug_unit,*) 'gmresr(l) failed to converge to ',eps*errden,& & ' in ',itx,' iterations ' End If diff --git a/krylov/psb_zcgs.f90 b/krylov/psb_zcgs.f90 index be036dc9..615904bb 100644 --- a/krylov/psb_zcgs.f90 +++ b/krylov/psb_zcgs.f90 @@ -75,16 +75,21 @@ ! performed. ! iter - integer(optional) Output: how many iterations have been ! performed. -! err - real (optional) Output: error estimate on exit +! performed. +! err - real (optional) Output: error estimate on exit. If the +! denominator of the estimate is exactly +! 0, it is changed into 1. ! itrace - integer(optional) Input: print an informational message ! with the error estimate every itrace ! iterations ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/|b| -! 2: err = |r|/(|a||x|+|b|) +! 1: err = |r|/|b|; here the iteration is +! stopped when |r| <= eps * |b| +! 2: err = |r|/(|a||x|+|b|); here the iteration is +! stopped when |r| <= eps * (|a||x|+|b|) ! where r is the (preconditioned, recursive -! estimate of) residual +! estimate of) residual. ! Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_base_mod @@ -106,7 +111,6 @@ Subroutine psb_zcgs(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 :: ww(:), q(:),& & r(:), p(:), v(:), s(:), z(:), f(:), rt(:),qt(:),uv(:) - Real(Kind(1.d0)) :: rerr Integer :: litmax, naux, mglob, it, itrace_,int_err(5),& & np,me, n_row, n_col,istop_, err_act Logical, Parameter :: exchange=.True., noexchange=.False. @@ -115,6 +119,7 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) integer :: debug_level, debug_unit Real(Kind(1.d0)) :: rni, xni, bni, ani,bn2 complex(Kind(1.d0)) :: alpha, beta, rho, rho_old, sigma + real(kind(1.d0)) :: errnum, errden character(len=20) :: name info = 0 @@ -209,6 +214,8 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) else if (istop_ == 2) then bn2 = psb_genrm2(b,desc_a,info) endif + errnum = dzero + errden = done if(info/=0)then info=4011 call psb_errpush(info,name) @@ -238,10 +245,12 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) if (istop_ == 1) then rni = psb_geamax(r,desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) else if (istop_ == 2) then rni = psb_genrm2(r,desc_a,info) - rerr = rni/bn2 + errnum = rni + errden = bn2 endif if(info/=0)then info=4011 @@ -249,12 +258,12 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) goto 9999 end if - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden end If iteration: Do @@ -316,13 +325,15 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) if (istop_ == 1) then rni = psb_geamax(r,desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) else if (istop_ == 2) then rni = psb_genrm2(r,desc_a,info) - rerr = rni/bn2 + errnum = rni + errden = bn2 endif - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If @@ -330,21 +341,26 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden end If End Do iteration End Do restart If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden end If - - If (Present(err)) err=rerr + If (Present(err)) then + if (errden /= dzero) then + err = errnum/errden + else + err = errnum + end if + end If If (Present(iter)) iter = itx - If (rerr>eps) Then - write(debug_unit,*) 'cgs failed to converge to ',eps,& + If ((errnum > eps*errden).and.(me==0)) Then + write(debug_unit,*) 'cgs failed to converge to ',eps*errden,& & ' in ',itx,' iterations ' End If diff --git a/krylov/psb_zcgstab.f90 b/krylov/psb_zcgstab.f90 index 7dda3936..1a4515d2 100644 --- a/krylov/psb_zcgstab.f90 +++ b/krylov/psb_zcgstab.f90 @@ -76,16 +76,21 @@ ! performed. ! iter - integer(optional) Output: how many iterations have been ! performed. -! err - real (optional) Output: error estimate on exit +! performed. +! err - real (optional) Output: error estimate on exit. If the +! denominator of the estimate is exactly +! 0, it is changed into 1. ! itrace - integer(optional) Input: print an informational message ! with the error estimate every itrace ! iterations ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/|b| -! 2: err = |r|/(|a||x|+|b|) +! 1: err = |r|/|b|; here the iteration is +! stopped when |r| <= eps * |b| +! 2: err = |r|/(|a||x|+|b|); here the iteration is +! stopped when |r| <= eps * (|a||x|+|b|) ! where r is the (preconditioned, recursive -! estimate of) residual +! estimate of) residual. ! Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_base_mod @@ -106,7 +111,6 @@ 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(:) - Real(Kind(1.d0)) :: rerr Integer :: litmax, naux, mglob, it,itrace_,& & np,me, n_row, n_col integer :: debug_level, debug_unit @@ -115,7 +119,8 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) Integer :: itx, isvch, ictxt, err_act, int_err(5) Integer :: istop_ complex(Kind(1.d0)) :: alpha, beta, rho, rho_old, sigma, omega, tau - Real(Kind(1.d0)) :: rni, xni, bni, ani, rn0, bn2 + Real(Kind(1.d0)) :: rni, xni, bni, ani, rn0, bn2 + real(kind(1.d0)) :: errnum, errden !!$ Integer istpb, istpe, ifctb, ifcte, imerr, irank, icomm,immb,imme !!$ Integer mpe_log_get_event_number,mpe_Describe_state,mpe_log_event character(len=20) :: name @@ -243,6 +248,8 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) Else If (istop_ == 2) Then rni = psb_genrm2(r,desc_a,info) Endif + errnum = dzero + errden = done if (info /= 0) Then info=4011 call psb_errpush(info,name) @@ -261,9 +268,11 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) If (istop_ == 1) Then xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) Else If (istop_ == 2) Then - rerr = rni/bn2 + errnum = rni + errden = bn2 Endif if (info /= 0) Then info=4011 @@ -271,13 +280,13 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) goto 9999 End If - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If If (itrace_ > 0) then if (((itx==0).or.(mod(itx,itrace_)==0)).and.(me == 0)) & - & write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',itx,errnum,eps*errden end If iteration: Do @@ -374,13 +383,15 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) If (istop_ == 1) Then rni = psb_geamax(r,desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) Else If (istop_ == 2) Then rni = psb_genrm2(r,desc_a,info) - rerr = rni/bn2 + errnum = rni + errden = bn2 Endif - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If @@ -389,19 +400,25 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0)) & & write(*,'(a,i4,3(2x,es10.4))') & - & 'bicgstab: ',itx,rerr + & 'bicgstab: ',itx,errnum,eps*errden Endif End Do iteration End Do restart If (itrace_ > 0) then - if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',itx,rerr + if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',itx,errnum,eps*errden Endif - If (Present(err)) err=rerr + If (Present(err)) then + if (errden /= dzero) then + err = errnum/errden + else + err = errnum + end if + end If If (Present(iter)) iter = itx - If (rerr>eps) Then - write(debug_unit,*) 'BI-CGSTAB failed to converge to ',EPS,& + If ((errnum > eps*errden).and.(me==0)) Then + write(debug_unit,*) 'BI-CGSTAB failed to converge to ',eps*errden,& & ' in ',ITX,' iterations. ' End If diff --git a/krylov/psb_zrgmres.f90 b/krylov/psb_zrgmres.f90 index db67b2b1..1bb8b234 100644 --- a/krylov/psb_zrgmres.f90 +++ b/krylov/psb_zrgmres.f90 @@ -87,18 +87,22 @@ ! performed. ! iter - integer(optional) Output: how many iterations have been ! performed. -! err - real (optional) Output: error estimate on exit +! performed. +! err - real (optional) Output: error estimate on exit. If the +! denominator of the estimate is exactly +! 0, it is changed into 1. ! itrace - integer(optional) Input: print an informational message ! with the error estimate every itrace ! iterations -! irst - integer(optional) Input: restart parameter -! ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/|b| -! 2: err = |r|/(|a||x|+|b|) +! 1: err = |r|/|b|; here the iteration is +! stopped when |r| <= eps * |b| +! 2: err = |r|/(|a||x|+|b|); here the iteration is +! stopped when |r| <= eps * (|a||x|+|b|) ! where r is the (preconditioned, recursive -! estimate of) residual +! estimate of) residual. +! irst - integer(optional) Input: restart parameter ! Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_base_mod @@ -119,7 +123,7 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist !!$ local data complex(Kind(1.d0)), allocatable, target :: aux(:),w(:),w1(:), v(:,:) complex(Kind(1.d0)), allocatable :: c(:),s(:), h(:,:), rs(:),rst(:),xt(:) - Real(Kind(1.d0)) :: rerr, tmp + Real(Kind(1.d0)) :: tmp complex(kind(1.d0)) :: rti, rti1, scal Integer ::litmax, naux, mglob, it,k, itrace_,& & np,me, n_row, n_col, nl, int_err(5) @@ -129,6 +133,7 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist integer :: debug_level, debug_unit Real(Kind(1.d0)) :: rni, xni, bni, ani,bn2 real(kind(1.d0)), external :: dznrm2 + real(kind(1.d0)) :: errnum, errden character(len=20) :: name info = 0 @@ -240,6 +245,8 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist else if (istop_ == 2) then bn2 = psb_genrm2(b,desc_a,info) endif + errnum = dzero + errden = done if (info.ne.0) Then info=4011 call psb_errpush(info,name) @@ -290,10 +297,12 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist if (istop_ == 1) then rni = psb_geamax(v(:,1),desc_a,info) xni = psb_geamax(x,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) else if (istop_ == 2) then rni = psb_genrm2(v(:,1),desc_a,info) - rerr = rni/bn2 + errnum = rni + errden = bn2 endif if (info.ne.0) Then info=4011 @@ -301,13 +310,13 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist goto 9999 End If - If (rerr<=eps) Then + If (errnum <= eps*errden) Then Exit restart End If If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden end If v(:,1) = v(:,1) * scal @@ -363,7 +372,8 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist call psb_spmm(-zone,a,xt,zone,w1,desc_a,info,work=aux) rni = psb_geamax(w1,desc_a,info) xni = psb_geamax(xt,desc_a,info) - rerr = rni/(ani*xni+bni) + errnum = rni + errden = (ani*xni+bni) ! else if (istop_ == 2) then @@ -372,11 +382,11 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist ! procedure of the least-squares problem ! rni = abs(rs(i+1)) - rerr = rni/bn2 - + errnum = rni + errden = bn2 endif - if (rerr < eps ) then + If (errnum <= eps*errden) Then if (istop_ == 1) then x = xt @@ -402,7 +412,7 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist If (itrace_ > 0) then if ((mod(itx,itrace_)==0).and.(me == 0))& - & write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,rerr + & write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden end If end Do inner @@ -427,13 +437,19 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist End Do restart If (itrace_ > 0) then - if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,rerr + if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden end If - If (Present(err)) err=rerr + If (Present(err)) then + if (errden /= dzero) then + err = errnum/errden + else + err = errnum + end if + end If If (Present(iter)) iter = itx - If ((rerr>eps).and. (me == 0)) Then - write(debug_unit,*) 'rgmres(l) failed to converge to ',eps,& + If ((errnum > eps*errden).and.(me==0)) Then + write(debug_unit,*) 'rgmres(l) failed to converge to ',eps*errden,& & ' in ',itx,' iterations ' End If