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

Refactor RGMRES step 1.

Browse Source
kinsol-stop
sfilippone 11 hours ago
parent 4ec9c5cfd3
commit a424441677
5 changed files with 332 additions and 383 deletions
Show Stats Download Patch File Download Diff File

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

@ -134,7 +134,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.
@ -202,9 +202,9 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@ -369,7 +369,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
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)
@ -391,26 +391,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
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)
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))
select case(istop_)
case(psb_istop_ani_)
@ -418,16 +400,9 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
! 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)
@ -466,35 +441,9 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
case(psb_istop_ani_)
call psb_geaxpby(cone,xt,czero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_)
!
! 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)
case(psb_istop_rrn2_)
!
! build x
!
call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(czero)
do k=1, i
call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(cone,w,cone,x,desc_a,info)
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) &
@ -512,38 +461,11 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
case(psb_istop_ani_)
call psb_geaxpby(cone,xt,czero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_)
!
! 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)
case(psb_istop_rrn2_)
!
! build x
!
call ctrsm('l','u','n','n',nl,1,cone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(czero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(cone,w,cone,x,desc_a,info)
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)
@ -573,6 +495,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

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

@ -134,7 +134,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.
@ -202,9 +202,9 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@ -369,7 +369,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
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)
@ -391,26 +391,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
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)
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))
select case(istop_)
case(psb_istop_ani_)
@ -418,16 +400,9 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
! 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)
@ -466,35 +441,9 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
case(psb_istop_ani_)
call psb_geaxpby(done,xt,dzero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_)
!
! 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)
case(psb_istop_rrn2_)
!
! build x
!
call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(dzero)
do k=1, i
call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(done,w,done,x,desc_a,info)
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) &
@ -512,38 +461,11 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
case(psb_istop_ani_)
call psb_geaxpby(done,xt,dzero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_)
!
! 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)
case(psb_istop_rrn2_)
!
! build x
!
call dtrsm('l','u','n','n',nl,1,done,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(dzero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(done,w,done,x,desc_a,info)
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)
@ -573,6 +495,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

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

@ -134,7 +134,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.
@ -202,9 +202,9 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@ -369,7 +369,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
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)
@ -391,26 +391,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
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)
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))
select case(istop_)
case(psb_istop_ani_)
@ -418,16 +400,9 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
! 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)
@ -466,35 +441,9 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
case(psb_istop_ani_)
call psb_geaxpby(sone,xt,szero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_)
!
! 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)
case(psb_istop_rrn2_)
!
! build x
!
call strsm('l','u','n','n',i,1,sone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(szero)
do k=1, i
call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(sone,w,sone,x,desc_a,info)
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) &
@ -512,38 +461,11 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
case(psb_istop_ani_)
call psb_geaxpby(sone,xt,szero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_)
!
! 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)
case(psb_istop_rrn2_)
!
! build x
!
call strsm('l','u','n','n',nl,1,sone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(szero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(sone,w,sone,x,desc_a,info)
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)
@ -573,6 +495,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

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

@ -134,7 +134,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.
@ -202,9 +202,9 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@ -369,7 +369,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
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)
@ -391,26 +391,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
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)
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))
select case(istop_)
case(psb_istop_ani_)
@ -418,16 +400,9 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
! 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)
@ -466,35 +441,9 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
case(psb_istop_ani_)
call psb_geaxpby(zone,xt,zzero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_)
!
! 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)
case(psb_istop_rrn2_)
!
! build x
!
call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(zzero)
do k=1, i
call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(zone,w,zone,x,desc_a,info)
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) &
@ -512,38 +461,11 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
case(psb_istop_ani_)
call psb_geaxpby(zone,xt,zzero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_)
!
! 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)
case(psb_istop_rrn2_)
!
! build x
!
call ztrsm('l','u','n','n',nl,1,zone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(zzero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(zone,w,zone,x,desc_a,info)
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)
@ -573,6 +495,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

47
util/psb_d_mmio_impl.f90
Unescape Escape View File

@ -39,8 +39,10 @@ subroutine mm_dvet_read(b, info, iunit, filename)
integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile
character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024
logical :: opened
info = psb_success_
opened = .false.
if (present(filename)) then
if (filename == '-') then
infile=5
@ -51,6 +53,7 @@ subroutine mm_dvet_read(b, info, iunit, filename)
infile=99
endif
open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif
else
if (present(iunit)) then
@ -83,7 +86,7 @@ subroutine mm_dvet_read(b, info, iunit, filename)
end do
end if ! read right hand sides
if (infile /= 5) close(infile)
if (opened) close(infile)
return
! open failed
@ -111,8 +114,10 @@ subroutine mm_dvet2_read(b, info, iunit, filename)
integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile
character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024
logical :: opened
info = psb_success_
opened = .false.
if (present(filename)) then
if (filename == '-') then
infile=5
@ -123,6 +128,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename)
infile=99
endif
open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif
else
if (present(iunit)) then
@ -153,7 +159,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename)
read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol)
end if ! read right hand sides
if (infile /= 5) close(infile)
if (opened) close(infile)
return
! open failed
@ -182,8 +188,10 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename)
integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile
character(len=80) :: frmtv
logical :: opened
info = psb_success_
opened = .false.
if (present(filename)) then
if (filename == '-') then
outfile=6
@ -194,6 +202,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename)
outfile=99
endif
open(outfile,file=filename, err=901, action='WRITE')
opened = .true.
endif
else
if (present(iunit)) then
@ -212,7 +221,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename)
write(outfile,fmt='(es26.18,1x)') ((b(i,j), i=1,nrow),j=1,ncol)
if (outfile /= 6) close(outfile)
if (opened) close(outfile)
return
! open failed
@ -234,8 +243,10 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename)
integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile
character(len=80) :: frmtv
logical :: opened
info = psb_success_
opened = .false.
if (present(filename)) then
if (filename == '-') then
outfile=6
@ -246,6 +257,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename)
outfile=99
endif
open(outfile,file=filename, err=901, action='WRITE')
opened = .true.
endif
else
if (present(iunit)) then
@ -268,7 +280,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename)
write(outfile,frmtv) b(i)
end do
if (outfile /= 6) close(outfile)
if (opened) close(outfile)
return
! open failed
@ -324,9 +336,10 @@ subroutine dmm_mat_read(a, info, iunit, filename)
integer(psb_ipk_) :: nrow, ncol, nnzero
integer(psb_ipk_) :: ircode, i,nzr,infile
type(psb_d_coo_sparse_mat), allocatable :: acoo
logical :: opened
info = psb_success_
opened = .false.
if (present(filename)) then
if (filename == '-') then
infile=5
@ -337,6 +350,7 @@ subroutine dmm_mat_read(a, info, iunit, filename)
infile=99
endif
open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif
else
if (present(iunit)) then
@ -421,7 +435,7 @@ subroutine dmm_mat_read(a, info, iunit, filename)
call a%mv_from(acoo)
end if
if (infile /= 5) close(infile)
if (opened) close(infile)
return
! open failed
@ -449,10 +463,11 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename)
integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename
integer(psb_ipk_) :: iout
logical :: opened
info = psb_success_
opened = .false.
if (present(filename)) then
if (filename == '-') then
iout=6
@ -463,6 +478,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename)
iout=99
endif
open(iout,file=filename, err=901, action='WRITE')
opened = .true.
endif
else
if (present(iunit)) then
@ -474,7 +490,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename)
call a%print(iout,head=mtitle)
if (iout /= 6) close(iout)
if (opened) close(iout)
return
@ -498,9 +514,10 @@ subroutine ldmm_mat_read(a, info, iunit, filename)
integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr
integer(psb_ipk_) :: ircode, infile
type(psb_ld_coo_sparse_mat), allocatable :: acoo
logical :: opened
info = psb_success_
opened = .false.
if (present(filename)) then
if (filename == '-') then
infile=5
@ -511,6 +528,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename)
infile=99
endif
open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif
else
if (present(iunit)) then
@ -595,7 +613,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename)
call a%mv_from(acoo)
end if
if (infile /= 5) close(infile)
if (opened) close(infile)
return
! open failed
@ -623,10 +641,10 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename
integer(psb_ipk_) :: iout
logical :: opened
info = psb_success_
opened = .false.
if (present(filename)) then
if (filename == '-') then
iout=6
@ -637,6 +655,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
iout=99
endif
open(iout,file=filename, err=901, action='WRITE')
opened = .true.
endif
else
if (present(iunit)) then
@ -648,7 +667,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
call a%print(iout,head=mtitle)
if (iout /= 6) close(iout)
if (opened) close(iout)
return
@ -658,5 +677,3 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
write(psb_err_unit,*) 'Error while opening ',filename
return
end subroutine ldmm_mat_write

Write Preview
Loading…
Cancel
Save