fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Linsolve mods

Browse Source
merge-kinsol-maint
sfilippone 1 week ago
parent 50607ad9bf
commit 53cca9a94c
42 changed files with 569 additions and 470 deletions
Show Stats Download Patch File Download Diff File

5
linsolve/Makefile
Unescape Escape View File

@ -9,6 +9,8 @@ MODDIR=../modules
MODOBJS= psb_base_linsolve_conv_mod.o \
psb_s_linsolve_conv_mod.o psb_c_linsolve_conv_mod.o \
psb_d_linsolve_conv_mod.o psb_z_linsolve_conv_mod.o \
psb_s_linsolve_mod.o psb_c_linsolve_mod.o \
psb_d_linsolve_mod.o psb_z_linsolve_mod.o \
psb_linsolve_mod.o
OBJS=$(MODOBJS)
@ -37,6 +39,9 @@ impld: $(OBJS)
psb_s_linsolve_conv_mod.o psb_c_linsolve_conv_mod.o psb_d_linsolve_conv_mod.o psb_z_linsolve_conv_mod.o: psb_base_linsolve_conv_mod.o
psb_linsolve_conv_mod.o: psb_s_linsolve_conv_mod.o psb_c_linsolve_conv_mod.o psb_d_linsolve_conv_mod.o psb_z_linsolve_conv_mod.o
psb_linsolve_mod.o: psb_s_linsolve_mod.o psb_c_linsolve_mod.o psb_d_linsolve_mod.o psb_z_linsolve_mod.o
$(F90OBJS): $(MODOBJS)
$(OBJS): $(MODDIR)/$(PRECMODNAME)$(.mod) $(MODDIR)/$(BASEMODNAME)$(.mod)

10
linsolve/impl/psb_cbicg.f90
Unescape Escape 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(:)
@ -236,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
@ -262,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
@ -316,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

9
linsolve/impl/psb_ccg.F90
Unescape Escape 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(:)
@ -245,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
@ -289,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

10
linsolve/impl/psb_ccgs.f90
Unescape Escape 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(:)
@ -223,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
@ -246,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
@ -320,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

10
linsolve/impl/psb_ccgstab.f90
Unescape Escape 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
@ -235,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
@ -252,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_
@ -372,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 ')

10
linsolve/impl/psb_ccgstabl.f90
Unescape Escape 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(:)
@ -267,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
@ -305,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
@ -409,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

9
linsolve/impl/psb_cfcg.F90
Unescape Escape 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
@ -227,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
@ -246,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
@ -302,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

10
linsolve/impl/psb_cgcr.f90
Unescape Escape 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(:)
@ -253,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
@ -276,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
@ -307,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

30
linsolve/impl/psb_ckrylov.f90
Unescape Escape 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,37 +184,37 @@ 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('MINRES','PMINRES')
call psb_cminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
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)

52
linsolve/impl/psb_crgmres.f90
Unescape Escape View File

@ -102,13 +102,15 @@
! 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
@ -124,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(:)
@ -268,9 +271,20 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(cone,s1,b,czero,v(1),desc_a,info)
bni = psb_geamax(v(1),desc_a,info)
else
bni = psb_geamax(b,desc_a,info)
end if
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(cone,s1,b,czero,v(1),desc_a,info)
bn2 = psb_genrm2(v(1),desc_a,info)
else
bn2 = psb_genrm2(b,desc_a,info)
end if
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
@ -282,6 +296,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
end if
call psb_spmm(-cone,a,x,cone,v(1),desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
@ -323,7 +338,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
@ -378,8 +394,14 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(cone,a,w1,czero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
@ -391,10 +413,11 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
rst = rs
call psb_geaxpby(cone,x,czero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
call psb_geaxpby(cone,b,czero,w1,desc_a,info)
call psb_spmm(-cone,a,xt,cone,w1,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
rni = psb_geamax(w1,desc_a,info)
xni = psb_geamax(xt,desc_a,info)
errnum = rni
@ -432,7 +455,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(cone,xt,czero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
!
end select
@ -452,7 +476,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(cone,xt,czero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) !
end select
if (itx >= itmax_) then
@ -523,11 +547,12 @@ contains
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
complex(psb_spk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_c_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_cprec_type) :: prec
type(psb_c_vect_type), intent(inout), optional :: s2
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
@ -539,12 +564,13 @@ contains
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
call w%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
call psb_geaxpby(rs(k),v(k),cone,w,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(cone,w,cone,x,desc_a,info)
if (present(s2)) call psb_gediv(s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
call psb_geaxpby(cone,w1,cone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_crgmres_vect

10
linsolve/impl/psb_dbicg.f90
Unescape Escape 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(:)
@ -236,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
@ -262,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
@ -316,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

9
linsolve/impl/psb_dcg.F90
Unescape Escape 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(:)
@ -253,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
@ -306,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

10
linsolve/impl/psb_dcgs.f90
Unescape Escape 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(:)
@ -223,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
@ -246,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
@ -320,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

10
linsolve/impl/psb_dcgstab.f90
Unescape Escape 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
@ -235,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
@ -252,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_
@ -372,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 ')

10
linsolve/impl/psb_dcgstabl.f90
Unescape Escape 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(:)
@ -267,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
@ -305,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
@ -409,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

9
linsolve/impl/psb_dfcg.F90
Unescape Escape 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
@ -227,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
@ -246,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
@ -302,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

10
linsolve/impl/psb_dgcr.f90
Unescape Escape 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(:)
@ -253,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
@ -276,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
@ -307,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

30
linsolve/impl/psb_dkrylov.f90
Unescape Escape 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,37 +184,37 @@ 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('MINRES','PMINRES')
call psb_dminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
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)

52
linsolve/impl/psb_drgmres.f90
Unescape Escape View File

@ -102,13 +102,15 @@
! 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
@ -124,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(:)
@ -268,9 +271,20 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(done,s1,b,dzero,v(1),desc_a,info)
bni = psb_geamax(v(1),desc_a,info)
else
bni = psb_geamax(b,desc_a,info)
end if
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(done,s1,b,dzero,v(1),desc_a,info)
bn2 = psb_genrm2(v(1),desc_a,info)
else
bn2 = psb_genrm2(b,desc_a,info)
end if
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
@ -282,6 +296,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
end if
call psb_spmm(-done,a,x,done,v(1),desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
@ -323,7 +338,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
@ -378,8 +394,14 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(done,a,w1,dzero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
@ -391,10 +413,11 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
rst = rs
call psb_geaxpby(done,x,dzero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
call psb_geaxpby(done,b,dzero,w1,desc_a,info)
call psb_spmm(-done,a,xt,done,w1,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
rni = psb_geamax(w1,desc_a,info)
xni = psb_geamax(xt,desc_a,info)
errnum = rni
@ -432,7 +455,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(done,xt,dzero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
!
end select
@ -452,7 +476,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(done,xt,dzero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) !
end select
if (itx >= itmax_) then
@ -523,11 +547,12 @@ contains
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
real(psb_dpk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_d_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_dprec_type) :: prec
type(psb_d_vect_type), intent(inout), optional :: s2
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
@ -539,12 +564,13 @@ contains
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
call w%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
call psb_geaxpby(rs(k),v(k),done,w,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(done,w,done,x,desc_a,info)
if (present(s2)) call psb_gediv(s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
call psb_geaxpby(done,w1,done,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_drgmres_vect

10
linsolve/impl/psb_sbicg.f90
Unescape Escape 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(:)
@ -236,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
@ -262,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
@ -316,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

9
linsolve/impl/psb_scg.F90
Unescape Escape 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(:)
@ -253,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
@ -306,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

10
linsolve/impl/psb_scgs.f90
Unescape Escape 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(:)
@ -223,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
@ -246,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
@ -320,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

10
linsolve/impl/psb_scgstab.f90
Unescape Escape 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
@ -235,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
@ -252,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_
@ -372,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 ')

10
linsolve/impl/psb_scgstabl.f90
Unescape Escape 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(:)
@ -267,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
@ -305,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
@ -409,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

9
linsolve/impl/psb_sfcg.F90
Unescape Escape 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
@ -227,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
@ -246,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
@ -302,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

10
linsolve/impl/psb_sgcr.f90
Unescape Escape 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(:)
@ -253,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
@ -276,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
@ -307,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

30
linsolve/impl/psb_skrylov.f90
Unescape Escape 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,37 +184,37 @@ 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('MINRES','PMINRES')
call psb_sminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
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)

52
linsolve/impl/psb_srgmres.f90
Unescape Escape View File

@ -102,13 +102,15 @@
! 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
@ -124,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(:)
@ -268,9 +271,20 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(sone,s1,b,szero,v(1),desc_a,info)
bni = psb_geamax(v(1),desc_a,info)
else
bni = psb_geamax(b,desc_a,info)
end if
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(sone,s1,b,szero,v(1),desc_a,info)
bn2 = psb_genrm2(v(1),desc_a,info)
else
bn2 = psb_genrm2(b,desc_a,info)
end if
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
@ -282,6 +296,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
end if
call psb_spmm(-sone,a,x,sone,v(1),desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
@ -323,7 +338,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
@ -378,8 +394,14 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(sone,a,w1,szero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
@ -391,10 +413,11 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
rst = rs
call psb_geaxpby(sone,x,szero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
call psb_geaxpby(sone,b,szero,w1,desc_a,info)
call psb_spmm(-sone,a,xt,sone,w1,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
rni = psb_geamax(w1,desc_a,info)
xni = psb_geamax(xt,desc_a,info)
errnum = rni
@ -432,7 +455,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(sone,xt,szero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
!
end select
@ -452,7 +476,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(sone,xt,szero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) !
end select
if (itx >= itmax_) then
@ -523,11 +547,12 @@ contains
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
real(psb_spk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_s_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_sprec_type) :: prec
type(psb_s_vect_type), intent(inout), optional :: s2
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
@ -539,12 +564,13 @@ contains
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
call w%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
call psb_geaxpby(rs(k),v(k),sone,w,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(sone,w,sone,x,desc_a,info)
if (present(s2)) call psb_gediv(s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
call psb_geaxpby(sone,w1,sone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_srgmres_vect

10
linsolve/impl/psb_zbicg.f90
Unescape Escape 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(:)
@ -236,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
@ -262,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
@ -316,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

9
linsolve/impl/psb_zcg.F90
Unescape Escape 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(:)
@ -245,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
@ -289,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

10
linsolve/impl/psb_zcgs.f90
Unescape Escape 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(:)
@ -223,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
@ -246,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
@ -320,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

10
linsolve/impl/psb_zcgstab.f90
Unescape Escape 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
@ -235,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
@ -252,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_
@ -372,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 ')

10
linsolve/impl/psb_zcgstabl.f90
Unescape Escape 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(:)
@ -267,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
@ -305,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
@ -409,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

9
linsolve/impl/psb_zfcg.F90
Unescape Escape 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
@ -227,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
@ -246,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
@ -302,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

10
linsolve/impl/psb_zgcr.f90
Unescape Escape 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(:)
@ -253,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
@ -276,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
@ -307,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

30
linsolve/impl/psb_zkrylov.f90
Unescape Escape 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,37 +184,37 @@ 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('MINRES','PMINRES')
call psb_zminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
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)

52
linsolve/impl/psb_zrgmres.f90
Unescape Escape View File

@ -102,13 +102,15 @@
! 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
@ -124,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(:)
@ -268,9 +271,20 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(zone,s1,b,zzero,v(1),desc_a,info)
bni = psb_geamax(v(1),desc_a,info)
else
bni = psb_geamax(b,desc_a,info)
end if
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(zone,s1,b,zzero,v(1),desc_a,info)
bn2 = psb_genrm2(v(1),desc_a,info)
else
bn2 = psb_genrm2(b,desc_a,info)
end if
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
@ -282,6 +296,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
end if
call psb_spmm(-zone,a,x,zone,v(1),desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
@ -323,7 +338,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
@ -378,8 +394,14 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(zone,a,w1,zzero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
@ -391,10 +413,11 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
rst = rs
call psb_geaxpby(zone,x,zzero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
call psb_geaxpby(zone,b,zzero,w1,desc_a,info)
call psb_spmm(-zone,a,xt,zone,w1,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
rni = psb_geamax(w1,desc_a,info)
xni = psb_geamax(xt,desc_a,info)
errnum = rni
@ -432,7 +455,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(zone,xt,zzero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
!
end select
@ -452,7 +476,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(zone,xt,zzero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) !
end select
if (itx >= itmax_) then
@ -523,11 +547,12 @@ contains
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
complex(psb_dpk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_z_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_zprec_type) :: prec
type(psb_z_vect_type), intent(inout), optional :: s2
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
@ -539,12 +564,13 @@ contains
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
call w%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
call psb_geaxpby(rs(k),v(k),zone,w,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(zone,w,zone,x,desc_a,info)
if (present(s2)) call psb_gediv(s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
call psb_geaxpby(zone,w1,zone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_zrgmres_vect

68
linsolve/psb_c_linsolve_conv_mod.f90
Unescape Escape 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_)

60
linsolve/psb_d_linsolve_conv_mod.f90
Unescape Escape 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_)

178
linsolve/psb_linsolve_mod.f90
Unescape Escape View File

@ -36,179 +36,9 @@
Module psb_linsolve_mod
use psb_const_mod
public
interface psb_krylov
Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_sspmat_type, &
& psb_spk_, psb_s_vect_type
use psb_prec_mod, only : psb_sprec_type
character(len=*) :: method
Type(psb_sspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_sprec_type), intent(inout) :: prec
type(psb_s_vect_type), Intent(inout) :: b
type(psb_s_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: 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
end Subroutine psb_skrylov_vect
Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_cspmat_type, &
& psb_spk_, psb_c_vect_type
use psb_prec_mod, only : psb_cprec_type
character(len=*) :: method
Type(psb_cspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_cprec_type), intent(inout) :: prec
type(psb_c_vect_type), Intent(inout) :: b
type(psb_c_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: 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
end Subroutine psb_ckrylov_vect
Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_dspmat_type, &
& psb_dpk_, psb_d_vect_type
use psb_prec_mod, only : psb_dprec_type
character(len=*) :: method
Type(psb_dspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_dprec_type), intent(inout) :: prec
type(psb_d_vect_type), Intent(inout) :: b
type(psb_d_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: 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
end Subroutine psb_dkrylov_vect
Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_zspmat_type, &
& psb_dpk_, psb_z_vect_type
use psb_prec_mod, only : psb_zprec_type
character(len=*) :: method
Type(psb_zspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_zprec_type), intent(inout) :: prec
type(psb_z_vect_type), Intent(inout) :: b
type(psb_z_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: 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
end Subroutine psb_zkrylov_vect
end interface
interface psb_richardson
Subroutine psb_srichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_sspmat_type, &
& psb_spk_, psb_s_vect_type
use psb_prec_mod, only : psb_sprec_type
Type(psb_sspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_sprec_type), intent(inout) :: prec
type(psb_s_vect_type), Intent(inout) :: b
type(psb_s_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
end Subroutine psb_srichardson_vect
Subroutine psb_crichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_cspmat_type, &
& psb_spk_, psb_c_vect_type
use psb_prec_mod, only : psb_cprec_type
Type(psb_cspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_cprec_type), intent(inout) :: prec
type(psb_c_vect_type), Intent(inout) :: b
type(psb_c_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
end Subroutine psb_crichardson_vect
Subroutine psb_drichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_dspmat_type, &
& psb_dpk_, psb_d_vect_type
use psb_prec_mod, only : psb_dprec_type
Type(psb_dspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_dprec_type), intent(inout) :: prec
type(psb_d_vect_type), Intent(inout) :: b
type(psb_d_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
end Subroutine psb_drichardson_vect
Subroutine psb_zrichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_zspmat_type, &
& psb_dpk_, psb_z_vect_type
use psb_prec_mod, only : psb_zprec_type
Type(psb_zspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_zprec_type), intent(inout) :: prec
type(psb_z_vect_type), Intent(inout) :: b
type(psb_z_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
end Subroutine psb_zrichardson_vect
end interface
use psb_s_linsolve_mod
use psb_d_linsolve_mod
use psb_c_linsolve_mod
use psb_z_linsolve_mod
end module psb_linsolve_mod

68
linsolve/psb_s_linsolve_conv_mod.f90
Unescape Escape 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_)

60
linsolve/psb_z_linsolve_conv_mod.f90
Unescape Escape 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_)

Write Preview
Loading…
Cancel
Save