diff --git a/krylov/psb_base_krylov_conv_mod.f90 b/krylov/psb_base_krylov_conv_mod.f90 index a00b20c8..9eb1a366 100644 --- a/krylov/psb_base_krylov_conv_mod.f90 +++ b/krylov/psb_base_krylov_conv_mod.f90 @@ -42,8 +42,9 @@ Module psb_base_krylov_conv_mod end interface 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_xn2_=6 - integer(psb_ipk_), parameter :: psb_ik_errnum_=7, psb_ik_errden_=8, psb_ik_eps_=9, psb_ik_rn2_=10 + 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_stopc_=1, psb_ik_trace_=2, psb_ik_itmax_=3 integer(psb_ipk_), parameter :: psb_ik_ivsz_=16 type psb_itconv_type diff --git a/krylov/psb_c_krylov_conv_mod.f90 b/krylov/psb_c_krylov_conv_mod.f90 index 632e84b7..fbca1cbb 100644 --- a/krylov/psb_c_krylov_conv_mod.f90 +++ b/krylov/psb_c_krylov_conv_mod.f90 @@ -48,13 +48,14 @@ Module psb_c_krylov_conv_mod contains - subroutine psb_c_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + subroutine psb_c_init_conv(methdname,stopc,trace,itmax,a,x,b,eps,& + & desc_a,stopdat,info) use psb_base_mod implicit none character(len=*), intent(in) :: methdname integer(psb_ipk_), intent(in) :: stopc, trace, itmax type(psb_cspmat_type), intent(in) :: a - complex(psb_spk_), intent(in) :: b(:) + complex(psb_spk_), intent(inout) :: b(:), x(:) real(psb_spk_), intent(in) :: eps type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat @@ -62,6 +63,7 @@ contains integer(psb_ipk_) :: ictxt, me, np, err_act, ierr(5) character(len=20) :: name + complex(psb_spk_), allocatable :: r(:) info = psb_success_ name = 'psb_init_conv' @@ -88,6 +90,12 @@ contains case (2) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + case (3) + 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) + stopdat%values(psb_ik_r0n2_) = psb_genrm2(r,desc_a,info) + call psb_gefree(r,desc_a,info) case default info=psb_err_invalid_istop_ ierr(1) = stopc @@ -152,6 +160,11 @@ contains stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + case(3) + 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_) + case default info=psb_err_internal_error_ call psb_errpush(info,name,a_err="Control data in stopdat messed up!") @@ -188,20 +201,21 @@ contains end function psb_c_check_conv - subroutine psb_c_init_conv_vect(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + subroutine psb_c_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info) use psb_base_mod implicit none character(len=*), intent(in) :: methdname integer(psb_ipk_), intent(in) :: stopc, trace,itmax type(psb_cspmat_type), intent(in) :: a real(psb_spk_), intent(in) :: eps - type(psb_c_vect_type), intent(inout) :: b + type(psb_c_vect_type), intent(inout) :: x, b type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: ictxt, me, np, err_act, ierr(5) character(len=20) :: name + type(psb_c_vect_type) :: r info = psb_success_ name = 'psb_init_conv' @@ -228,6 +242,12 @@ contains case (2) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + case (3) + 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) + stopdat%values(psb_ik_r0n2_) = psb_genrm2(r,desc_a,info) + call psb_gefree(r,desc_a,info) case default info=psb_err_invalid_istop_ ierr(1) = stopc @@ -293,6 +313,11 @@ contains stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + case(3) + 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_) + case default info=psb_err_internal_error_ call psb_errpush(info,name,a_err="Control data in stopdat messed up!") diff --git a/krylov/psb_cbicg.f90 b/krylov/psb_cbicg.f90 index d15674a8..12d91f66 100644 --- a/krylov/psb_cbicg.f90 +++ b/krylov/psb_cbicg.f90 @@ -226,7 +226,7 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_ccg.F90 b/krylov/psb_ccg.F90 index 07017652..c9787f82 100644 --- a/krylov/psb_ccg.F90 +++ b/krylov/psb_ccg.F90 @@ -222,7 +222,7 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,& rho = czero - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_ccgs.f90 b/krylov/psb_ccgs.f90 index 0a52bc5f..1b35ef98 100644 --- a/krylov/psb_ccgs.f90 +++ b/krylov/psb_ccgs.f90 @@ -200,7 +200,7 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_ccgstab.f90 b/krylov/psb_ccgstab.f90 index b5ae0f02..4834cf2e 100644 --- a/krylov/psb_ccgstab.f90 +++ b/krylov/psb_ccgstab.f90 @@ -215,7 +215,7 @@ 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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_ccgstabl.f90 b/krylov/psb_ccgstabl.f90 index 9382e312..ca2bcd63 100644 --- a/krylov/psb_ccgstabl.f90 +++ b/krylov/psb_ccgstabl.f90 @@ -245,7 +245,7 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,& rt0 => wwrk(10) - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_cfcg.F90 b/krylov/psb_cfcg.F90 index b47fb7e1..72ea8a60 100644 --- a/krylov/psb_cfcg.F90 +++ b/krylov/psb_cfcg.F90 @@ -209,7 +209,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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) itx=0 restart: do diff --git a/krylov/psb_cgcr.f90 b/krylov/psb_cgcr.f90 index 76dba4d6..933652fb 100644 --- a/krylov/psb_cgcr.f90 +++ b/krylov/psb_cgcr.f90 @@ -246,7 +246,7 @@ 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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) restart: do if (itx>= itmax_) exit restart h = czero diff --git a/krylov/psb_crgmres.f90 b/krylov/psb_crgmres.f90 index 13875c1d..505f6609 100644 --- a/krylov/psb_crgmres.f90 +++ b/krylov/psb_crgmres.f90 @@ -137,7 +137,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_) :: itx, i, istop_, err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: ictxt, np, me - Real(psb_spk_) :: rni, xni, bni, ani,bn2, dt + Real(psb_spk_) :: rni, xni, bni, ani,bn2, dt, r0n2 real(psb_dpk_) :: errnum, errden, deps, derr character(len=20) :: name character(len=*), parameter :: methdname='RGMRES' @@ -258,6 +258,21 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,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 + call psb_geaxpby(cone,b,czero,v(1),desc_a,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + call psb_spmm(-cone,a,x,cone,v(1),desc_a,info,work=aux) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + r0n2 = psb_genrm2(v(1),desc_a,info) endif errnum = czero errden = cone @@ -274,7 +289,6 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& ! compute r0 = b-ax0 ! check convergence - ! compute v1 = r0/||r0||_2 if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),& @@ -319,6 +333,10 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 + else if (istop_ == 3) then + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = r0n2 endif if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -397,6 +415,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 + ! + ! compute the residual 2-norm as byproduct of the solution + ! procedure of the least-squares problem + ! + rni = abs(rs(i+1)) + errnum = rni + errden = r0n2 endif if (errnum <= eps*errden) then diff --git a/krylov/psb_d_krylov_conv_mod.f90 b/krylov/psb_d_krylov_conv_mod.f90 index 2f839f7a..1a28f0a8 100644 --- a/krylov/psb_d_krylov_conv_mod.f90 +++ b/krylov/psb_d_krylov_conv_mod.f90 @@ -48,13 +48,14 @@ Module psb_d_krylov_conv_mod contains - subroutine psb_d_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + subroutine psb_d_init_conv(methdname,stopc,trace,itmax,a,x,b,eps,& + & desc_a,stopdat,info) use psb_base_mod implicit none character(len=*), intent(in) :: methdname integer(psb_ipk_), intent(in) :: stopc, trace, itmax type(psb_dspmat_type), intent(in) :: a - real(psb_dpk_), intent(in) :: b(:) + real(psb_dpk_), intent(inout) :: b(:), x(:) real(psb_dpk_), intent(in) :: eps type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat @@ -62,6 +63,7 @@ contains integer(psb_ipk_) :: ictxt, me, np, err_act, ierr(5) character(len=20) :: name + real(psb_dpk_), allocatable :: r(:) info = psb_success_ name = 'psb_init_conv' @@ -88,6 +90,12 @@ contains case (2) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + case (3) + 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) + stopdat%values(psb_ik_r0n2_) = psb_genrm2(r,desc_a,info) + call psb_gefree(r,desc_a,info) case default info=psb_err_invalid_istop_ ierr(1) = stopc @@ -152,6 +160,11 @@ contains stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + case(3) + 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_) + case default info=psb_err_internal_error_ call psb_errpush(info,name,a_err="Control data in stopdat messed up!") @@ -188,20 +201,21 @@ contains end function psb_d_check_conv - subroutine psb_d_init_conv_vect(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + subroutine psb_d_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info) use psb_base_mod implicit none character(len=*), intent(in) :: methdname integer(psb_ipk_), intent(in) :: stopc, trace,itmax type(psb_dspmat_type), intent(in) :: a real(psb_dpk_), intent(in) :: eps - type(psb_d_vect_type), intent(inout) :: b + type(psb_d_vect_type), intent(inout) :: x, b type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: ictxt, me, np, err_act, ierr(5) character(len=20) :: name + type(psb_d_vect_type) :: r info = psb_success_ name = 'psb_init_conv' @@ -228,6 +242,12 @@ contains case (2) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + case (3) + 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) + stopdat%values(psb_ik_r0n2_) = psb_genrm2(r,desc_a,info) + call psb_gefree(r,desc_a,info) case default info=psb_err_invalid_istop_ ierr(1) = stopc @@ -293,6 +313,11 @@ contains stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + case(3) + 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_) + case default info=psb_err_internal_error_ call psb_errpush(info,name,a_err="Control data in stopdat messed up!") diff --git a/krylov/psb_dbicg.f90 b/krylov/psb_dbicg.f90 index a347ddf1..55133945 100644 --- a/krylov/psb_dbicg.f90 +++ b/krylov/psb_dbicg.f90 @@ -226,7 +226,7 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_dcg.F90 b/krylov/psb_dcg.F90 index d42ff2bf..152293f6 100644 --- a/krylov/psb_dcg.F90 +++ b/krylov/psb_dcg.F90 @@ -230,7 +230,7 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,& rho = dzero - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_dcgs.f90 b/krylov/psb_dcgs.f90 index 83b8c4a7..675c0881 100644 --- a/krylov/psb_dcgs.f90 +++ b/krylov/psb_dcgs.f90 @@ -200,7 +200,7 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_dcgstab.f90 b/krylov/psb_dcgstab.f90 index 272b2f07..5c879cd7 100644 --- a/krylov/psb_dcgstab.f90 +++ b/krylov/psb_dcgstab.f90 @@ -215,7 +215,7 @@ 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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_dcgstabl.f90 b/krylov/psb_dcgstabl.f90 index 738ef03f..b45811bc 100644 --- a/krylov/psb_dcgstabl.f90 +++ b/krylov/psb_dcgstabl.f90 @@ -245,7 +245,7 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,& rt0 => wwrk(10) - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_dfcg.F90 b/krylov/psb_dfcg.F90 index a6456327..b3d8d0b3 100644 --- a/krylov/psb_dfcg.F90 +++ b/krylov/psb_dfcg.F90 @@ -209,7 +209,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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) itx=0 restart: do diff --git a/krylov/psb_dgcr.f90 b/krylov/psb_dgcr.f90 index 4a6c0e19..3c9c83f5 100644 --- a/krylov/psb_dgcr.f90 +++ b/krylov/psb_dgcr.f90 @@ -246,7 +246,7 @@ 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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) restart: do if (itx>= itmax_) exit restart h = dzero diff --git a/krylov/psb_drgmres.f90 b/krylov/psb_drgmres.f90 index 1aa2b1f2..27a62ed9 100644 --- a/krylov/psb_drgmres.f90 +++ b/krylov/psb_drgmres.f90 @@ -137,7 +137,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_) :: itx, i, istop_, err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: ictxt, np, me - Real(psb_dpk_) :: rni, xni, bni, ani,bn2, dt + Real(psb_dpk_) :: rni, xni, bni, ani,bn2, dt, r0n2 real(psb_dpk_) :: errnum, errden, deps, derr character(len=20) :: name character(len=*), parameter :: methdname='RGMRES' @@ -258,6 +258,21 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,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 + call psb_geaxpby(done,b,dzero,v(1),desc_a,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + call psb_spmm(-done,a,x,done,v(1),desc_a,info,work=aux) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + r0n2 = psb_genrm2(v(1),desc_a,info) endif errnum = dzero errden = done @@ -274,7 +289,6 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& ! compute r0 = b-ax0 ! check convergence - ! compute v1 = r0/||r0||_2 if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),& @@ -319,6 +333,10 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 + else if (istop_ == 3) then + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = r0n2 endif if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -397,6 +415,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 + ! + ! compute the residual 2-norm as byproduct of the solution + ! procedure of the least-squares problem + ! + rni = abs(rs(i+1)) + errnum = rni + errden = r0n2 endif if (errnum <= eps*errden) then diff --git a/krylov/psb_s_krylov_conv_mod.f90 b/krylov/psb_s_krylov_conv_mod.f90 index daff12ed..aac4736f 100644 --- a/krylov/psb_s_krylov_conv_mod.f90 +++ b/krylov/psb_s_krylov_conv_mod.f90 @@ -48,13 +48,14 @@ Module psb_s_krylov_conv_mod contains - subroutine psb_s_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + subroutine psb_s_init_conv(methdname,stopc,trace,itmax,a,x,b,eps,& + & desc_a,stopdat,info) use psb_base_mod implicit none character(len=*), intent(in) :: methdname integer(psb_ipk_), intent(in) :: stopc, trace, itmax type(psb_sspmat_type), intent(in) :: a - real(psb_spk_), intent(in) :: b(:) + real(psb_spk_), intent(inout) :: b(:), x(:) real(psb_spk_), intent(in) :: eps type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat @@ -62,6 +63,7 @@ contains integer(psb_ipk_) :: ictxt, me, np, err_act, ierr(5) character(len=20) :: name + real(psb_spk_), allocatable :: r(:) info = psb_success_ name = 'psb_init_conv' @@ -88,6 +90,12 @@ contains case (2) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + case (3) + 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) + stopdat%values(psb_ik_r0n2_) = psb_genrm2(r,desc_a,info) + call psb_gefree(r,desc_a,info) case default info=psb_err_invalid_istop_ ierr(1) = stopc @@ -152,6 +160,11 @@ contains stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + case(3) + 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_) + case default info=psb_err_internal_error_ call psb_errpush(info,name,a_err="Control data in stopdat messed up!") @@ -188,20 +201,21 @@ contains end function psb_s_check_conv - subroutine psb_s_init_conv_vect(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + subroutine psb_s_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info) use psb_base_mod implicit none character(len=*), intent(in) :: methdname integer(psb_ipk_), intent(in) :: stopc, trace,itmax type(psb_sspmat_type), intent(in) :: a real(psb_spk_), intent(in) :: eps - type(psb_s_vect_type), intent(inout) :: b + type(psb_s_vect_type), intent(inout) :: x, b type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: ictxt, me, np, err_act, ierr(5) character(len=20) :: name + type(psb_s_vect_type) :: r info = psb_success_ name = 'psb_init_conv' @@ -228,6 +242,12 @@ contains case (2) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + case (3) + 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) + stopdat%values(psb_ik_r0n2_) = psb_genrm2(r,desc_a,info) + call psb_gefree(r,desc_a,info) case default info=psb_err_invalid_istop_ ierr(1) = stopc @@ -293,6 +313,11 @@ contains stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + case(3) + 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_) + case default info=psb_err_internal_error_ call psb_errpush(info,name,a_err="Control data in stopdat messed up!") diff --git a/krylov/psb_sbicg.f90 b/krylov/psb_sbicg.f90 index b47b463e..3f6470af 100644 --- a/krylov/psb_sbicg.f90 +++ b/krylov/psb_sbicg.f90 @@ -226,7 +226,7 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_scg.F90 b/krylov/psb_scg.F90 index 74e2b6dd..2bee8573 100644 --- a/krylov/psb_scg.F90 +++ b/krylov/psb_scg.F90 @@ -230,7 +230,7 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,& rho = szero - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_scgs.f90 b/krylov/psb_scgs.f90 index 9a566a44..64f11ee0 100644 --- a/krylov/psb_scgs.f90 +++ b/krylov/psb_scgs.f90 @@ -200,7 +200,7 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_scgstab.f90 b/krylov/psb_scgstab.f90 index d1a4ff81..30c0e4c9 100644 --- a/krylov/psb_scgstab.f90 +++ b/krylov/psb_scgstab.f90 @@ -215,7 +215,7 @@ 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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_scgstabl.f90 b/krylov/psb_scgstabl.f90 index f3d65ea3..99c675f6 100644 --- a/krylov/psb_scgstabl.f90 +++ b/krylov/psb_scgstabl.f90 @@ -245,7 +245,7 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,& rt0 => wwrk(10) - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_sfcg.F90 b/krylov/psb_sfcg.F90 index 78dba287..8109e46f 100644 --- a/krylov/psb_sfcg.F90 +++ b/krylov/psb_sfcg.F90 @@ -209,7 +209,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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) itx=0 restart: do diff --git a/krylov/psb_sgcr.f90 b/krylov/psb_sgcr.f90 index d1a7ec9c..22e3bdad 100644 --- a/krylov/psb_sgcr.f90 +++ b/krylov/psb_sgcr.f90 @@ -246,7 +246,7 @@ 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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) restart: do if (itx>= itmax_) exit restart h = szero diff --git a/krylov/psb_srgmres.f90 b/krylov/psb_srgmres.f90 index 1d30b744..5bfb0cfd 100644 --- a/krylov/psb_srgmres.f90 +++ b/krylov/psb_srgmres.f90 @@ -137,7 +137,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_) :: itx, i, istop_, err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: ictxt, np, me - Real(psb_spk_) :: rni, xni, bni, ani,bn2, dt + Real(psb_spk_) :: rni, xni, bni, ani,bn2, dt, r0n2 real(psb_dpk_) :: errnum, errden, deps, derr character(len=20) :: name character(len=*), parameter :: methdname='RGMRES' @@ -258,6 +258,21 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,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 + call psb_geaxpby(sone,b,szero,v(1),desc_a,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + call psb_spmm(-sone,a,x,sone,v(1),desc_a,info,work=aux) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + r0n2 = psb_genrm2(v(1),desc_a,info) endif errnum = szero errden = sone @@ -274,7 +289,6 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& ! compute r0 = b-ax0 ! check convergence - ! compute v1 = r0/||r0||_2 if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),& @@ -319,6 +333,10 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 + else if (istop_ == 3) then + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = r0n2 endif if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -397,6 +415,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 + ! + ! compute the residual 2-norm as byproduct of the solution + ! procedure of the least-squares problem + ! + rni = abs(rs(i+1)) + errnum = rni + errden = r0n2 endif if (errnum <= eps*errden) then diff --git a/krylov/psb_z_krylov_conv_mod.f90 b/krylov/psb_z_krylov_conv_mod.f90 index f6e49dd5..46915450 100644 --- a/krylov/psb_z_krylov_conv_mod.f90 +++ b/krylov/psb_z_krylov_conv_mod.f90 @@ -48,13 +48,14 @@ Module psb_z_krylov_conv_mod contains - subroutine psb_z_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + subroutine psb_z_init_conv(methdname,stopc,trace,itmax,a,x,b,eps,& + & desc_a,stopdat,info) use psb_base_mod implicit none character(len=*), intent(in) :: methdname integer(psb_ipk_), intent(in) :: stopc, trace, itmax type(psb_zspmat_type), intent(in) :: a - complex(psb_dpk_), intent(in) :: b(:) + complex(psb_dpk_), intent(inout) :: b(:), x(:) real(psb_dpk_), intent(in) :: eps type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat @@ -62,6 +63,7 @@ contains integer(psb_ipk_) :: ictxt, me, np, err_act, ierr(5) character(len=20) :: name + complex(psb_dpk_), allocatable :: r(:) info = psb_success_ name = 'psb_init_conv' @@ -88,6 +90,12 @@ contains case (2) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + case (3) + 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) + stopdat%values(psb_ik_r0n2_) = psb_genrm2(r,desc_a,info) + call psb_gefree(r,desc_a,info) case default info=psb_err_invalid_istop_ ierr(1) = stopc @@ -152,6 +160,11 @@ contains stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + case(3) + 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_) + case default info=psb_err_internal_error_ call psb_errpush(info,name,a_err="Control data in stopdat messed up!") @@ -188,20 +201,21 @@ contains end function psb_z_check_conv - subroutine psb_z_init_conv_vect(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + subroutine psb_z_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info) use psb_base_mod implicit none character(len=*), intent(in) :: methdname integer(psb_ipk_), intent(in) :: stopc, trace,itmax type(psb_zspmat_type), intent(in) :: a real(psb_dpk_), intent(in) :: eps - type(psb_z_vect_type), intent(inout) :: b + type(psb_z_vect_type), intent(inout) :: x, b type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: ictxt, me, np, err_act, ierr(5) character(len=20) :: name + type(psb_z_vect_type) :: r info = psb_success_ name = 'psb_init_conv' @@ -228,6 +242,12 @@ contains case (2) stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + case (3) + 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) + stopdat%values(psb_ik_r0n2_) = psb_genrm2(r,desc_a,info) + call psb_gefree(r,desc_a,info) case default info=psb_err_invalid_istop_ ierr(1) = stopc @@ -293,6 +313,11 @@ contains stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + case(3) + 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_) + case default info=psb_err_internal_error_ call psb_errpush(info,name,a_err="Control data in stopdat messed up!") diff --git a/krylov/psb_zbicg.f90 b/krylov/psb_zbicg.f90 index 57b81bb1..7d17a042 100644 --- a/krylov/psb_zbicg.f90 +++ b/krylov/psb_zbicg.f90 @@ -226,7 +226,7 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_zcg.F90 b/krylov/psb_zcg.F90 index 23e25141..87f4d43a 100644 --- a/krylov/psb_zcg.F90 +++ b/krylov/psb_zcg.F90 @@ -222,7 +222,7 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,& rho = zzero - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_zcgs.f90 b/krylov/psb_zcgs.f90 index ba68ddb3..fe179ef2 100644 --- a/krylov/psb_zcgs.f90 +++ b/krylov/psb_zcgs.f90 @@ -200,7 +200,7 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,& itx = 0 - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_zcgstab.f90 b/krylov/psb_zcgstab.f90 index 599b8429..cedc6cd5 100644 --- a/krylov/psb_zcgstab.f90 +++ b/krylov/psb_zcgstab.f90 @@ -215,7 +215,7 @@ 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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (psb_errstatus_fatal()) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_zcgstabl.f90 b/krylov/psb_zcgstabl.f90 index eca398e7..bd78e9f9 100644 --- a/krylov/psb_zcgstabl.f90 +++ b/krylov/psb_zcgstabl.f90 @@ -245,7 +245,7 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,& rt0 => wwrk(10) - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) if (info /= psb_success_) Then call psb_errpush(psb_err_from_subroutine_non_,name) goto 9999 diff --git a/krylov/psb_zfcg.F90 b/krylov/psb_zfcg.F90 index 8a30c591..4380ee20 100644 --- a/krylov/psb_zfcg.F90 +++ b/krylov/psb_zfcg.F90 @@ -209,7 +209,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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) itx=0 restart: do diff --git a/krylov/psb_zgcr.f90 b/krylov/psb_zgcr.f90 index a19a7b4f..4dbfe3f8 100644 --- a/krylov/psb_zgcr.f90 +++ b/krylov/psb_zgcr.f90 @@ -246,7 +246,7 @@ 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,b,eps,desc_a,stopdat,info) + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) restart: do if (itx>= itmax_) exit restart h = zzero diff --git a/krylov/psb_zrgmres.f90 b/krylov/psb_zrgmres.f90 index 68d25f85..26209935 100644 --- a/krylov/psb_zrgmres.f90 +++ b/krylov/psb_zrgmres.f90 @@ -137,7 +137,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_) :: itx, i, istop_, err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: ictxt, np, me - Real(psb_dpk_) :: rni, xni, bni, ani,bn2, dt + Real(psb_dpk_) :: rni, xni, bni, ani,bn2, dt, r0n2 real(psb_dpk_) :: errnum, errden, deps, derr character(len=20) :: name character(len=*), parameter :: methdname='RGMRES' @@ -258,6 +258,21 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,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 + call psb_geaxpby(zone,b,zzero,v(1),desc_a,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + call psb_spmm(-zone,a,x,zone,v(1),desc_a,info,work=aux) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + r0n2 = psb_genrm2(v(1),desc_a,info) endif errnum = zzero errden = zone @@ -274,7 +289,6 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& ! compute r0 = b-ax0 ! check convergence - ! compute v1 = r0/||r0||_2 if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),& @@ -319,6 +333,10 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 + else if (istop_ == 3) then + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = r0n2 endif if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -397,6 +415,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 + ! + ! compute the residual 2-norm as byproduct of the solution + ! procedure of the least-squares problem + ! + rni = abs(rs(i+1)) + errnum = rni + errden = r0n2 endif if (errnum <= eps*errden) then