fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-06 22:55:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

psblas:

Browse Source 
krylov/psb_base_krylov_conv_mod.f90
 krylov/psb_c_krylov_conv_mod.f90
 krylov/psb_cbicg.f90
 krylov/psb_ccg.F90
 krylov/psb_ccgs.f90
 krylov/psb_ccgstab.f90
 krylov/psb_ccgstabl.f90
 krylov/psb_cfcg.F90
 krylov/psb_cgcr.f90
 krylov/psb_crgmres.f90
 krylov/psb_d_krylov_conv_mod.f90
 krylov/psb_dbicg.f90
 krylov/psb_dcg.F90
 krylov/psb_dcgs.f90
 krylov/psb_dcgstab.f90
 krylov/psb_dcgstabl.f90
 krylov/psb_dfcg.F90
 krylov/psb_dgcr.f90
 krylov/psb_drgmres.f90
 krylov/psb_s_krylov_conv_mod.f90
 krylov/psb_sbicg.f90
 krylov/psb_scg.F90
 krylov/psb_scgs.f90
 krylov/psb_scgstab.f90
 krylov/psb_scgstabl.f90
 krylov/psb_sfcg.F90
 krylov/psb_sgcr.f90
 krylov/psb_srgmres.f90
 krylov/psb_z_krylov_conv_mod.f90
 krylov/psb_zbicg.f90
 krylov/psb_zcg.F90
 krylov/psb_zcgs.f90
 krylov/psb_zcgstab.f90
 krylov/psb_zcgstabl.f90
 krylov/psb_zfcg.F90
 krylov/psb_zgcr.f90
 krylov/psb_zrgmres.f90

New stopping criterion.
This commit is contained in:
Salvatore Filippone
2017-02-11 19:12:52 +00:00
parent 0cb97675d9
commit 9ab0539881
37 changed files with 259 additions and 54 deletions
Download Patch File Download Diff File Expand all files Collapse all files
krylov/psb_base_krylov_conv_mod.f90
+3 -2
View File
@@ -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
krylov/psb_c_krylov_conv_mod.f90
+29 -4
View File
@@ -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!")
krylov/psb_cbicg.f90
+1 -1
View File
@@ -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
krylov/psb_ccg.F90
+1 -1
View File
@@ -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
krylov/psb_ccgs.f90
+1 -1
View File
@@ -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
krylov/psb_ccgstab.f90
+1 -1
View File
@@ -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
krylov/psb_ccgstabl.f90
+1 -1
View File
@@ -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
krylov/psb_cfcg.F90
+1 -1
View File
@@ -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
krylov/psb_cgcr.f90
+1 -1
View File
@@ -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
krylov/psb_crgmres.f90
+28 -2
View File
@@ -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
krylov/psb_d_krylov_conv_mod.f90
+29 -4
View File
@@ -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!")
krylov/psb_dbicg.f90
+1 -1
View File
@@ -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
krylov/psb_dcg.F90
+1 -1
View File
@@ -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
krylov/psb_dcgs.f90
+1 -1
View File
@@ -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
krylov/psb_dcgstab.f90
+1 -1
View File
@@ -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
krylov/psb_dcgstabl.f90
+1 -1
View File
@@ -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
krylov/psb_dfcg.F90
+1 -1
View File
@@ -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
krylov/psb_dgcr.f90
+1 -1
View File
@@ -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
krylov/psb_drgmres.f90
+28 -2
View File
@@ -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
krylov/psb_s_krylov_conv_mod.f90
+29 -4
View File
@@ -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!")
krylov/psb_sbicg.f90
+1 -1
View File
@@ -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
krylov/psb_scg.F90
+1 -1
View File
@@ -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
krylov/psb_scgs.f90
+1 -1
View File
@@ -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
krylov/psb_scgstab.f90
+1 -1
View File
@@ -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
krylov/psb_scgstabl.f90
+1 -1
View File
@@ -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
krylov/psb_sfcg.F90
+1 -1
View File
@@ -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
krylov/psb_sgcr.f90
+1 -1
View File
@@ -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
krylov/psb_srgmres.f90
+28 -2
View File
@@ -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
krylov/psb_z_krylov_conv_mod.f90
+29 -4
View File
@@ -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!")
krylov/psb_zbicg.f90
+1 -1
View File
@@ -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
krylov/psb_zcg.F90
+1 -1
View File
@@ -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
krylov/psb_zcgs.f90
+1 -1
View File
@@ -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
krylov/psb_zcgstab.f90
+1 -1
View File
@@ -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
krylov/psb_zcgstabl.f90
+1 -1
View File
@@ -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
krylov/psb_zfcg.F90
+1 -1
View File
@@ -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
krylov/psb_zgcr.f90
+1 -1
View File
@@ -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
krylov/psb_zrgmres.f90
+28 -2
View File
@@ -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