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

New stop criteria

Browse Source
This commit is contained in:
sfilippone
2026-05-22 14:56:12 +02:00
parent 9e1c7b775e
commit eabdeda1af
45 changed files with 1452 additions and 489 deletions
Download Patch File Download Diff File Expand all files Collapse all files
linsolve/impl/psb_cbicg.f90
+24 -12
View File
@@ -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
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
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_) 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
linsolve/impl/psb_ccg.F90
+28 -4
View File
@@ -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
linsolve/impl/psb_ccgs.f90
+28 -5
View File
@@ -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
linsolve/impl/psb_ccgstab.f90
+29 -9
View File
@@ -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 ')
linsolve/impl/psb_ccgstabl.f90
+28 -5
View File
@@ -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
linsolve/impl/psb_cfcg.F90
+27 -6
View File
@@ -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
linsolve/impl/psb_cgcr.f90
+25 -15
View File
@@ -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
linsolve/impl/psb_ckrylov.f90
+17 -13
View File
@@ -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)
linsolve/impl/psb_crgmres.f90
+96 -31
View File
@@ -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
linsolve/impl/psb_crichardson.f90
+10 -1
View File
@@ -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
linsolve/impl/psb_dbicg.f90
+24 -12
View File
@@ -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
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
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_) 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
linsolve/impl/psb_dcg.F90
+28 -4
View File
@@ -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
linsolve/impl/psb_dcgs.f90
+28 -5
View File
@@ -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
linsolve/impl/psb_dcgstab.f90
+29 -9
View File
@@ -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 ')
linsolve/impl/psb_dcgstabl.f90
+28 -5
View File
@@ -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
linsolve/impl/psb_dfcg.F90
+27 -6
View File
@@ -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
linsolve/impl/psb_dgcr.f90
+25 -15
View File
@@ -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
linsolve/impl/psb_dkrylov.f90
+17 -13
View File
@@ -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)
linsolve/impl/psb_drgmres.f90
+96 -31
View File
@@ -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
linsolve/impl/psb_drichardson.f90
+10 -1
View File
@@ -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
linsolve/impl/psb_sbicg.f90
+24 -12
View File
@@ -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
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
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_) 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
linsolve/impl/psb_scg.F90
+28 -4
View File
@@ -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
linsolve/impl/psb_scgs.f90
+28 -5
View File
@@ -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
linsolve/impl/psb_scgstab.f90
+29 -9
View File
@@ -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 ')
linsolve/impl/psb_scgstabl.f90
+28 -5
View File
@@ -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
linsolve/impl/psb_sfcg.F90
+27 -6
View File
@@ -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
linsolve/impl/psb_sgcr.f90
+25 -15
View File
@@ -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
linsolve/impl/psb_skrylov.f90
+17 -13
View File
@@ -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)
linsolve/impl/psb_srgmres.f90
+96 -31
View File
@@ -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
linsolve/impl/psb_srichardson.f90
+10 -1
View File
@@ -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
linsolve/impl/psb_zbicg.f90
+24 -12
View File
@@ -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
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
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_) 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
linsolve/impl/psb_zcg.F90
+28 -4
View File
@@ -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
linsolve/impl/psb_zcgs.f90
+28 -5
View File
@@ -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
linsolve/impl/psb_zcgstab.f90
+29 -9
View File
@@ -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 ')
linsolve/impl/psb_zcgstabl.f90
+28 -5
View File
@@ -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
linsolve/impl/psb_zfcg.F90
+27 -6
View File
@@ -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
linsolve/impl/psb_zgcr.f90
+25 -15
View File
@@ -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
linsolve/impl/psb_zkrylov.f90
+17 -13
View File
@@ -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)
linsolve/impl/psb_zrgmres.f90
+96 -31
View File
@@ -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
linsolve/impl/psb_zrichardson.f90
+10 -1
View File
@@ -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
linsolve/psb_base_linsolve_conv_mod.f90
+32 -1
View File
@@ -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
linsolve/psb_c_linsolve_conv_mod.f90
+45 -23
View File
@@ -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_)
linsolve/psb_d_linsolve_conv_mod.f90
+41 -19
View File
@@ -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_)
linsolve/psb_s_linsolve_conv_mod.f90
+45 -23
View File
@@ -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_)
linsolve/psb_z_linsolve_conv_mod.f90
+41 -19
View File
@@ -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_)