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

Change handling of stopping criterion, refactor GMRES

Browse Source
development
sfilippone 12 hours ago
parent 014abc941c
commit d99e759fd1
33 changed files with 1152 additions and 469 deletions
Show Stats Download Patch File Download Diff File

24
linsolve/impl/psb_cbicg.f90
Unescape Escape View File

@ -160,19 +160,29 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
! istop_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=psb_err_invalid_istop_
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
endif
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then

23
linsolve/impl/psb_ccg.F90
Unescape Escape View File

@ -159,8 +159,29 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&

23
linsolve/impl/psb_ccgs.f90
Unescape Escape View File

@ -154,8 +154,29 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
istop_ = psb_get_istop_default()
Endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)

28
linsolve/impl/psb_ccgstab.f90
Unescape Escape View File

@ -156,13 +156,31 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
Endif
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b|| norm 2
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
! = if (.not.same_type_as(x,b)) then
! = write(0,*) 'Warning: different dynamic types for X and B '
! = end if

23
linsolve/impl/psb_ccgstabl.f90
Unescape Escape View File

@ -172,8 +172,29 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
itmax_ = itmax

24
linsolve/impl/psb_cfcg.F90
Unescape Escape View File

@ -164,9 +164,29 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&

30
linsolve/impl/psb_cgcr.f90
Unescape Escape View File

@ -167,22 +167,30 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)

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

@ -97,7 +97,8 @@
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! 1: err = |r|/(|a||x|+|b|); here
! the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
@ -130,7 +131,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_c_vect_type) :: w, w1, xt
real(psb_spk_) :: tmp
complex(psb_spk_) :: scal, gm, rti, rti1
integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
& n_row, n_col, nl
integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@ -172,24 +173,35 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_,psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@ -253,12 +265,15 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
& w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then
bn2 = psb_genrm2(b,desc_a,info)
else if (istop_ == 3) then
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(cone,b,czero,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
@ -273,7 +288,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
r0n2 = psb_genrm2(v(1),desc_a,info)
endif
end select
errnum = czero
errden = cone
deps = eps
@ -324,27 +340,32 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
! check convergence
!
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = sone
case(psb_istop_rrn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = r0n2
endif
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart
if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@ -360,42 +381,18 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
call prec%apply(v(i),w1,desc_a,info)
call psb_spmm(cone,a,w1,czero,w,desc_a,info,work=aux)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
do k = 1, i
h(k,i) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,i),v(k),cone,w,desc_a,info)
end do
h(i+1,i) = psb_genrm2(w,desc_a,info)
scal=cone/h(i+1,i)
call psb_geaxpby(scal,w,czero,v(i+1),desc_a,info)
do k=2,i
call crot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_spk_),s(k-1))
enddo
rti = h(i,i)
rti1 = h(i+1,i)
call crotg(rti,rti1,tmp,s(i))
c(i) = cmplx(tmp,szero)
call crot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_spk_),s(i))
h(i+1,i) = czero
call crot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_spk_),s(i))
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
!
! build x and then compute the residual and its infinity norm
!
rst = rs
call w1%set(czero)
call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rst,size(rst,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rst(1:i)
do k=1, i
call psb_geaxpby(rst(k),v(k),cone,xt,desc_a,info)
end do
call prec%apply(xt,desc_a,info)
call psb_geaxpby(cone,x,cone,xt,desc_a,info)
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 psb_geaxpby(cone,b,czero,w1,desc_a,info)
call psb_spmm(-cone,a,xt,cone,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info)
@ -403,8 +400,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
errnum = rni
errden = (ani*xni+bni)
!
else if (istop_ == 2) then
case(psb_istop_bn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@ -412,7 +408,14 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = abs(rs(i+1))
errnum = rni
errden = sone
case(psb_istop_rrn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@ -420,28 +423,18 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = r0n2
endif
end select
if (errnum <= eps*errden) then
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(cone,xt,czero,x,desc_a,info)
! = x = xt
else if (istop_ == 2) then
!
! build x
!
call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(czero)
do k=1, i
call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(cone,w,cone,x,desc_a,info)
end if
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) !
end select
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@ -454,25 +447,15 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
end do inner
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(cone,xt,czero,x,desc_a,info)! x = xt
else if (istop_ == 2) then
!
! build x
!
call ctrsm('l','u','n','n',nl,1,cone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(czero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(cone,w,cone,x,desc_a,info)
end if
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) !
end select
if (itx >= litmax) then
if (itx >= itmax_) then
if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@ -502,6 +485,67 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
9999 call psb_error_handler(err_act)
return
contains
!
! Advance one step the Gram-Schmidt process, and apply
! Givens rotations to keep H triangular
!
subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
complex(psb_spk_) :: c(:), s(:), h(:,:), rs(:)
type(psb_c_vect_type) :: v(:), w
type(psb_desc_type) :: desc_a
complex(psb_spk_) :: scal, gm, rti, rti1
real(psb_spk_) :: tmp
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
do k = 1, n
h(k,n) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,n),v(k),cone,w,desc_a,info)
end do
h(n+1,n) = psb_genrm2(w,desc_a,info)
scal=cone/h(n+1,n)
call psb_geaxpby(scal,w,czero,v(n+1),desc_a,info)
do k=2,n
call crot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_spk_),s(k-1))
enddo
rti = h(n,n)
rti1 = h(n+1,n)
call crotg(rti,rti1,tmp,s(n))
c(n) = cmplx(tmp,szero)
call crot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_spk_),s(n))
call crot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_spk_),s(n))
end subroutine mgs
!
! 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)
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
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
!
! build x
!
call ctrsm('l','u','n','n',n,1,cone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(cone,w,cone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_crgmres_vect

24
linsolve/impl/psb_dbicg.f90
Unescape Escape View File

@ -160,19 +160,29 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
! istop_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=psb_err_invalid_istop_
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
endif
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then

23
linsolve/impl/psb_dcg.F90
Unescape Escape View File

@ -159,8 +159,29 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&

23
linsolve/impl/psb_dcgs.f90
Unescape Escape View File

@ -154,8 +154,29 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
istop_ = psb_get_istop_default()
Endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)

28
linsolve/impl/psb_dcgstab.f90
Unescape Escape View File

@ -156,13 +156,31 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
Endif
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b|| norm 2
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
! = if (.not.same_type_as(x,b)) then
! = write(0,*) 'Warning: different dynamic types for X and B '
! = end if

23
linsolve/impl/psb_dcgstabl.f90
Unescape Escape View File

@ -172,8 +172,29 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
itmax_ = itmax

24
linsolve/impl/psb_dfcg.F90
Unescape Escape View File

@ -164,9 +164,29 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&

30
linsolve/impl/psb_dgcr.f90
Unescape Escape View File

@ -167,22 +167,30 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)

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

@ -97,7 +97,8 @@
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! 1: err = |r|/(|a||x|+|b|); here
! the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
@ -130,7 +131,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_d_vect_type) :: w, w1, xt
real(psb_dpk_) :: tmp
real(psb_dpk_) :: scal, gm, rti, rti1
integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
& n_row, n_col, nl
integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@ -172,24 +173,35 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_,psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@ -253,12 +265,15 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
& w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then
bn2 = psb_genrm2(b,desc_a,info)
else if (istop_ == 3) then
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(done,b,dzero,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
@ -273,7 +288,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
r0n2 = psb_genrm2(v(1),desc_a,info)
endif
end select
errnum = dzero
errden = done
deps = eps
@ -324,27 +340,32 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
! check convergence
!
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = done
case(psb_istop_rrn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = r0n2
endif
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart
if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@ -360,42 +381,18 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
call prec%apply(v(i),w1,desc_a,info)
call psb_spmm(done,a,w1,dzero,w,desc_a,info,work=aux)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
do k = 1, i
h(k,i) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,i),v(k),done,w,desc_a,info)
end do
h(i+1,i) = psb_genrm2(w,desc_a,info)
scal=done/h(i+1,i)
call psb_geaxpby(scal,w,dzero,v(i+1),desc_a,info)
do k=2,i
call drot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_dpk_),s(k-1))
enddo
rti = h(i,i)
rti1 = h(i+1,i)
call drotg(rti,rti1,tmp,s(i))
c(i) = cmplx(tmp,dzero)
call drot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_dpk_),s(i))
h(i+1,i) = dzero
call drot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_dpk_),s(i))
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
!
! build x and then compute the residual and its infinity norm
!
rst = rs
call w1%set(dzero)
call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rst,size(rst,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rst(1:i)
do k=1, i
call psb_geaxpby(rst(k),v(k),done,xt,desc_a,info)
end do
call prec%apply(xt,desc_a,info)
call psb_geaxpby(done,x,done,xt,desc_a,info)
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 psb_geaxpby(done,b,dzero,w1,desc_a,info)
call psb_spmm(-done,a,xt,done,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info)
@ -403,8 +400,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
errnum = rni
errden = (ani*xni+bni)
!
else if (istop_ == 2) then
case(psb_istop_bn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@ -412,7 +408,14 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = abs(rs(i+1))
errnum = rni
errden = done
case(psb_istop_rrn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@ -420,28 +423,18 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = r0n2
endif
end select
if (errnum <= eps*errden) then
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(done,xt,dzero,x,desc_a,info)
! = x = xt
else if (istop_ == 2) then
!
! build x
!
call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(dzero)
do k=1, i
call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(done,w,done,x,desc_a,info)
end if
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) !
end select
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@ -454,25 +447,15 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
end do inner
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(done,xt,dzero,x,desc_a,info)! x = xt
else if (istop_ == 2) then
!
! build x
!
call dtrsm('l','u','n','n',nl,1,done,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(dzero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(done,w,done,x,desc_a,info)
end if
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) !
end select
if (itx >= litmax) then
if (itx >= itmax_) then
if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@ -502,6 +485,67 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
9999 call psb_error_handler(err_act)
return
contains
!
! Advance one step the Gram-Schmidt process, and apply
! Givens rotations to keep H triangular
!
subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
real(psb_dpk_) :: c(:), s(:), h(:,:), rs(:)
type(psb_d_vect_type) :: v(:), w
type(psb_desc_type) :: desc_a
real(psb_dpk_) :: scal, gm, rti, rti1
real(psb_dpk_) :: tmp
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
do k = 1, n
h(k,n) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,n),v(k),done,w,desc_a,info)
end do
h(n+1,n) = psb_genrm2(w,desc_a,info)
scal=done/h(n+1,n)
call psb_geaxpby(scal,w,dzero,v(n+1),desc_a,info)
do k=2,n
call drot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_dpk_),s(k-1))
enddo
rti = h(n,n)
rti1 = h(n+1,n)
call drotg(rti,rti1,tmp,s(n))
c(n) = cmplx(tmp,dzero)
call drot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_dpk_),s(n))
call drot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_dpk_),s(n))
end subroutine mgs
!
! 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)
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
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
!
! build x
!
call dtrsm('l','u','n','n',n,1,done,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(done,w,done,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_drgmres_vect

24
linsolve/impl/psb_sbicg.f90
Unescape Escape View File

@ -160,19 +160,29 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
! istop_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=psb_err_invalid_istop_
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
endif
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then

23
linsolve/impl/psb_scg.F90
Unescape Escape View File

@ -159,8 +159,29 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&

23
linsolve/impl/psb_scgs.f90
Unescape Escape View File

@ -154,8 +154,29 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
istop_ = psb_get_istop_default()
Endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)

28
linsolve/impl/psb_scgstab.f90
Unescape Escape View File

@ -156,13 +156,31 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
Endif
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b|| norm 2
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
! = if (.not.same_type_as(x,b)) then
! = write(0,*) 'Warning: different dynamic types for X and B '
! = end if

23
linsolve/impl/psb_scgstabl.f90
Unescape Escape View File

@ -172,8 +172,29 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
itmax_ = itmax

24
linsolve/impl/psb_sfcg.F90
Unescape Escape View File

@ -164,9 +164,29 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&

30
linsolve/impl/psb_sgcr.f90
Unescape Escape View File

@ -167,22 +167,30 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)

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

@ -97,7 +97,8 @@
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! 1: err = |r|/(|a||x|+|b|); here
! the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
@ -130,7 +131,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_s_vect_type) :: w, w1, xt
real(psb_spk_) :: tmp
real(psb_spk_) :: scal, gm, rti, rti1
integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
& n_row, n_col, nl
integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@ -172,24 +173,35 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_,psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@ -253,12 +265,15 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
& w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then
bn2 = psb_genrm2(b,desc_a,info)
else if (istop_ == 3) then
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(sone,b,szero,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
@ -273,7 +288,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
r0n2 = psb_genrm2(v(1),desc_a,info)
endif
end select
errnum = szero
errden = sone
deps = eps
@ -324,27 +340,32 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
! check convergence
!
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = sone
case(psb_istop_rrn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = r0n2
endif
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart
if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@ -360,42 +381,18 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
call prec%apply(v(i),w1,desc_a,info)
call psb_spmm(sone,a,w1,szero,w,desc_a,info,work=aux)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
do k = 1, i
h(k,i) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,i),v(k),sone,w,desc_a,info)
end do
h(i+1,i) = psb_genrm2(w,desc_a,info)
scal=sone/h(i+1,i)
call psb_geaxpby(scal,w,szero,v(i+1),desc_a,info)
do k=2,i
call srot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_spk_),s(k-1))
enddo
rti = h(i,i)
rti1 = h(i+1,i)
call srotg(rti,rti1,tmp,s(i))
c(i) = cmplx(tmp,szero)
call srot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_spk_),s(i))
h(i+1,i) = szero
call srot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_spk_),s(i))
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
!
! build x and then compute the residual and its infinity norm
!
rst = rs
call w1%set(szero)
call strsm('l','u','n','n',i,1,sone,h,size(h,1),rst,size(rst,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rst(1:i)
do k=1, i
call psb_geaxpby(rst(k),v(k),sone,xt,desc_a,info)
end do
call prec%apply(xt,desc_a,info)
call psb_geaxpby(sone,x,sone,xt,desc_a,info)
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 psb_geaxpby(sone,b,szero,w1,desc_a,info)
call psb_spmm(-sone,a,xt,sone,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info)
@ -403,8 +400,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
errnum = rni
errden = (ani*xni+bni)
!
else if (istop_ == 2) then
case(psb_istop_bn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@ -412,7 +408,14 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = abs(rs(i+1))
errnum = rni
errden = sone
case(psb_istop_rrn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@ -420,28 +423,18 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = r0n2
endif
end select
if (errnum <= eps*errden) then
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(sone,xt,szero,x,desc_a,info)
! = x = xt
else if (istop_ == 2) then
!
! build x
!
call strsm('l','u','n','n',i,1,sone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(szero)
do k=1, i
call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(sone,w,sone,x,desc_a,info)
end if
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) !
end select
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@ -454,25 +447,15 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
end do inner
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(sone,xt,szero,x,desc_a,info)! x = xt
else if (istop_ == 2) then
!
! build x
!
call strsm('l','u','n','n',nl,1,sone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(szero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(sone,w,sone,x,desc_a,info)
end if
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) !
end select
if (itx >= litmax) then
if (itx >= itmax_) then
if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@ -502,6 +485,67 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
9999 call psb_error_handler(err_act)
return
contains
!
! Advance one step the Gram-Schmidt process, and apply
! Givens rotations to keep H triangular
!
subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
real(psb_spk_) :: c(:), s(:), h(:,:), rs(:)
type(psb_s_vect_type) :: v(:), w
type(psb_desc_type) :: desc_a
real(psb_spk_) :: scal, gm, rti, rti1
real(psb_spk_) :: tmp
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
do k = 1, n
h(k,n) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,n),v(k),sone,w,desc_a,info)
end do
h(n+1,n) = psb_genrm2(w,desc_a,info)
scal=sone/h(n+1,n)
call psb_geaxpby(scal,w,szero,v(n+1),desc_a,info)
do k=2,n
call srot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_spk_),s(k-1))
enddo
rti = h(n,n)
rti1 = h(n+1,n)
call srotg(rti,rti1,tmp,s(n))
c(n) = cmplx(tmp,szero)
call srot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_spk_),s(n))
call srot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_spk_),s(n))
end subroutine mgs
!
! 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)
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
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
!
! build x
!
call strsm('l','u','n','n',n,1,sone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(sone,w,sone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_srgmres_vect

24
linsolve/impl/psb_zbicg.f90
Unescape Escape View File

@ -160,19 +160,29 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
! istop_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=psb_err_invalid_istop_
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
endif
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then

23
linsolve/impl/psb_zcg.F90
Unescape Escape View File

@ -159,8 +159,29 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&

23
linsolve/impl/psb_zcgs.f90
Unescape Escape View File

@ -154,8 +154,29 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
istop_ = psb_get_istop_default()
Endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)

28
linsolve/impl/psb_zcgstab.f90
Unescape Escape View File

@ -156,13 +156,31 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
Endif
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b|| norm 2
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
! = if (.not.same_type_as(x,b)) then
! = write(0,*) 'Warning: different dynamic types for X and B '
! = end if

23
linsolve/impl/psb_zcgstabl.f90
Unescape Escape View File

@ -172,8 +172,29 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
itmax_ = itmax

24
linsolve/impl/psb_zfcg.F90
Unescape Escape View File

@ -164,9 +164,29 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&

30
linsolve/impl/psb_zgcr.f90
Unescape Escape View File

@ -167,22 +167,30 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)

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

@ -97,7 +97,8 @@
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! 1: err = |r|/(|a||x|+|b|); here
! the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
@ -130,7 +131,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_z_vect_type) :: w, w1, xt
real(psb_dpk_) :: tmp
complex(psb_dpk_) :: scal, gm, rti, rti1
integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
& n_row, n_col, nl
integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@ -172,24 +173,35 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_,psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@ -253,12 +265,15 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
& w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then
bn2 = psb_genrm2(b,desc_a,info)
else if (istop_ == 3) then
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(zone,b,zzero,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
@ -273,7 +288,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
r0n2 = psb_genrm2(v(1),desc_a,info)
endif
end select
errnum = zzero
errden = zone
deps = eps
@ -324,27 +340,32 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
! check convergence
!
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = done
case(psb_istop_rrn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = r0n2
endif
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart
if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@ -360,42 +381,18 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
call prec%apply(v(i),w1,desc_a,info)
call psb_spmm(zone,a,w1,zzero,w,desc_a,info,work=aux)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
do k = 1, i
h(k,i) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,i),v(k),zone,w,desc_a,info)
end do
h(i+1,i) = psb_genrm2(w,desc_a,info)
scal=zone/h(i+1,i)
call psb_geaxpby(scal,w,zzero,v(i+1),desc_a,info)
do k=2,i
call zrot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_dpk_),s(k-1))
enddo
rti = h(i,i)
rti1 = h(i+1,i)
call zrotg(rti,rti1,tmp,s(i))
c(i) = cmplx(tmp,dzero)
call zrot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_dpk_),s(i))
h(i+1,i) = zzero
call zrot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_dpk_),s(i))
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
!
! build x and then compute the residual and its infinity norm
!
rst = rs
call w1%set(zzero)
call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rst,size(rst,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rst(1:i)
do k=1, i
call psb_geaxpby(rst(k),v(k),zone,xt,desc_a,info)
end do
call prec%apply(xt,desc_a,info)
call psb_geaxpby(zone,x,zone,xt,desc_a,info)
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 psb_geaxpby(zone,b,zzero,w1,desc_a,info)
call psb_spmm(-zone,a,xt,zone,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info)
@ -403,8 +400,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
errnum = rni
errden = (ani*xni+bni)
!
else if (istop_ == 2) then
case(psb_istop_bn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@ -412,7 +408,14 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = abs(rs(i+1))
errnum = rni
errden = done
case(psb_istop_rrn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@ -420,28 +423,18 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = r0n2
endif
end select
if (errnum <= eps*errden) then
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(zone,xt,zzero,x,desc_a,info)
! = x = xt
else if (istop_ == 2) then
!
! build x
!
call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(zzero)
do k=1, i
call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(zone,w,zone,x,desc_a,info)
end if
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) !
end select
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@ -454,25 +447,15 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
end do inner
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(zone,xt,zzero,x,desc_a,info)! x = xt
else if (istop_ == 2) then
!
! build x
!
call ztrsm('l','u','n','n',nl,1,zone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(zzero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(zone,w,zone,x,desc_a,info)
end if
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) !
end select
if (itx >= litmax) then
if (itx >= itmax_) then
if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@ -502,6 +485,67 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
9999 call psb_error_handler(err_act)
return
contains
!
! Advance one step the Gram-Schmidt process, and apply
! Givens rotations to keep H triangular
!
subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
complex(psb_dpk_) :: c(:), s(:), h(:,:), rs(:)
type(psb_z_vect_type) :: v(:), w
type(psb_desc_type) :: desc_a
complex(psb_dpk_) :: scal, gm, rti, rti1
real(psb_dpk_) :: tmp
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
do k = 1, n
h(k,n) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,n),v(k),zone,w,desc_a,info)
end do
h(n+1,n) = psb_genrm2(w,desc_a,info)
scal=zone/h(n+1,n)
call psb_geaxpby(scal,w,zzero,v(n+1),desc_a,info)
do k=2,n
call zrot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_dpk_),s(k-1))
enddo
rti = h(n,n)
rti1 = h(n+1,n)
call zrotg(rti,rti1,tmp,s(n))
c(n) = cmplx(tmp,dzero)
call zrot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_dpk_),s(n))
call zrot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_dpk_),s(n))
end subroutine mgs
!
! 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)
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
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
!
! build x
!
call ztrsm('l','u','n','n',n,1,zone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(zone,w,zone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_zrgmres_vect

33
linsolve/psb_base_linsolve_conv_mod.f90
Unescape Escape View File

@ -41,10 +41,20 @@ Module psb_base_linsolve_conv_mod
module procedure psb_d_end_conv
end interface
!
integer(psb_ipk_), parameter :: psb_istop_min_ = 1
integer(psb_ipk_), parameter :: psb_istop_ani_ = 1
integer(psb_ipk_), parameter :: psb_istop_bn2_ = 2
integer(psb_ipk_), parameter :: psb_istop_rn2_abs_ = 3
integer(psb_ipk_), parameter :: psb_istop_rrn2_ = 4
integer(psb_ipk_), parameter :: psb_istop_scbn2_ = 5
integer(psb_ipk_), parameter :: psb_istop_max_ = 3
! Fields in controls and values
integer(psb_ipk_), parameter :: psb_ik_bni_=1, psb_ik_rni_=2, psb_ik_ani_=3
integer(psb_ipk_), parameter :: psb_ik_xni_=4, psb_ik_bn2_=5, psb_ik_r0n2_=6
integer(psb_ipk_), parameter :: psb_ik_xn2_=7, psb_ik_errnum_=8, psb_ik_errden_=9
integer(psb_ipk_), parameter :: psb_ik_eps_=10, psb_ik_rn2_=11
integer(psb_ipk_), parameter :: psb_ik_eps_=10, psb_ik_rn2_=11, psb_ik_rn2_abs_=12
integer(psb_ipk_), parameter :: psb_ik_stopc_=1, psb_ik_trace_=2, psb_ik_itmax_=3
integer(psb_ipk_), parameter :: psb_ik_ivsz_=16
type psb_itconv_type
@ -52,8 +62,29 @@ Module psb_base_linsolve_conv_mod
real(psb_dpk_) :: values(psb_ik_ivsz_)
end type psb_itconv_type
integer(psb_ipk_), save :: psb_istop_default = psb_istop_bn2_
contains
function psb_is_valid_istop(istop) result(res)
integer(psb_ipk_) :: istop
logical :: res
res = ((psb_istop_min_<=istop).and.(istop<=psb_istop_max_))
end function psb_is_valid_istop
function psb_get_istop_default() result(res)
integer(psb_ipk_) :: res
res = psb_istop_default
end function psb_get_istop_default
subroutine psb_set_istop_default(val)
integer(psb_ipk_) :: val
if ((psb_istop_min_<=val).and.(val<=psb_istop_max_)) &
& psb_istop_default = val
end subroutine psb_set_istop_default
subroutine log_header(methdname)
!use psb_base_mod
implicit none

Write Preview
Loading…
Cancel
Save