From eabdeda1afed8b3c0bd066ff0a859df6b96ecd18 Mon Sep 17 00:00:00 2001 From: sfilippone Date: Fri, 22 May 2026 14:56:12 +0200 Subject: [PATCH] New stop criteria --- linsolve/impl/psb_cbicg.f90 | 34 +++++-- linsolve/impl/psb_ccg.F90 | 32 +++++- linsolve/impl/psb_ccgs.f90 | 33 +++++- linsolve/impl/psb_ccgstab.f90 | 38 +++++-- linsolve/impl/psb_ccgstabl.f90 | 33 +++++- linsolve/impl/psb_cfcg.F90 | 33 ++++-- linsolve/impl/psb_cgcr.f90 | 40 +++++--- linsolve/impl/psb_ckrylov.f90 | 30 +++--- linsolve/impl/psb_crgmres.f90 | 127 ++++++++++++++++++------ linsolve/impl/psb_crichardson.f90 | 11 +- linsolve/impl/psb_dbicg.f90 | 34 +++++-- linsolve/impl/psb_dcg.F90 | 32 +++++- linsolve/impl/psb_dcgs.f90 | 33 +++++- linsolve/impl/psb_dcgstab.f90 | 38 +++++-- linsolve/impl/psb_dcgstabl.f90 | 33 +++++- linsolve/impl/psb_dfcg.F90 | 33 ++++-- linsolve/impl/psb_dgcr.f90 | 40 +++++--- linsolve/impl/psb_dkrylov.f90 | 30 +++--- linsolve/impl/psb_drgmres.f90 | 127 ++++++++++++++++++------ linsolve/impl/psb_drichardson.f90 | 11 +- linsolve/impl/psb_sbicg.f90 | 34 +++++-- linsolve/impl/psb_scg.F90 | 32 +++++- linsolve/impl/psb_scgs.f90 | 33 +++++- linsolve/impl/psb_scgstab.f90 | 38 +++++-- linsolve/impl/psb_scgstabl.f90 | 33 +++++- linsolve/impl/psb_sfcg.F90 | 33 ++++-- linsolve/impl/psb_sgcr.f90 | 40 +++++--- linsolve/impl/psb_skrylov.f90 | 30 +++--- linsolve/impl/psb_srgmres.f90 | 127 ++++++++++++++++++------ linsolve/impl/psb_srichardson.f90 | 11 +- linsolve/impl/psb_zbicg.f90 | 34 +++++-- linsolve/impl/psb_zcg.F90 | 32 +++++- linsolve/impl/psb_zcgs.f90 | 33 +++++- linsolve/impl/psb_zcgstab.f90 | 38 +++++-- linsolve/impl/psb_zcgstabl.f90 | 33 +++++- linsolve/impl/psb_zfcg.F90 | 33 ++++-- linsolve/impl/psb_zgcr.f90 | 40 +++++--- linsolve/impl/psb_zkrylov.f90 | 30 +++--- linsolve/impl/psb_zrgmres.f90 | 127 ++++++++++++++++++------ linsolve/impl/psb_zrichardson.f90 | 11 +- linsolve/psb_base_linsolve_conv_mod.f90 | 33 +++++- linsolve/psb_c_linsolve_conv_mod.f90 | 68 ++++++++----- linsolve/psb_d_linsolve_conv_mod.f90 | 60 +++++++---- linsolve/psb_s_linsolve_conv_mod.f90 | 68 ++++++++----- linsolve/psb_z_linsolve_conv_mod.f90 | 60 +++++++---- 45 files changed, 1448 insertions(+), 485 deletions(-) diff --git a/linsolve/impl/psb_cbicg.f90 b/linsolve/impl/psb_cbicg.f90 index da9acad76..30d150aa6 100644 --- a/linsolve/impl/psb_cbicg.f90 +++ b/linsolve/impl/psb_cbicg.f90 @@ -95,7 +95,7 @@ ! subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop) + & itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_c_linsolve_conv_mod @@ -111,6 +111,7 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), optional, intent(in) :: itmax, itrace, istop integer(psb_ipk_), optional, intent(out) :: iter real(psb_spk_), optional, intent(out) :: err + type(psb_c_vect_type), intent(inout), optional :: s1, s2 ! !$ local data complex(psb_spk_), allocatable, target :: aux(:) type(psb_c_vect_type), allocatable, target :: wwrk(:) @@ -160,19 +161,29 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if ! ! istop_ = 1: normwise backward error, infinity norm - ! istop_ = 2: ||r||/||b|| norm 2 + ! istop_ = 2: ||r||/||b|| norm 2 ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then - info=psb_err_invalid_istop_ + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ err=info - call psb_errpush(info,name,i_err=(/istop_/)) + call psb_errpush(info,name,a_err="invalid istop_") goto 9999 - endif + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if(info /= psb_success_) then @@ -226,7 +237,8 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -252,7 +264,7 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,& rho = czero ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -306,7 +318,7 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(-alpha,q,cone,r,desc_a,info) call psb_geaxpby(-alpha,qt,cone,rt,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_ccg.F90 b/linsolve/impl/psb_ccg.F90 index 493fa3872..cf1ea1bea 100644 --- a/linsolve/impl/psb_ccg.F90 +++ b/linsolve/impl/psb_ccg.F90 @@ -96,7 +96,7 @@ ! ! subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop,cond) + & itmax,iter,err,itrace,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod use psb_c_linsolve_conv_mod @@ -112,6 +112,8 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err,cond + type(psb_c_vect_type), intent(inout), optional :: s1, s2 + ! = Local data complex(psb_spk_), allocatable, target :: aux(:),td(:),tu(:),eig(:),ewrk(:) integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:) @@ -159,8 +161,29 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& @@ -224,7 +247,8 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,& rho = czero - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + &desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -268,7 +292,7 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(alpha,p,cone,x,desc_a,info) call psb_geaxpby(-alpha,q,cone,r,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_ccgs.f90 b/linsolve/impl/psb_ccgs.f90 index 0c6e0673e..8783a0dad 100644 --- a/linsolve/impl/psb_ccgs.f90 +++ b/linsolve/impl/psb_ccgs.f90 @@ -93,7 +93,7 @@ ! estimate of) residual. ! Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop) + & itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_c_linsolve_conv_mod @@ -109,6 +109,7 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err + type(psb_c_vect_type), intent(inout), optional :: s1, s2 ! = local data complex(psb_spk_), allocatable, target :: aux(:) type(psb_c_vect_type), allocatable, target :: wwrk(:) @@ -154,8 +155,29 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,& If (Present(istop)) Then istop_ = istop Else - istop_ = 2 + istop_ = psb_get_istop_default() Endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) @@ -202,7 +224,8 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -225,7 +248,7 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,& ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -299,7 +322,7 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_ccgstab.f90 b/linsolve/impl/psb_ccgstab.f90 index 088d752c5..20c0a22a7 100644 --- a/linsolve/impl/psb_ccgstab.f90 +++ b/linsolve/impl/psb_ccgstab.f90 @@ -93,7 +93,7 @@ ! where r is the (preconditioned, recursive ! estimate of) residual. ! -Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) +Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_c_linsolve_conv_mod @@ -109,6 +109,7 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err + type(psb_c_vect_type), intent(inout), optional :: s1, s2 ! = Local data complex(psb_spk_), allocatable, target :: aux(:),wwrk(:,:) type(psb_c_vect_type) :: q, r, p, v, s, t, z, f @@ -156,13 +157,31 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist If (Present(istop)) Then istop_ = istop - Else - istop_ = 2 - Endif + else + istop_ = psb_get_istop_default() + endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if ! - ! ISTOP_ = 1: Normwise backward error, infinity norm - ! ISTOP_ = 2: ||r||/||b|| norm 2 + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select + ! = if (.not.same_type_as(x,b)) then ! = write(0,*) 'Warning: different dynamic types for X and B ' ! = end if @@ -217,7 +236,8 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist End If itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -234,7 +254,7 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist call psb_geaxpby(cone,r,czero,q,desc_a,info) ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (psb_errstatus_fatal()) then info=psb_err_from_subroutine_ @@ -354,7 +374,7 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist call psb_geaxpby(omega,z,cone,x,desc_a,info) call psb_geaxpby(cone,s,czero,r,desc_a,info) call psb_geaxpby(-omega,t,cone,r,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ') diff --git a/linsolve/impl/psb_ccgstabl.f90 b/linsolve/impl/psb_ccgstabl.f90 index ac713cac5..b8898c3aa 100644 --- a/linsolve/impl/psb_ccgstabl.f90 +++ b/linsolve/impl/psb_ccgstabl.f90 @@ -104,7 +104,7 @@ ! ! Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop) + & itmax,iter,err,itrace,irst,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_c_linsolve_conv_mod @@ -120,6 +120,7 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err + type(psb_c_vect_type), intent(inout), optional :: s1, s2 ! = local data complex(psb_spk_), allocatable, target :: aux(:), gamma(:),& & gamma1(:), gamma2(:), taum(:,:), sigma(:) @@ -172,8 +173,29 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then itmax_ = itmax @@ -246,7 +268,8 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,& rt0 => wwrk(10) - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -284,7 +307,7 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,& & write(debug_unit,*) me,' ',trim(name),& & ' on entry to amax: b: ',b%get_nrows() - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -388,7 +411,7 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(-gamma1(j),rh(j),cone,rh(0),desc_a,info) enddo - if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_cfcg.F90 b/linsolve/impl/psb_cfcg.F90 index c77191d09..252e19e8e 100644 --- a/linsolve/impl/psb_cfcg.F90 +++ b/linsolve/impl/psb_cfcg.F90 @@ -104,7 +104,7 @@ ! ! subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop,cond) + & itmax,iter,err,itrace,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod use psb_c_linsolve_conv_mod @@ -120,6 +120,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_spk_), Optional, Intent(out) :: err,cond + type(psb_c_vect_type), intent(inout), optional :: s1, s2 ! = Local data type(psb_c_vect_type) :: v, w, d , q, r complex(psb_spk_) :: alpha, beta, delta, gamma, theta @@ -164,9 +165,29 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif - + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& @@ -207,7 +228,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,& & scratch=.true.,mold=x%v) call psb_init_conv(methdname,istop_,itrace_,itmax_,& - & a,x,b,eps,desc_a,stopdat,info) + & a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2) itx = 0 restart: do @@ -226,7 +247,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,& end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from restart' exit restart end if @@ -282,7 +303,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,& itx = itx + 1 - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from iteration' exit restart end if diff --git a/linsolve/impl/psb_cgcr.f90 b/linsolve/impl/psb_cgcr.f90 index 7201df51d..2552003c6 100644 --- a/linsolve/impl/psb_cgcr.f90 +++ b/linsolve/impl/psb_cgcr.f90 @@ -106,7 +106,7 @@ ! subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace, irst, istop) + & itmax,iter,err,itrace, irst, istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_c_linsolve_conv_mod @@ -124,6 +124,7 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_spk_), Optional, Intent(out) :: err + type(psb_c_vect_type), intent(inout), optional :: s1, s2 ! = local data complex(psb_spk_), allocatable :: alpha(:), h(:,:) type(psb_c_vect_type), allocatable :: z(:), c(:), c_scale(:) @@ -167,22 +168,30 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif - - ! - ! ISTOP_ = 1: Normwise backward error, infinity norm - ! ISTOP_ = 2: ||r||/||b||, 2-norm - ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif - - + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select + call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& & call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) @@ -245,7 +254,8 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 nrst = -1 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) restart: do if (itx>= itmax_) exit restart h = czero @@ -268,7 +278,7 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart nrst = nrst + 1 @@ -299,7 +309,7 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(cone, r, czero, r, desc_a, info) call psb_geaxpby(-alpha(j), c_scale(j), cone, r, desc_a, info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (j >= irst) exit iteration diff --git a/linsolve/impl/psb_ckrylov.f90 b/linsolve/impl/psb_ckrylov.f90 index 22712da05..a465eaae1 100644 --- a/linsolve/impl/psb_ckrylov.f90 +++ b/linsolve/impl/psb_ckrylov.f90 @@ -80,7 +80,7 @@ ! estimate of) residual ! Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop,cond) + & itmax,iter,err,itrace,irst,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod,only : psb_cprec_type @@ -97,11 +97,12 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err,cond + type(psb_c_vect_type), intent(inout), optional :: s1, s2 abstract interface subroutine psb_ckryl_vect(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) + & desc_a,info,itmax,iter,err,itrace,istop,s1,s2) import :: psb_ipk_, psb_spk_, psb_desc_type, & & psb_cspmat_type, psb_cprec_type, psb_c_vect_type type(psb_cspmat_type), intent(in) :: a @@ -114,9 +115,10 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), optional, intent(in) :: itmax, itrace,istop integer(psb_ipk_), optional, intent(out) :: iter real(psb_spk_), optional, intent(out) :: err + type(psb_c_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_ckryl_vect Subroutine psb_ckryl_rest_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,irst,istop) + &itmax,iter,err, itrace,irst,istop,s1,s2) import :: psb_ipk_, psb_spk_, psb_desc_type, & & psb_cspmat_type, psb_cprec_type, psb_c_vect_type Type(psb_cspmat_type), Intent(in) :: a @@ -129,9 +131,10 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err + type(psb_c_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_ckryl_rest_vect Subroutine psb_ckryl_cond_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,istop,cond) + &itmax,iter,err, itrace,istop,cond,s1,s2) import :: psb_ipk_, psb_spk_, psb_desc_type, & & psb_cspmat_type, psb_cprec_type, psb_c_vect_type Type(psb_cspmat_type), Intent(in) :: a @@ -144,6 +147,7 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err, cond + type(psb_c_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_ckryl_cond_vect end interface @@ -180,34 +184,34 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,& select case(psb_toupper(method)) case('CG') call psb_ccg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond) + &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2) case('FCG') call psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond) + &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2) case('GCR') call psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('CGS') call psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('BICG') call psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('BICGSTAB') call psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('RGMRES','GMRES') call psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop) + &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2) case('BICGSTABL') call psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop) + &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2) case default if (me == 0) write(psb_err_unit,*) trim(name),& & ': Warning: Unknown method ',method,& & ', defaulting to BiCGSTAB' call psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) end select if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) diff --git a/linsolve/impl/psb_crgmres.f90 b/linsolve/impl/psb_crgmres.f90 index 673d13994..dc2fe6e2f 100644 --- a/linsolve/impl/psb_crgmres.f90 +++ b/linsolve/impl/psb_crgmres.f90 @@ -97,17 +97,20 @@ ! iterations ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/(|a||x|+|b|); here the iteration is +! 1: err = |r|/(|a||x|+|b|); here +! the iteration is ! stopped when |r| <= eps * (|a||x|+|b|) ! 2: err = |r|/|b|; here the iteration is ! stopped when |r| <= eps * |b| +! 3: Same as 2 but with X and B scaled +! by s1 and s2 ! where r is the (preconditioned, recursive ! estimate of) residual. ! irst - integer(optional) Input: restart parameter ! subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop) + & itmax,iter,err,itrace,irst,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_c_linsolve_conv_mod @@ -123,6 +126,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err + type(psb_c_vect_type), intent(inout), optional :: s1, s2 ! = local data complex(psb_spk_), allocatable :: aux(:) complex(psb_spk_), allocatable :: c(:), s(:), h(:,:), rs(:), rst(:) @@ -172,19 +176,30 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_,psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then litmax = itmax @@ -253,12 +268,15 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& & w%get_nrows(),w1%get_nrows() - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ani = psb_spnrmi(a,desc_a,info) bni = psb_geamax(b,desc_a,info) - else if (istop_ == 2) then - bn2 = psb_genrm2(b,desc_a,info) - else if (istop_ == 3) then + case(psb_istop_bn2_) + bn2 = psb_genrm2(b,desc_a,info) + case(psb_istop_rn2_abs_) + ! do nothing + case(psb_istop_rrn2_) call psb_geaxpby(cone,b,czero,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -273,7 +291,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if r0n2 = psb_genrm2(v(1),desc_a,info) - endif + end select + errnum = czero errden = cone deps = eps @@ -307,7 +326,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 end if - + if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info) + rs(1) = psb_genrm2(v(1),desc_a,info) rs(2:) = czero if (info /= psb_success_) then @@ -324,20 +344,25 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& ! ! check convergence ! - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) rni = psb_geamax(v(1),desc_a,info) xni = psb_geamax(x,desc_a,info) errnum = rni errden = (ani*xni+bni) - else if (istop_ == 2) then + case(psb_istop_bn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 - else if (istop_ == 3) then + case(psb_istop_rn2_abs_) + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = sone + case(psb_istop_rrn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = r0n2 - endif + end select if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) @@ -381,7 +406,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& h(i+1,i) = czero call crot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_spk_),s(i)) - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ! ! build x and then compute the residual and its infinity norm ! @@ -403,8 +429,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& errnum = rni errden = (ani*xni+bni) ! - - else if (istop_ == 2) then + case(psb_istop_bn2_) ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -412,7 +437,14 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = bn2 - else if (istop_ == 3) then + + case(psb_istop_rn2_abs_) + rni = abs(rs(i+1)) + errnum = rni + errden = sone + + case(psb_istop_rrn2_) + ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -420,14 +452,15 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = r0n2 - endif + end select if (errnum <= eps*errden) then - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(cone,xt,czero,x,desc_a,info) ! = x = xt - else if (istop_ == 2) then + case(psb_istop_bn2_, psb_istop_rn2_abs_) ! ! build x ! @@ -441,7 +474,22 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& end do call prec%apply(w1,w,desc_a,info) call psb_geaxpby(cone,w,cone,x,desc_a,info) - end if + + case(psb_istop_rrn2_) + ! + ! build x + ! + call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:i) + call w1%set(czero) + do k=1, i + call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(cone,w,cone,x,desc_a,info) + end select if (itrace_ > 0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -454,9 +502,11 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& end do inner - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(cone,xt,czero,x,desc_a,info)! x = xt - else if (istop_ == 2) then + + case(psb_istop_bn2_, psb_istop_rn2_abs_) ! ! build x ! @@ -470,7 +520,22 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& end do call prec%apply(w1,w,desc_a,info) call psb_geaxpby(cone,w,cone,x,desc_a,info) - end if + + case(psb_istop_rrn2_) + ! + ! build x + ! + call ctrsm('l','u','n','n',nl,1,cone,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:nl) + call w1%set(czero) + do k=1, nl + call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(cone,w,cone,x,desc_a,info) + end select if (itx >= litmax) then if (itrace_ > 0) then diff --git a/linsolve/impl/psb_crichardson.f90 b/linsolve/impl/psb_crichardson.f90 index 2ff18eb68..270b63ad1 100644 --- a/linsolve/impl/psb_crichardson.f90 +++ b/linsolve/impl/psb_crichardson.f90 @@ -120,8 +120,17 @@ Subroutine psb_crichardson_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + + if (present(itmax)) then itmax_ = itmax else diff --git a/linsolve/impl/psb_dbicg.f90 b/linsolve/impl/psb_dbicg.f90 index fef0fbc2e..dc8018800 100644 --- a/linsolve/impl/psb_dbicg.f90 +++ b/linsolve/impl/psb_dbicg.f90 @@ -95,7 +95,7 @@ ! subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop) + & itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_d_linsolve_conv_mod @@ -111,6 +111,7 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), optional, intent(in) :: itmax, itrace, istop integer(psb_ipk_), optional, intent(out) :: iter real(psb_dpk_), optional, intent(out) :: err + type(psb_d_vect_type), intent(inout), optional :: s1, s2 ! !$ local data real(psb_dpk_), allocatable, target :: aux(:) type(psb_d_vect_type), allocatable, target :: wwrk(:) @@ -160,19 +161,29 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if ! ! istop_ = 1: normwise backward error, infinity norm - ! istop_ = 2: ||r||/||b|| norm 2 + ! istop_ = 2: ||r||/||b|| norm 2 ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then - info=psb_err_invalid_istop_ + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ err=info - call psb_errpush(info,name,i_err=(/istop_/)) + call psb_errpush(info,name,a_err="invalid istop_") goto 9999 - endif + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if(info /= psb_success_) then @@ -226,7 +237,8 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -252,7 +264,7 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,& rho = dzero ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -306,7 +318,7 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(-alpha,q,done,r,desc_a,info) call psb_geaxpby(-alpha,qt,done,rt,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_dcg.F90 b/linsolve/impl/psb_dcg.F90 index 1fea599ef..a60742b58 100644 --- a/linsolve/impl/psb_dcg.F90 +++ b/linsolve/impl/psb_dcg.F90 @@ -96,7 +96,7 @@ ! ! subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop,cond) + & itmax,iter,err,itrace,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod use psb_d_linsolve_conv_mod @@ -112,6 +112,8 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err,cond + type(psb_d_vect_type), intent(inout), optional :: s1, s2 + ! = Local data real(psb_dpk_), allocatable, target :: aux(:),td(:),tu(:),eig(:),ewrk(:) integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:) @@ -159,8 +161,29 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& @@ -232,7 +255,8 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,& rho = dzero - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + &desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -285,7 +309,7 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(alpha,p,done,x,desc_a,info) call psb_geaxpby(-alpha,q,done,r,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_dcgs.f90 b/linsolve/impl/psb_dcgs.f90 index 37dde891a..60173a057 100644 --- a/linsolve/impl/psb_dcgs.f90 +++ b/linsolve/impl/psb_dcgs.f90 @@ -93,7 +93,7 @@ ! estimate of) residual. ! Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop) + & itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_d_linsolve_conv_mod @@ -109,6 +109,7 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err + type(psb_d_vect_type), intent(inout), optional :: s1, s2 ! = local data real(psb_dpk_), allocatable, target :: aux(:) type(psb_d_vect_type), allocatable, target :: wwrk(:) @@ -154,8 +155,29 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,& If (Present(istop)) Then istop_ = istop Else - istop_ = 2 + istop_ = psb_get_istop_default() Endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) @@ -202,7 +224,8 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -225,7 +248,7 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,& ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -299,7 +322,7 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_dcgstab.f90 b/linsolve/impl/psb_dcgstab.f90 index 20c5e31d2..95633379b 100644 --- a/linsolve/impl/psb_dcgstab.f90 +++ b/linsolve/impl/psb_dcgstab.f90 @@ -93,7 +93,7 @@ ! where r is the (preconditioned, recursive ! estimate of) residual. ! -Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) +Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_d_linsolve_conv_mod @@ -109,6 +109,7 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err + type(psb_d_vect_type), intent(inout), optional :: s1, s2 ! = Local data real(psb_dpk_), allocatable, target :: aux(:),wwrk(:,:) type(psb_d_vect_type) :: q, r, p, v, s, t, z, f @@ -156,13 +157,31 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist If (Present(istop)) Then istop_ = istop - Else - istop_ = 2 - Endif + else + istop_ = psb_get_istop_default() + endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if ! - ! ISTOP_ = 1: Normwise backward error, infinity norm - ! ISTOP_ = 2: ||r||/||b|| norm 2 + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select + ! = if (.not.same_type_as(x,b)) then ! = write(0,*) 'Warning: different dynamic types for X and B ' ! = end if @@ -217,7 +236,8 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist End If itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -234,7 +254,7 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist call psb_geaxpby(done,r,dzero,q,desc_a,info) ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (psb_errstatus_fatal()) then info=psb_err_from_subroutine_ @@ -354,7 +374,7 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist call psb_geaxpby(omega,z,done,x,desc_a,info) call psb_geaxpby(done,s,dzero,r,desc_a,info) call psb_geaxpby(-omega,t,done,r,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ') diff --git a/linsolve/impl/psb_dcgstabl.f90 b/linsolve/impl/psb_dcgstabl.f90 index c5a9d157f..152dffdc8 100644 --- a/linsolve/impl/psb_dcgstabl.f90 +++ b/linsolve/impl/psb_dcgstabl.f90 @@ -104,7 +104,7 @@ ! ! Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop) + & itmax,iter,err,itrace,irst,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_d_linsolve_conv_mod @@ -120,6 +120,7 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err + type(psb_d_vect_type), intent(inout), optional :: s1, s2 ! = local data real(psb_dpk_), allocatable, target :: aux(:), gamma(:),& & gamma1(:), gamma2(:), taum(:,:), sigma(:) @@ -172,8 +173,29 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then itmax_ = itmax @@ -246,7 +268,8 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,& rt0 => wwrk(10) - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -284,7 +307,7 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,& & write(debug_unit,*) me,' ',trim(name),& & ' on entry to amax: b: ',b%get_nrows() - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -388,7 +411,7 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(-gamma1(j),rh(j),done,rh(0),desc_a,info) enddo - if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_dfcg.F90 b/linsolve/impl/psb_dfcg.F90 index 157e7f76a..9e0d62213 100644 --- a/linsolve/impl/psb_dfcg.F90 +++ b/linsolve/impl/psb_dfcg.F90 @@ -104,7 +104,7 @@ ! ! subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop,cond) + & itmax,iter,err,itrace,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod use psb_d_linsolve_conv_mod @@ -120,6 +120,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_dpk_), Optional, Intent(out) :: err,cond + type(psb_d_vect_type), intent(inout), optional :: s1, s2 ! = Local data type(psb_d_vect_type) :: v, w, d , q, r real(psb_dpk_) :: alpha, beta, delta, gamma, theta @@ -164,9 +165,29 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif - + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& @@ -207,7 +228,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,& & scratch=.true.,mold=x%v) call psb_init_conv(methdname,istop_,itrace_,itmax_,& - & a,x,b,eps,desc_a,stopdat,info) + & a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2) itx = 0 restart: do @@ -226,7 +247,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,& end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from restart' exit restart end if @@ -282,7 +303,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,& itx = itx + 1 - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from iteration' exit restart end if diff --git a/linsolve/impl/psb_dgcr.f90 b/linsolve/impl/psb_dgcr.f90 index 8d59d83ab..43400a9aa 100644 --- a/linsolve/impl/psb_dgcr.f90 +++ b/linsolve/impl/psb_dgcr.f90 @@ -106,7 +106,7 @@ ! subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace, irst, istop) + & itmax,iter,err,itrace, irst, istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_d_linsolve_conv_mod @@ -124,6 +124,7 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_dpk_), Optional, Intent(out) :: err + type(psb_d_vect_type), intent(inout), optional :: s1, s2 ! = local data real(psb_dpk_), allocatable :: alpha(:), h(:,:) type(psb_d_vect_type), allocatable :: z(:), c(:), c_scale(:) @@ -167,22 +168,30 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif - - ! - ! ISTOP_ = 1: Normwise backward error, infinity norm - ! ISTOP_ = 2: ||r||/||b||, 2-norm - ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif - - + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select + call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& & call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) @@ -245,7 +254,8 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 nrst = -1 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) restart: do if (itx>= itmax_) exit restart h = dzero @@ -268,7 +278,7 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart nrst = nrst + 1 @@ -299,7 +309,7 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(done, r, dzero, r, desc_a, info) call psb_geaxpby(-alpha(j), c_scale(j), done, r, desc_a, info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (j >= irst) exit iteration diff --git a/linsolve/impl/psb_dkrylov.f90 b/linsolve/impl/psb_dkrylov.f90 index 5963f5a8a..3652cc934 100644 --- a/linsolve/impl/psb_dkrylov.f90 +++ b/linsolve/impl/psb_dkrylov.f90 @@ -80,7 +80,7 @@ ! estimate of) residual ! Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop,cond) + & itmax,iter,err,itrace,irst,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod,only : psb_dprec_type @@ -97,11 +97,12 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err,cond + type(psb_d_vect_type), intent(inout), optional :: s1, s2 abstract interface subroutine psb_dkryl_vect(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) + & desc_a,info,itmax,iter,err,itrace,istop,s1,s2) import :: psb_ipk_, psb_dpk_, psb_desc_type, & & psb_dspmat_type, psb_dprec_type, psb_d_vect_type type(psb_dspmat_type), intent(in) :: a @@ -114,9 +115,10 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), optional, intent(in) :: itmax, itrace,istop integer(psb_ipk_), optional, intent(out) :: iter real(psb_dpk_), optional, intent(out) :: err + type(psb_d_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_dkryl_vect Subroutine psb_dkryl_rest_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,irst,istop) + &itmax,iter,err, itrace,irst,istop,s1,s2) import :: psb_ipk_, psb_dpk_, psb_desc_type, & & psb_dspmat_type, psb_dprec_type, psb_d_vect_type Type(psb_dspmat_type), Intent(in) :: a @@ -129,9 +131,10 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err + type(psb_d_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_dkryl_rest_vect Subroutine psb_dkryl_cond_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,istop,cond) + &itmax,iter,err, itrace,istop,cond,s1,s2) import :: psb_ipk_, psb_dpk_, psb_desc_type, & & psb_dspmat_type, psb_dprec_type, psb_d_vect_type Type(psb_dspmat_type), Intent(in) :: a @@ -144,6 +147,7 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err, cond + type(psb_d_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_dkryl_cond_vect end interface @@ -180,34 +184,34 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& select case(psb_toupper(method)) case('CG') call psb_dcg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond) + &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2) case('FCG') call psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond) + &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2) case('GCR') call psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('CGS') call psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('BICG') call psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('BICGSTAB') call psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('RGMRES','GMRES') call psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop) + &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2) case('BICGSTABL') call psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop) + &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2) case default if (me == 0) write(psb_err_unit,*) trim(name),& & ': Warning: Unknown method ',method,& & ', defaulting to BiCGSTAB' call psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) end select if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) diff --git a/linsolve/impl/psb_drgmres.f90 b/linsolve/impl/psb_drgmres.f90 index e03e043f7..712ddebc8 100644 --- a/linsolve/impl/psb_drgmres.f90 +++ b/linsolve/impl/psb_drgmres.f90 @@ -97,17 +97,20 @@ ! iterations ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/(|a||x|+|b|); here the iteration is +! 1: err = |r|/(|a||x|+|b|); here +! the iteration is ! stopped when |r| <= eps * (|a||x|+|b|) ! 2: err = |r|/|b|; here the iteration is ! stopped when |r| <= eps * |b| +! 3: Same as 2 but with X and B scaled +! by s1 and s2 ! where r is the (preconditioned, recursive ! estimate of) residual. ! irst - integer(optional) Input: restart parameter ! subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop) + & itmax,iter,err,itrace,irst,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_d_linsolve_conv_mod @@ -123,6 +126,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err + type(psb_d_vect_type), intent(inout), optional :: s1, s2 ! = local data real(psb_dpk_), allocatable :: aux(:) real(psb_dpk_), allocatable :: c(:), s(:), h(:,:), rs(:), rst(:) @@ -172,19 +176,30 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_,psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then litmax = itmax @@ -253,12 +268,15 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& & w%get_nrows(),w1%get_nrows() - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ani = psb_spnrmi(a,desc_a,info) bni = psb_geamax(b,desc_a,info) - else if (istop_ == 2) then - bn2 = psb_genrm2(b,desc_a,info) - else if (istop_ == 3) then + case(psb_istop_bn2_) + bn2 = psb_genrm2(b,desc_a,info) + case(psb_istop_rn2_abs_) + ! do nothing + case(psb_istop_rrn2_) call psb_geaxpby(done,b,dzero,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -273,7 +291,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if r0n2 = psb_genrm2(v(1),desc_a,info) - endif + end select + errnum = dzero errden = done deps = eps @@ -307,7 +326,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 end if - + if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info) + rs(1) = psb_genrm2(v(1),desc_a,info) rs(2:) = dzero if (info /= psb_success_) then @@ -324,20 +344,25 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& ! ! check convergence ! - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) rni = psb_geamax(v(1),desc_a,info) xni = psb_geamax(x,desc_a,info) errnum = rni errden = (ani*xni+bni) - else if (istop_ == 2) then + case(psb_istop_bn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 - else if (istop_ == 3) then + case(psb_istop_rn2_abs_) + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = done + case(psb_istop_rrn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = r0n2 - endif + end select if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) @@ -381,7 +406,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& h(i+1,i) = dzero call drot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_dpk_),s(i)) - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ! ! build x and then compute the residual and its infinity norm ! @@ -403,8 +429,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& errnum = rni errden = (ani*xni+bni) ! - - else if (istop_ == 2) then + case(psb_istop_bn2_) ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -412,7 +437,14 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = bn2 - else if (istop_ == 3) then + + case(psb_istop_rn2_abs_) + rni = abs(rs(i+1)) + errnum = rni + errden = done + + case(psb_istop_rrn2_) + ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -420,14 +452,15 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = r0n2 - endif + end select if (errnum <= eps*errden) then - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(done,xt,dzero,x,desc_a,info) ! = x = xt - else if (istop_ == 2) then + case(psb_istop_bn2_, psb_istop_rn2_abs_) ! ! build x ! @@ -441,7 +474,22 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& end do call prec%apply(w1,w,desc_a,info) call psb_geaxpby(done,w,done,x,desc_a,info) - end if + + case(psb_istop_rrn2_) + ! + ! build x + ! + call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:i) + call w1%set(dzero) + do k=1, i + call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(done,w,done,x,desc_a,info) + end select if (itrace_ > 0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -454,9 +502,11 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& end do inner - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(done,xt,dzero,x,desc_a,info)! x = xt - else if (istop_ == 2) then + + case(psb_istop_bn2_, psb_istop_rn2_abs_) ! ! build x ! @@ -470,7 +520,22 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& end do call prec%apply(w1,w,desc_a,info) call psb_geaxpby(done,w,done,x,desc_a,info) - end if + + case(psb_istop_rrn2_) + ! + ! build x + ! + call dtrsm('l','u','n','n',nl,1,done,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:nl) + call w1%set(dzero) + do k=1, nl + call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(done,w,done,x,desc_a,info) + end select if (itx >= litmax) then if (itrace_ > 0) then diff --git a/linsolve/impl/psb_drichardson.f90 b/linsolve/impl/psb_drichardson.f90 index 70b60e650..cd350cf38 100644 --- a/linsolve/impl/psb_drichardson.f90 +++ b/linsolve/impl/psb_drichardson.f90 @@ -120,8 +120,17 @@ Subroutine psb_drichardson_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + + if (present(itmax)) then itmax_ = itmax else diff --git a/linsolve/impl/psb_sbicg.f90 b/linsolve/impl/psb_sbicg.f90 index 46cde1319..88d81d977 100644 --- a/linsolve/impl/psb_sbicg.f90 +++ b/linsolve/impl/psb_sbicg.f90 @@ -95,7 +95,7 @@ ! subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop) + & itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_s_linsolve_conv_mod @@ -111,6 +111,7 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), optional, intent(in) :: itmax, itrace, istop integer(psb_ipk_), optional, intent(out) :: iter real(psb_spk_), optional, intent(out) :: err + type(psb_s_vect_type), intent(inout), optional :: s1, s2 ! !$ local data real(psb_spk_), allocatable, target :: aux(:) type(psb_s_vect_type), allocatable, target :: wwrk(:) @@ -160,19 +161,29 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if ! ! istop_ = 1: normwise backward error, infinity norm - ! istop_ = 2: ||r||/||b|| norm 2 + ! istop_ = 2: ||r||/||b|| norm 2 ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then - info=psb_err_invalid_istop_ + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ err=info - call psb_errpush(info,name,i_err=(/istop_/)) + call psb_errpush(info,name,a_err="invalid istop_") goto 9999 - endif + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if(info /= psb_success_) then @@ -226,7 +237,8 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -252,7 +264,7 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,& rho = szero ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -306,7 +318,7 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(-alpha,q,sone,r,desc_a,info) call psb_geaxpby(-alpha,qt,sone,rt,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_scg.F90 b/linsolve/impl/psb_scg.F90 index b199284f0..b93603660 100644 --- a/linsolve/impl/psb_scg.F90 +++ b/linsolve/impl/psb_scg.F90 @@ -96,7 +96,7 @@ ! ! subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop,cond) + & itmax,iter,err,itrace,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod use psb_s_linsolve_conv_mod @@ -112,6 +112,8 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err,cond + type(psb_s_vect_type), intent(inout), optional :: s1, s2 + ! = Local data real(psb_spk_), allocatable, target :: aux(:),td(:),tu(:),eig(:),ewrk(:) integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:) @@ -159,8 +161,29 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& @@ -232,7 +255,8 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,& rho = szero - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + &desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -285,7 +309,7 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(alpha,p,sone,x,desc_a,info) call psb_geaxpby(-alpha,q,sone,r,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_scgs.f90 b/linsolve/impl/psb_scgs.f90 index 084b1f743..1e13d39a4 100644 --- a/linsolve/impl/psb_scgs.f90 +++ b/linsolve/impl/psb_scgs.f90 @@ -93,7 +93,7 @@ ! estimate of) residual. ! Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop) + & itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_s_linsolve_conv_mod @@ -109,6 +109,7 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err + type(psb_s_vect_type), intent(inout), optional :: s1, s2 ! = local data real(psb_spk_), allocatable, target :: aux(:) type(psb_s_vect_type), allocatable, target :: wwrk(:) @@ -154,8 +155,29 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,& If (Present(istop)) Then istop_ = istop Else - istop_ = 2 + istop_ = psb_get_istop_default() Endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) @@ -202,7 +224,8 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -225,7 +248,7 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,& ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -299,7 +322,7 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_scgstab.f90 b/linsolve/impl/psb_scgstab.f90 index 813bb2b5b..f1e72b73b 100644 --- a/linsolve/impl/psb_scgstab.f90 +++ b/linsolve/impl/psb_scgstab.f90 @@ -93,7 +93,7 @@ ! where r is the (preconditioned, recursive ! estimate of) residual. ! -Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) +Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_s_linsolve_conv_mod @@ -109,6 +109,7 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err + type(psb_s_vect_type), intent(inout), optional :: s1, s2 ! = Local data real(psb_spk_), allocatable, target :: aux(:),wwrk(:,:) type(psb_s_vect_type) :: q, r, p, v, s, t, z, f @@ -156,13 +157,31 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist If (Present(istop)) Then istop_ = istop - Else - istop_ = 2 - Endif + else + istop_ = psb_get_istop_default() + endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if ! - ! ISTOP_ = 1: Normwise backward error, infinity norm - ! ISTOP_ = 2: ||r||/||b|| norm 2 + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select + ! = if (.not.same_type_as(x,b)) then ! = write(0,*) 'Warning: different dynamic types for X and B ' ! = end if @@ -217,7 +236,8 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist End If itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -234,7 +254,7 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist call psb_geaxpby(sone,r,szero,q,desc_a,info) ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (psb_errstatus_fatal()) then info=psb_err_from_subroutine_ @@ -354,7 +374,7 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist call psb_geaxpby(omega,z,sone,x,desc_a,info) call psb_geaxpby(sone,s,szero,r,desc_a,info) call psb_geaxpby(-omega,t,sone,r,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ') diff --git a/linsolve/impl/psb_scgstabl.f90 b/linsolve/impl/psb_scgstabl.f90 index adf304a6f..a9d5e783a 100644 --- a/linsolve/impl/psb_scgstabl.f90 +++ b/linsolve/impl/psb_scgstabl.f90 @@ -104,7 +104,7 @@ ! ! Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop) + & itmax,iter,err,itrace,irst,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_s_linsolve_conv_mod @@ -120,6 +120,7 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err + type(psb_s_vect_type), intent(inout), optional :: s1, s2 ! = local data real(psb_spk_), allocatable, target :: aux(:), gamma(:),& & gamma1(:), gamma2(:), taum(:,:), sigma(:) @@ -172,8 +173,29 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then itmax_ = itmax @@ -246,7 +268,8 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,& rt0 => wwrk(10) - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -284,7 +307,7 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,& & write(debug_unit,*) me,' ',trim(name),& & ' on entry to amax: b: ',b%get_nrows() - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -388,7 +411,7 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(-gamma1(j),rh(j),sone,rh(0),desc_a,info) enddo - if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_sfcg.F90 b/linsolve/impl/psb_sfcg.F90 index 9453a0651..e55cd6a8b 100644 --- a/linsolve/impl/psb_sfcg.F90 +++ b/linsolve/impl/psb_sfcg.F90 @@ -104,7 +104,7 @@ ! ! subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop,cond) + & itmax,iter,err,itrace,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod use psb_s_linsolve_conv_mod @@ -120,6 +120,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_spk_), Optional, Intent(out) :: err,cond + type(psb_s_vect_type), intent(inout), optional :: s1, s2 ! = Local data type(psb_s_vect_type) :: v, w, d , q, r real(psb_spk_) :: alpha, beta, delta, gamma, theta @@ -164,9 +165,29 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif - + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& @@ -207,7 +228,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,& & scratch=.true.,mold=x%v) call psb_init_conv(methdname,istop_,itrace_,itmax_,& - & a,x,b,eps,desc_a,stopdat,info) + & a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2) itx = 0 restart: do @@ -226,7 +247,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,& end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from restart' exit restart end if @@ -282,7 +303,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,& itx = itx + 1 - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from iteration' exit restart end if diff --git a/linsolve/impl/psb_sgcr.f90 b/linsolve/impl/psb_sgcr.f90 index bee8bdff7..8966c280d 100644 --- a/linsolve/impl/psb_sgcr.f90 +++ b/linsolve/impl/psb_sgcr.f90 @@ -106,7 +106,7 @@ ! subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace, irst, istop) + & itmax,iter,err,itrace, irst, istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_s_linsolve_conv_mod @@ -124,6 +124,7 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_spk_), Optional, Intent(out) :: err + type(psb_s_vect_type), intent(inout), optional :: s1, s2 ! = local data real(psb_spk_), allocatable :: alpha(:), h(:,:) type(psb_s_vect_type), allocatable :: z(:), c(:), c_scale(:) @@ -167,22 +168,30 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif - - ! - ! ISTOP_ = 1: Normwise backward error, infinity norm - ! ISTOP_ = 2: ||r||/||b||, 2-norm - ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif - - + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select + call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& & call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) @@ -245,7 +254,8 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 nrst = -1 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) restart: do if (itx>= itmax_) exit restart h = szero @@ -268,7 +278,7 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart nrst = nrst + 1 @@ -299,7 +309,7 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(sone, r, szero, r, desc_a, info) call psb_geaxpby(-alpha(j), c_scale(j), sone, r, desc_a, info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (j >= irst) exit iteration diff --git a/linsolve/impl/psb_skrylov.f90 b/linsolve/impl/psb_skrylov.f90 index e24d7d05f..985c50c18 100644 --- a/linsolve/impl/psb_skrylov.f90 +++ b/linsolve/impl/psb_skrylov.f90 @@ -80,7 +80,7 @@ ! estimate of) residual ! Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop,cond) + & itmax,iter,err,itrace,irst,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod,only : psb_sprec_type @@ -97,11 +97,12 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err,cond + type(psb_s_vect_type), intent(inout), optional :: s1, s2 abstract interface subroutine psb_skryl_vect(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) + & desc_a,info,itmax,iter,err,itrace,istop,s1,s2) import :: psb_ipk_, psb_spk_, psb_desc_type, & & psb_sspmat_type, psb_sprec_type, psb_s_vect_type type(psb_sspmat_type), intent(in) :: a @@ -114,9 +115,10 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), optional, intent(in) :: itmax, itrace,istop integer(psb_ipk_), optional, intent(out) :: iter real(psb_spk_), optional, intent(out) :: err + type(psb_s_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_skryl_vect Subroutine psb_skryl_rest_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,irst,istop) + &itmax,iter,err, itrace,irst,istop,s1,s2) import :: psb_ipk_, psb_spk_, psb_desc_type, & & psb_sspmat_type, psb_sprec_type, psb_s_vect_type Type(psb_sspmat_type), Intent(in) :: a @@ -129,9 +131,10 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err + type(psb_s_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_skryl_rest_vect Subroutine psb_skryl_cond_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,istop,cond) + &itmax,iter,err, itrace,istop,cond,s1,s2) import :: psb_ipk_, psb_spk_, psb_desc_type, & & psb_sspmat_type, psb_sprec_type, psb_s_vect_type Type(psb_sspmat_type), Intent(in) :: a @@ -144,6 +147,7 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err, cond + type(psb_s_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_skryl_cond_vect end interface @@ -180,34 +184,34 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,& select case(psb_toupper(method)) case('CG') call psb_scg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond) + &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2) case('FCG') call psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond) + &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2) case('GCR') call psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('CGS') call psb_scgs_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('BICG') call psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('BICGSTAB') call psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('RGMRES','GMRES') call psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop) + &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2) case('BICGSTABL') call psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop) + &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2) case default if (me == 0) write(psb_err_unit,*) trim(name),& & ': Warning: Unknown method ',method,& & ', defaulting to BiCGSTAB' call psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) end select if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) diff --git a/linsolve/impl/psb_srgmres.f90 b/linsolve/impl/psb_srgmres.f90 index cf455385e..5f90f9a43 100644 --- a/linsolve/impl/psb_srgmres.f90 +++ b/linsolve/impl/psb_srgmres.f90 @@ -97,17 +97,20 @@ ! iterations ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/(|a||x|+|b|); here the iteration is +! 1: err = |r|/(|a||x|+|b|); here +! the iteration is ! stopped when |r| <= eps * (|a||x|+|b|) ! 2: err = |r|/|b|; here the iteration is ! stopped when |r| <= eps * |b| +! 3: Same as 2 but with X and B scaled +! by s1 and s2 ! where r is the (preconditioned, recursive ! estimate of) residual. ! irst - integer(optional) Input: restart parameter ! subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop) + & itmax,iter,err,itrace,irst,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_s_linsolve_conv_mod @@ -123,6 +126,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err + type(psb_s_vect_type), intent(inout), optional :: s1, s2 ! = local data real(psb_spk_), allocatable :: aux(:) real(psb_spk_), allocatable :: c(:), s(:), h(:,:), rs(:), rst(:) @@ -172,19 +176,30 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_,psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then litmax = itmax @@ -253,12 +268,15 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& & w%get_nrows(),w1%get_nrows() - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ani = psb_spnrmi(a,desc_a,info) bni = psb_geamax(b,desc_a,info) - else if (istop_ == 2) then - bn2 = psb_genrm2(b,desc_a,info) - else if (istop_ == 3) then + case(psb_istop_bn2_) + bn2 = psb_genrm2(b,desc_a,info) + case(psb_istop_rn2_abs_) + ! do nothing + case(psb_istop_rrn2_) call psb_geaxpby(sone,b,szero,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -273,7 +291,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if r0n2 = psb_genrm2(v(1),desc_a,info) - endif + end select + errnum = szero errden = sone deps = eps @@ -307,7 +326,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 end if - + if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info) + rs(1) = psb_genrm2(v(1),desc_a,info) rs(2:) = szero if (info /= psb_success_) then @@ -324,20 +344,25 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& ! ! check convergence ! - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) rni = psb_geamax(v(1),desc_a,info) xni = psb_geamax(x,desc_a,info) errnum = rni errden = (ani*xni+bni) - else if (istop_ == 2) then + case(psb_istop_bn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 - else if (istop_ == 3) then + case(psb_istop_rn2_abs_) + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = sone + case(psb_istop_rrn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = r0n2 - endif + end select if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) @@ -381,7 +406,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& h(i+1,i) = szero call srot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_spk_),s(i)) - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ! ! build x and then compute the residual and its infinity norm ! @@ -403,8 +429,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& errnum = rni errden = (ani*xni+bni) ! - - else if (istop_ == 2) then + case(psb_istop_bn2_) ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -412,7 +437,14 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = bn2 - else if (istop_ == 3) then + + case(psb_istop_rn2_abs_) + rni = abs(rs(i+1)) + errnum = rni + errden = sone + + case(psb_istop_rrn2_) + ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -420,14 +452,15 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = r0n2 - endif + end select if (errnum <= eps*errden) then - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(sone,xt,szero,x,desc_a,info) ! = x = xt - else if (istop_ == 2) then + case(psb_istop_bn2_, psb_istop_rn2_abs_) ! ! build x ! @@ -441,7 +474,22 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& end do call prec%apply(w1,w,desc_a,info) call psb_geaxpby(sone,w,sone,x,desc_a,info) - end if + + case(psb_istop_rrn2_) + ! + ! build x + ! + call strsm('l','u','n','n',i,1,sone,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:i) + call w1%set(szero) + do k=1, i + call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(sone,w,sone,x,desc_a,info) + end select if (itrace_ > 0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -454,9 +502,11 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& end do inner - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(sone,xt,szero,x,desc_a,info)! x = xt - else if (istop_ == 2) then + + case(psb_istop_bn2_, psb_istop_rn2_abs_) ! ! build x ! @@ -470,7 +520,22 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& end do call prec%apply(w1,w,desc_a,info) call psb_geaxpby(sone,w,sone,x,desc_a,info) - end if + + case(psb_istop_rrn2_) + ! + ! build x + ! + call strsm('l','u','n','n',nl,1,sone,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:nl) + call w1%set(szero) + do k=1, nl + call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(sone,w,sone,x,desc_a,info) + end select if (itx >= litmax) then if (itrace_ > 0) then diff --git a/linsolve/impl/psb_srichardson.f90 b/linsolve/impl/psb_srichardson.f90 index c42390839..73422792d 100644 --- a/linsolve/impl/psb_srichardson.f90 +++ b/linsolve/impl/psb_srichardson.f90 @@ -120,8 +120,17 @@ Subroutine psb_srichardson_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + + if (present(itmax)) then itmax_ = itmax else diff --git a/linsolve/impl/psb_zbicg.f90 b/linsolve/impl/psb_zbicg.f90 index 7f9829599..55e482423 100644 --- a/linsolve/impl/psb_zbicg.f90 +++ b/linsolve/impl/psb_zbicg.f90 @@ -95,7 +95,7 @@ ! subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop) + & itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_z_linsolve_conv_mod @@ -111,6 +111,7 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), optional, intent(in) :: itmax, itrace, istop integer(psb_ipk_), optional, intent(out) :: iter real(psb_dpk_), optional, intent(out) :: err + type(psb_z_vect_type), intent(inout), optional :: s1, s2 ! !$ local data complex(psb_dpk_), allocatable, target :: aux(:) type(psb_z_vect_type), allocatable, target :: wwrk(:) @@ -160,19 +161,29 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if ! ! istop_ = 1: normwise backward error, infinity norm - ! istop_ = 2: ||r||/||b|| norm 2 + ! istop_ = 2: ||r||/||b|| norm 2 ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then - info=psb_err_invalid_istop_ + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ err=info - call psb_errpush(info,name,i_err=(/istop_/)) + call psb_errpush(info,name,a_err="invalid istop_") goto 9999 - endif + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if(info /= psb_success_) then @@ -226,7 +237,8 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -252,7 +264,7 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,& rho = zzero ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -306,7 +318,7 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(-alpha,q,zone,r,desc_a,info) call psb_geaxpby(-alpha,qt,zone,rt,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_zcg.F90 b/linsolve/impl/psb_zcg.F90 index ceb5948e5..cb34f85d1 100644 --- a/linsolve/impl/psb_zcg.F90 +++ b/linsolve/impl/psb_zcg.F90 @@ -96,7 +96,7 @@ ! ! subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop,cond) + & itmax,iter,err,itrace,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod use psb_z_linsolve_conv_mod @@ -112,6 +112,8 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err,cond + type(psb_z_vect_type), intent(inout), optional :: s1, s2 + ! = Local data complex(psb_dpk_), allocatable, target :: aux(:),td(:),tu(:),eig(:),ewrk(:) integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:) @@ -159,8 +161,29 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& @@ -224,7 +247,8 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,& rho = zzero - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + &desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -268,7 +292,7 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(alpha,p,zone,x,desc_a,info) call psb_geaxpby(-alpha,q,zone,r,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_zcgs.f90 b/linsolve/impl/psb_zcgs.f90 index b90af9603..208976112 100644 --- a/linsolve/impl/psb_zcgs.f90 +++ b/linsolve/impl/psb_zcgs.f90 @@ -93,7 +93,7 @@ ! estimate of) residual. ! Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop) + & itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_z_linsolve_conv_mod @@ -109,6 +109,7 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err + type(psb_z_vect_type), intent(inout), optional :: s1, s2 ! = local data complex(psb_dpk_), allocatable, target :: aux(:) type(psb_z_vect_type), allocatable, target :: wwrk(:) @@ -154,8 +155,29 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,& If (Present(istop)) Then istop_ = istop Else - istop_ = 2 + istop_ = psb_get_istop_default() Endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) @@ -202,7 +224,8 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -225,7 +248,7 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,& ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -299,7 +322,7 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_zcgstab.f90 b/linsolve/impl/psb_zcgstab.f90 index 1784a418f..5777fa209 100644 --- a/linsolve/impl/psb_zcgstab.f90 +++ b/linsolve/impl/psb_zcgstab.f90 @@ -93,7 +93,7 @@ ! where r is the (preconditioned, recursive ! estimate of) residual. ! -Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) +Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_z_linsolve_conv_mod @@ -109,6 +109,7 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err + type(psb_z_vect_type), intent(inout), optional :: s1, s2 ! = Local data complex(psb_dpk_), allocatable, target :: aux(:),wwrk(:,:) type(psb_z_vect_type) :: q, r, p, v, s, t, z, f @@ -156,13 +157,31 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist If (Present(istop)) Then istop_ = istop - Else - istop_ = 2 - Endif + else + istop_ = psb_get_istop_default() + endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if ! - ! ISTOP_ = 1: Normwise backward error, infinity norm - ! ISTOP_ = 2: ||r||/||b|| norm 2 + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select + ! = if (.not.same_type_as(x,b)) then ! = write(0,*) 'Warning: different dynamic types for X and B ' ! = end if @@ -217,7 +236,8 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist End If itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -234,7 +254,7 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist call psb_geaxpby(zone,r,zzero,q,desc_a,info) ! Perhaps we already satisfy the convergence criterion... - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (psb_errstatus_fatal()) then info=psb_err_from_subroutine_ @@ -354,7 +374,7 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist call psb_geaxpby(omega,z,zone,x,desc_a,info) call psb_geaxpby(zone,s,zzero,r,desc_a,info) call psb_geaxpby(-omega,t,zone,r,desc_a,info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ') diff --git a/linsolve/impl/psb_zcgstabl.f90 b/linsolve/impl/psb_zcgstabl.f90 index 679d02ce6..c99cf7d20 100644 --- a/linsolve/impl/psb_zcgstabl.f90 +++ b/linsolve/impl/psb_zcgstabl.f90 @@ -104,7 +104,7 @@ ! ! Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop) + & itmax,iter,err,itrace,irst,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_z_linsolve_conv_mod @@ -120,6 +120,7 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err + type(psb_z_vect_type), intent(inout), optional :: s1, s2 ! = local data complex(psb_dpk_), allocatable, target :: aux(:), gamma(:),& & gamma1(:), gamma2(:), taum(:,:), sigma(:) @@ -172,8 +173,29 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then itmax_ = itmax @@ -246,7 +268,8 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,& rt0 => wwrk(10) - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -284,7 +307,7 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,& & write(debug_unit,*) me,' ',trim(name),& & ' on entry to amax: b: ',b%get_nrows() - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 @@ -388,7 +411,7 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(-gamma1(j),rh(j),zone,rh(0),desc_a,info) enddo - if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/linsolve/impl/psb_zfcg.F90 b/linsolve/impl/psb_zfcg.F90 index 550740ff2..7f2806129 100644 --- a/linsolve/impl/psb_zfcg.F90 +++ b/linsolve/impl/psb_zfcg.F90 @@ -104,7 +104,7 @@ ! ! subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop,cond) + & itmax,iter,err,itrace,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod use psb_z_linsolve_conv_mod @@ -120,6 +120,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_dpk_), Optional, Intent(out) :: err,cond + type(psb_z_vect_type), intent(inout), optional :: s1, s2 ! = Local data type(psb_z_vect_type) :: v, w, d , q, r complex(psb_dpk_) :: alpha, beta, delta, gamma, theta @@ -164,9 +165,29 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif - + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& @@ -207,7 +228,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,& & scratch=.true.,mold=x%v) call psb_init_conv(methdname,istop_,itrace_,itmax_,& - & a,x,b,eps,desc_a,stopdat,info) + & a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2) itx = 0 restart: do @@ -226,7 +247,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,& end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from restart' exit restart end if @@ -282,7 +303,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,& itx = itx + 1 - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from iteration' exit restart end if diff --git a/linsolve/impl/psb_zgcr.f90 b/linsolve/impl/psb_zgcr.f90 index 2566b9e36..e2c64d91a 100644 --- a/linsolve/impl/psb_zgcr.f90 +++ b/linsolve/impl/psb_zgcr.f90 @@ -106,7 +106,7 @@ ! subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace, irst, istop) + & itmax,iter,err,itrace, irst, istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_z_linsolve_conv_mod @@ -124,6 +124,7 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_dpk_), Optional, Intent(out) :: err + type(psb_z_vect_type), intent(inout), optional :: s1, s2 ! = local data complex(psb_dpk_), allocatable :: alpha(:), h(:,:) type(psb_z_vect_type), allocatable :: z(:), c(:), c_scale(:) @@ -167,22 +168,30 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif - - ! - ! ISTOP_ = 1: Normwise backward error, infinity norm - ! ISTOP_ = 2: ||r||/||b||, 2-norm - ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif - - + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_, psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select + call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if (info == psb_success_)& & call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) @@ -245,7 +254,8 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 nrst = -1 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,& + & desc_a,stopdat,info,s1=s1,s2=s2) restart: do if (itx>= itmax_) exit restart h = zzero @@ -268,7 +278,7 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart nrst = nrst + 1 @@ -299,7 +309,7 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(zone, r, zzero, r, desc_a, info) call psb_geaxpby(-alpha(j), c_scale(j), zone, r, desc_a, info) - if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart if (j >= irst) exit iteration diff --git a/linsolve/impl/psb_zkrylov.f90 b/linsolve/impl/psb_zkrylov.f90 index 091271d15..acc41a240 100644 --- a/linsolve/impl/psb_zkrylov.f90 +++ b/linsolve/impl/psb_zkrylov.f90 @@ -80,7 +80,7 @@ ! estimate of) residual ! Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop,cond) + & itmax,iter,err,itrace,irst,istop,cond,s1,s2) use psb_base_mod use psb_prec_mod,only : psb_zprec_type @@ -97,11 +97,12 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err,cond + type(psb_z_vect_type), intent(inout), optional :: s1, s2 abstract interface subroutine psb_zkryl_vect(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) + & desc_a,info,itmax,iter,err,itrace,istop,s1,s2) import :: psb_ipk_, psb_dpk_, psb_desc_type, & & psb_zspmat_type, psb_zprec_type, psb_z_vect_type type(psb_zspmat_type), intent(in) :: a @@ -114,9 +115,10 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), optional, intent(in) :: itmax, itrace,istop integer(psb_ipk_), optional, intent(out) :: iter real(psb_dpk_), optional, intent(out) :: err + type(psb_z_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_zkryl_vect Subroutine psb_zkryl_rest_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,irst,istop) + &itmax,iter,err, itrace,irst,istop,s1,s2) import :: psb_ipk_, psb_dpk_, psb_desc_type, & & psb_zspmat_type, psb_zprec_type, psb_z_vect_type Type(psb_zspmat_type), Intent(in) :: a @@ -129,9 +131,10 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err + type(psb_z_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_zkryl_rest_vect Subroutine psb_zkryl_cond_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,istop,cond) + &itmax,iter,err, itrace,istop,cond,s1,s2) import :: psb_ipk_, psb_dpk_, psb_desc_type, & & psb_zspmat_type, psb_zprec_type, psb_z_vect_type Type(psb_zspmat_type), Intent(in) :: a @@ -144,6 +147,7 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err, cond + type(psb_z_vect_type), intent(inout), optional :: s1, s2 end subroutine psb_zkryl_cond_vect end interface @@ -180,34 +184,34 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& select case(psb_toupper(method)) case('CG') call psb_zcg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond) + &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2) case('FCG') call psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond) + &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2) case('GCR') call psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('CGS') call psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('BICG') call psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('BICGSTAB') call psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) case('RGMRES','GMRES') call psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop) + &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2) case('BICGSTABL') call psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop) + &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2) case default if (me == 0) write(psb_err_unit,*) trim(name),& & ': Warning: Unknown method ',method,& & ', defaulting to BiCGSTAB' call psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace=itrace_,istop=istop) + &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2) end select if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) diff --git a/linsolve/impl/psb_zrgmres.f90 b/linsolve/impl/psb_zrgmres.f90 index 1d2aa3484..e56626223 100644 --- a/linsolve/impl/psb_zrgmres.f90 +++ b/linsolve/impl/psb_zrgmres.f90 @@ -97,17 +97,20 @@ ! iterations ! istop - integer(optional) Input: stopping criterion, or how ! to estimate the error. -! 1: err = |r|/(|a||x|+|b|); here the iteration is +! 1: err = |r|/(|a||x|+|b|); here +! the iteration is ! stopped when |r| <= eps * (|a||x|+|b|) ! 2: err = |r|/|b|; here the iteration is ! stopped when |r| <= eps * |b| +! 3: Same as 2 but with X and B scaled +! by s1 and s2 ! where r is the (preconditioned, recursive ! estimate of) residual. ! irst - integer(optional) Input: restart parameter ! subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop) + & itmax,iter,err,itrace,irst,istop,s1,s2) use psb_base_mod use psb_prec_mod use psb_z_linsolve_conv_mod @@ -123,6 +126,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop integer(psb_ipk_), Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err + type(psb_z_vect_type), intent(inout), optional :: s1, s2 ! = local data complex(psb_dpk_), allocatable :: aux(:) complex(psb_dpk_), allocatable :: c(:), s(:), h(:,:), rs(:), rst(:) @@ -172,19 +176,30 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_,psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then litmax = itmax @@ -253,12 +268,15 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& & w%get_nrows(),w1%get_nrows() - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ani = psb_spnrmi(a,desc_a,info) bni = psb_geamax(b,desc_a,info) - else if (istop_ == 2) then - bn2 = psb_genrm2(b,desc_a,info) - else if (istop_ == 3) then + case(psb_istop_bn2_) + bn2 = psb_genrm2(b,desc_a,info) + case(psb_istop_rn2_abs_) + ! do nothing + case(psb_istop_rrn2_) call psb_geaxpby(zone,b,zzero,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -273,7 +291,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if r0n2 = psb_genrm2(v(1),desc_a,info) - endif + end select + errnum = zzero errden = zone deps = eps @@ -307,7 +326,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 end if - + if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info) + rs(1) = psb_genrm2(v(1),desc_a,info) rs(2:) = zzero if (info /= psb_success_) then @@ -324,20 +344,25 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& ! ! check convergence ! - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) rni = psb_geamax(v(1),desc_a,info) xni = psb_geamax(x,desc_a,info) errnum = rni errden = (ani*xni+bni) - else if (istop_ == 2) then + case(psb_istop_bn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 - else if (istop_ == 3) then + case(psb_istop_rn2_abs_) + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = done + case(psb_istop_rrn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = r0n2 - endif + end select if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) @@ -381,7 +406,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& h(i+1,i) = zzero call zrot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_dpk_),s(i)) - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ! ! build x and then compute the residual and its infinity norm ! @@ -403,8 +429,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& errnum = rni errden = (ani*xni+bni) ! - - else if (istop_ == 2) then + case(psb_istop_bn2_) ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -412,7 +437,14 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = bn2 - else if (istop_ == 3) then + + case(psb_istop_rn2_abs_) + rni = abs(rs(i+1)) + errnum = rni + errden = done + + case(psb_istop_rrn2_) + ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -420,14 +452,15 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = r0n2 - endif + end select if (errnum <= eps*errden) then - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(zone,xt,zzero,x,desc_a,info) ! = x = xt - else if (istop_ == 2) then + case(psb_istop_bn2_, psb_istop_rn2_abs_) ! ! build x ! @@ -441,7 +474,22 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& end do call prec%apply(w1,w,desc_a,info) call psb_geaxpby(zone,w,zone,x,desc_a,info) - end if + + case(psb_istop_rrn2_) + ! + ! build x + ! + call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:i) + call w1%set(zzero) + do k=1, i + call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(zone,w,zone,x,desc_a,info) + end select if (itrace_ > 0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -454,9 +502,11 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& end do inner - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(zone,xt,zzero,x,desc_a,info)! x = xt - else if (istop_ == 2) then + + case(psb_istop_bn2_, psb_istop_rn2_abs_) ! ! build x ! @@ -470,7 +520,22 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& end do call prec%apply(w1,w,desc_a,info) call psb_geaxpby(zone,w,zone,x,desc_a,info) - end if + + case(psb_istop_rrn2_) + ! + ! build x + ! + call ztrsm('l','u','n','n',nl,1,zone,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:nl) + call w1%set(zzero) + do k=1, nl + call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(zone,w,zone,x,desc_a,info) + end select if (itx >= litmax) then if (itrace_ > 0) then diff --git a/linsolve/impl/psb_zrichardson.f90 b/linsolve/impl/psb_zrichardson.f90 index 54c019a51..d03bb95af 100644 --- a/linsolve/impl/psb_zrichardson.f90 +++ b/linsolve/impl/psb_zrichardson.f90 @@ -120,8 +120,17 @@ Subroutine psb_zrichardson_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif + + if (.not.psb_is_valid_istop(istop_)) then + info=psb_err_invalid_istop_ + err=info + call psb_errpush(info,name,i_err=(/istop_/)) + goto 9999 + end if + + if (present(itmax)) then itmax_ = itmax else diff --git a/linsolve/psb_base_linsolve_conv_mod.f90 b/linsolve/psb_base_linsolve_conv_mod.f90 index 8e45614f4..0e45d11e9 100644 --- a/linsolve/psb_base_linsolve_conv_mod.f90 +++ b/linsolve/psb_base_linsolve_conv_mod.f90 @@ -41,10 +41,20 @@ Module psb_base_linsolve_conv_mod module procedure psb_d_end_conv end interface + ! + integer(psb_ipk_), parameter :: psb_istop_min_ = 1 + integer(psb_ipk_), parameter :: psb_istop_ani_ = 1 + integer(psb_ipk_), parameter :: psb_istop_bn2_ = 2 + integer(psb_ipk_), parameter :: psb_istop_rn2_abs_ = 3 + integer(psb_ipk_), parameter :: psb_istop_rrn2_ = 4 + integer(psb_ipk_), parameter :: psb_istop_scbn2_ = 5 + integer(psb_ipk_), parameter :: psb_istop_max_ = 3 + + ! Fields in controls and values integer(psb_ipk_), parameter :: psb_ik_bni_=1, psb_ik_rni_=2, psb_ik_ani_=3 integer(psb_ipk_), parameter :: psb_ik_xni_=4, psb_ik_bn2_=5, psb_ik_r0n2_=6 integer(psb_ipk_), parameter :: psb_ik_xn2_=7, psb_ik_errnum_=8, psb_ik_errden_=9 - integer(psb_ipk_), parameter :: psb_ik_eps_=10, psb_ik_rn2_=11 + integer(psb_ipk_), parameter :: psb_ik_eps_=10, psb_ik_rn2_=11, psb_ik_rn2_abs_=12 integer(psb_ipk_), parameter :: psb_ik_stopc_=1, psb_ik_trace_=2, psb_ik_itmax_=3 integer(psb_ipk_), parameter :: psb_ik_ivsz_=16 type psb_itconv_type @@ -52,8 +62,29 @@ Module psb_base_linsolve_conv_mod real(psb_dpk_) :: values(psb_ik_ivsz_) end type psb_itconv_type + integer(psb_ipk_), save :: psb_istop_default = psb_istop_bn2_ + contains + function psb_is_valid_istop(istop) result(res) + integer(psb_ipk_) :: istop + logical :: res + + res = ((psb_istop_min_<=istop).and.(istop<=psb_istop_max_)) + end function psb_is_valid_istop + + function psb_get_istop_default() result(res) + integer(psb_ipk_) :: res + + res = psb_istop_default + end function psb_get_istop_default + + subroutine psb_set_istop_default(val) + integer(psb_ipk_) :: val + if ((psb_istop_min_<=val).and.(val<=psb_istop_max_)) & + & psb_istop_default = val + end subroutine psb_set_istop_default + subroutine log_header(methdname) !use psb_base_mod implicit none diff --git a/linsolve/psb_c_linsolve_conv_mod.f90 b/linsolve/psb_c_linsolve_conv_mod.f90 index 474acfaa6..840410813 100644 --- a/linsolve/psb_c_linsolve_conv_mod.f90 +++ b/linsolve/psb_c_linsolve_conv_mod.f90 @@ -83,15 +83,18 @@ contains stopdat%controls(psb_ik_itmax_) = itmax select case(stopdat%controls(psb_ik_stopc_)) - case (1) + case (psb_istop_ani_) stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) if (info == psb_success_)& & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) - case (2) + case (psb_istop_bn2_) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) - - case (3) + + case (psb_istop_rn2_abs_) + ! Do nothing + + case (psb_istop_rrn2_) call psb_geall(r,desc_a,info) call psb_geaxpby(cone,b,czero,r,desc_a,info) call psb_spmm(-cone,a,x,cone,r,desc_a,info) @@ -108,8 +111,8 @@ contains end if stopdat%values(psb_ik_eps_) = eps - stopdat%values(psb_ik_errnum_) = dzero - stopdat%values(psb_ik_errden_) = done + stopdat%values(psb_ik_errnum_) = szero + stopdat%values(psb_ik_errden_) = sone if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))& & call log_header(methdname) @@ -123,7 +126,6 @@ contains end subroutine psb_c_init_conv - function psb_c_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res) use psb_base_mod implicit none @@ -149,19 +151,26 @@ contains res = .false. select case(stopdat%controls(psb_ik_stopc_)) - case(1) + case(psb_istop_ani_) stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + if (info == psb_success_) & + & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) stopdat%values(psb_ik_errden_) =& & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)& & +stopdat%values(psb_ik_bni_)) - case(2) + + case(psb_istop_bn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) - case(3) + case (psb_istop_rn2_abs_) + stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_) + stopdat%values(psb_ik_errden_) = sone + + case(psb_istop_rrn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_) @@ -201,8 +210,8 @@ contains end function psb_c_check_conv - - subroutine psb_c_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info) + subroutine psb_c_init_conv_vect(methdname,stopc,trace,itmax,& + & a,x,b,eps,desc_a,stopdat,info,s1,s2) use psb_base_mod implicit none character(len=*), intent(in) :: methdname @@ -213,6 +222,7 @@ contains type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat integer(psb_ipk_), intent(out) :: info + type(psb_c_vect_type), optional :: s1, s2 type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me, np, err_act @@ -236,15 +246,18 @@ contains stopdat%controls(psb_ik_itmax_) = itmax select case(stopdat%controls(psb_ik_stopc_)) - case (1) + case (psb_istop_ani_) stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) if (info == psb_success_)& & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) - case (2) + case (psb_istop_bn2_) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + + case (psb_istop_rn2_abs_) + ! Do nothing - case (3) + case (psb_istop_rrn2_) call psb_geasb(r,desc_a,info,scratch=.true.) call psb_geaxpby(cone,b,czero,r,desc_a,info) call psb_spmm(-cone,a,x,cone,r,desc_a,info) @@ -261,8 +274,8 @@ contains end if stopdat%values(psb_ik_eps_) = eps - stopdat%values(psb_ik_errnum_) = dzero - stopdat%values(psb_ik_errden_) = done + stopdat%values(psb_ik_errnum_) = szero + stopdat%values(psb_ik_errden_) = sone if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))& & call log_header(methdname) @@ -276,7 +289,8 @@ contains end subroutine psb_c_init_conv_vect - function psb_c_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res) + function psb_c_check_conv_vect(methdname,it,x,r,& + & desc_a,stopdat,info,s1,s2) result(res) use psb_base_mod implicit none character(len=*), intent(in) :: methdname @@ -286,6 +300,7 @@ contains type(psb_itconv_type) :: stopdat logical :: res integer(psb_ipk_), intent(out) :: info + type(psb_c_vect_type), optional :: s1, s2 type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me, np, err_act @@ -303,19 +318,26 @@ contains select case(stopdat%controls(psb_ik_stopc_)) - case(1) + case(psb_istop_ani_) stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + if (info == psb_success_) & + & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) stopdat%values(psb_ik_errden_) = & & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)& & +stopdat%values(psb_ik_bni_)) - case(2) + + case(psb_istop_bn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) - case(3) + case (psb_istop_rn2_abs_) + stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_) + stopdat%values(psb_ik_errden_) = sone + + case(psb_istop_rrn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_) diff --git a/linsolve/psb_d_linsolve_conv_mod.f90 b/linsolve/psb_d_linsolve_conv_mod.f90 index f92edce39..64e5c3010 100644 --- a/linsolve/psb_d_linsolve_conv_mod.f90 +++ b/linsolve/psb_d_linsolve_conv_mod.f90 @@ -83,15 +83,18 @@ contains stopdat%controls(psb_ik_itmax_) = itmax select case(stopdat%controls(psb_ik_stopc_)) - case (1) + case (psb_istop_ani_) stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) if (info == psb_success_)& & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) - case (2) + case (psb_istop_bn2_) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) - - case (3) + + case (psb_istop_rn2_abs_) + ! Do nothing + + case (psb_istop_rrn2_) call psb_geall(r,desc_a,info) call psb_geaxpby(done,b,dzero,r,desc_a,info) call psb_spmm(-done,a,x,done,r,desc_a,info) @@ -123,7 +126,6 @@ contains end subroutine psb_d_init_conv - function psb_d_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res) use psb_base_mod implicit none @@ -149,19 +151,26 @@ contains res = .false. select case(stopdat%controls(psb_ik_stopc_)) - case(1) + case(psb_istop_ani_) stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + if (info == psb_success_) & + & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) stopdat%values(psb_ik_errden_) =& & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)& & +stopdat%values(psb_ik_bni_)) - case(2) + + case(psb_istop_bn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) - case(3) + case (psb_istop_rn2_abs_) + stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_) + stopdat%values(psb_ik_errden_) = done + + case(psb_istop_rrn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_) @@ -201,8 +210,8 @@ contains end function psb_d_check_conv - - subroutine psb_d_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info) + subroutine psb_d_init_conv_vect(methdname,stopc,trace,itmax,& + & a,x,b,eps,desc_a,stopdat,info,s1,s2) use psb_base_mod implicit none character(len=*), intent(in) :: methdname @@ -213,6 +222,7 @@ contains type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat integer(psb_ipk_), intent(out) :: info + type(psb_d_vect_type), optional :: s1, s2 type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me, np, err_act @@ -236,15 +246,18 @@ contains stopdat%controls(psb_ik_itmax_) = itmax select case(stopdat%controls(psb_ik_stopc_)) - case (1) + case (psb_istop_ani_) stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) if (info == psb_success_)& & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) - case (2) + case (psb_istop_bn2_) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + + case (psb_istop_rn2_abs_) + ! Do nothing - case (3) + case (psb_istop_rrn2_) call psb_geasb(r,desc_a,info,scratch=.true.) call psb_geaxpby(done,b,dzero,r,desc_a,info) call psb_spmm(-done,a,x,done,r,desc_a,info) @@ -276,7 +289,8 @@ contains end subroutine psb_d_init_conv_vect - function psb_d_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res) + function psb_d_check_conv_vect(methdname,it,x,r,& + & desc_a,stopdat,info,s1,s2) result(res) use psb_base_mod implicit none character(len=*), intent(in) :: methdname @@ -286,6 +300,7 @@ contains type(psb_itconv_type) :: stopdat logical :: res integer(psb_ipk_), intent(out) :: info + type(psb_d_vect_type), optional :: s1, s2 type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me, np, err_act @@ -303,19 +318,26 @@ contains select case(stopdat%controls(psb_ik_stopc_)) - case(1) + case(psb_istop_ani_) stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + if (info == psb_success_) & + & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) stopdat%values(psb_ik_errden_) = & & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)& & +stopdat%values(psb_ik_bni_)) - case(2) + + case(psb_istop_bn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) - case(3) + case (psb_istop_rn2_abs_) + stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_) + stopdat%values(psb_ik_errden_) = done + + case(psb_istop_rrn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_) diff --git a/linsolve/psb_s_linsolve_conv_mod.f90 b/linsolve/psb_s_linsolve_conv_mod.f90 index 8cac66a41..3d7028052 100644 --- a/linsolve/psb_s_linsolve_conv_mod.f90 +++ b/linsolve/psb_s_linsolve_conv_mod.f90 @@ -83,15 +83,18 @@ contains stopdat%controls(psb_ik_itmax_) = itmax select case(stopdat%controls(psb_ik_stopc_)) - case (1) + case (psb_istop_ani_) stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) if (info == psb_success_)& & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) - case (2) + case (psb_istop_bn2_) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) - - case (3) + + case (psb_istop_rn2_abs_) + ! Do nothing + + case (psb_istop_rrn2_) call psb_geall(r,desc_a,info) call psb_geaxpby(sone,b,szero,r,desc_a,info) call psb_spmm(-sone,a,x,sone,r,desc_a,info) @@ -108,8 +111,8 @@ contains end if stopdat%values(psb_ik_eps_) = eps - stopdat%values(psb_ik_errnum_) = dzero - stopdat%values(psb_ik_errden_) = done + stopdat%values(psb_ik_errnum_) = szero + stopdat%values(psb_ik_errden_) = sone if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))& & call log_header(methdname) @@ -123,7 +126,6 @@ contains end subroutine psb_s_init_conv - function psb_s_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res) use psb_base_mod implicit none @@ -149,19 +151,26 @@ contains res = .false. select case(stopdat%controls(psb_ik_stopc_)) - case(1) + case(psb_istop_ani_) stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + if (info == psb_success_) & + & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) stopdat%values(psb_ik_errden_) =& & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)& & +stopdat%values(psb_ik_bni_)) - case(2) + + case(psb_istop_bn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) - case(3) + case (psb_istop_rn2_abs_) + stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_) + stopdat%values(psb_ik_errden_) = sone + + case(psb_istop_rrn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_) @@ -201,8 +210,8 @@ contains end function psb_s_check_conv - - subroutine psb_s_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info) + subroutine psb_s_init_conv_vect(methdname,stopc,trace,itmax,& + & a,x,b,eps,desc_a,stopdat,info,s1,s2) use psb_base_mod implicit none character(len=*), intent(in) :: methdname @@ -213,6 +222,7 @@ contains type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat integer(psb_ipk_), intent(out) :: info + type(psb_s_vect_type), optional :: s1, s2 type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me, np, err_act @@ -236,15 +246,18 @@ contains stopdat%controls(psb_ik_itmax_) = itmax select case(stopdat%controls(psb_ik_stopc_)) - case (1) + case (psb_istop_ani_) stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) if (info == psb_success_)& & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) - case (2) + case (psb_istop_bn2_) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + + case (psb_istop_rn2_abs_) + ! Do nothing - case (3) + case (psb_istop_rrn2_) call psb_geasb(r,desc_a,info,scratch=.true.) call psb_geaxpby(sone,b,szero,r,desc_a,info) call psb_spmm(-sone,a,x,sone,r,desc_a,info) @@ -261,8 +274,8 @@ contains end if stopdat%values(psb_ik_eps_) = eps - stopdat%values(psb_ik_errnum_) = dzero - stopdat%values(psb_ik_errden_) = done + stopdat%values(psb_ik_errnum_) = szero + stopdat%values(psb_ik_errden_) = sone if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))& & call log_header(methdname) @@ -276,7 +289,8 @@ contains end subroutine psb_s_init_conv_vect - function psb_s_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res) + function psb_s_check_conv_vect(methdname,it,x,r,& + & desc_a,stopdat,info,s1,s2) result(res) use psb_base_mod implicit none character(len=*), intent(in) :: methdname @@ -286,6 +300,7 @@ contains type(psb_itconv_type) :: stopdat logical :: res integer(psb_ipk_), intent(out) :: info + type(psb_s_vect_type), optional :: s1, s2 type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me, np, err_act @@ -303,19 +318,26 @@ contains select case(stopdat%controls(psb_ik_stopc_)) - case(1) + case(psb_istop_ani_) stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + if (info == psb_success_) & + & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) stopdat%values(psb_ik_errden_) = & & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)& & +stopdat%values(psb_ik_bni_)) - case(2) + + case(psb_istop_bn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) - case(3) + case (psb_istop_rn2_abs_) + stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_) + stopdat%values(psb_ik_errden_) = sone + + case(psb_istop_rrn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_) diff --git a/linsolve/psb_z_linsolve_conv_mod.f90 b/linsolve/psb_z_linsolve_conv_mod.f90 index d6082262c..efe60f425 100644 --- a/linsolve/psb_z_linsolve_conv_mod.f90 +++ b/linsolve/psb_z_linsolve_conv_mod.f90 @@ -83,15 +83,18 @@ contains stopdat%controls(psb_ik_itmax_) = itmax select case(stopdat%controls(psb_ik_stopc_)) - case (1) + case (psb_istop_ani_) stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) if (info == psb_success_)& & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) - case (2) + case (psb_istop_bn2_) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) - - case (3) + + case (psb_istop_rn2_abs_) + ! Do nothing + + case (psb_istop_rrn2_) call psb_geall(r,desc_a,info) call psb_geaxpby(zone,b,zzero,r,desc_a,info) call psb_spmm(-zone,a,x,zone,r,desc_a,info) @@ -123,7 +126,6 @@ contains end subroutine psb_z_init_conv - function psb_z_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res) use psb_base_mod implicit none @@ -149,19 +151,26 @@ contains res = .false. select case(stopdat%controls(psb_ik_stopc_)) - case(1) + case(psb_istop_ani_) stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + if (info == psb_success_) & + & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) stopdat%values(psb_ik_errden_) =& & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)& & +stopdat%values(psb_ik_bni_)) - case(2) + + case(psb_istop_bn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) - case(3) + case (psb_istop_rn2_abs_) + stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_) + stopdat%values(psb_ik_errden_) = done + + case(psb_istop_rrn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_) @@ -201,8 +210,8 @@ contains end function psb_z_check_conv - - subroutine psb_z_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info) + subroutine psb_z_init_conv_vect(methdname,stopc,trace,itmax,& + & a,x,b,eps,desc_a,stopdat,info,s1,s2) use psb_base_mod implicit none character(len=*), intent(in) :: methdname @@ -213,6 +222,7 @@ contains type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat integer(psb_ipk_), intent(out) :: info + type(psb_z_vect_type), optional :: s1, s2 type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me, np, err_act @@ -236,15 +246,18 @@ contains stopdat%controls(psb_ik_itmax_) = itmax select case(stopdat%controls(psb_ik_stopc_)) - case (1) + case (psb_istop_ani_) stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) if (info == psb_success_)& & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) - case (2) + case (psb_istop_bn2_) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + + case (psb_istop_rn2_abs_) + ! Do nothing - case (3) + case (psb_istop_rrn2_) call psb_geasb(r,desc_a,info,scratch=.true.) call psb_geaxpby(zone,b,zzero,r,desc_a,info) call psb_spmm(-zone,a,x,zone,r,desc_a,info) @@ -276,7 +289,8 @@ contains end subroutine psb_z_init_conv_vect - function psb_z_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res) + function psb_z_check_conv_vect(methdname,it,x,r,& + & desc_a,stopdat,info,s1,s2) result(res) use psb_base_mod implicit none character(len=*), intent(in) :: methdname @@ -286,6 +300,7 @@ contains type(psb_itconv_type) :: stopdat logical :: res integer(psb_ipk_), intent(out) :: info + type(psb_z_vect_type), optional :: s1, s2 type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me, np, err_act @@ -303,19 +318,26 @@ contains select case(stopdat%controls(psb_ik_stopc_)) - case(1) + case(psb_istop_ani_) stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + if (info == psb_success_) & + & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) stopdat%values(psb_ik_errden_) = & & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)& & +stopdat%values(psb_ik_bni_)) - case(2) + + case(psb_istop_bn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) - case(3) + case (psb_istop_rn2_abs_) + stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_) + stopdat%values(psb_ik_errden_) = done + + case(psb_istop_rrn2_) stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)