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

psblas:

Browse Source 
krylov/psb_ccgr.f90
 krylov/psb_dcgr.f90
 krylov/psb_scgr.f90
 krylov/psb_zcgr.f90

Cosmetic changes
psblas3-pattern
Salvatore Filippone 9 years ago
parent 5edd786cc6
commit c83985fe16
4 changed files with 396 additions and 436 deletions
Show Stats Download Patch File Download Diff File

98
krylov/psb_ccgr.f90
Unescape Escape View File

@ -124,12 +124,12 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop
integer(psb_ipk_), Optional, Intent(out) :: iter integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_spk_), Optional, Intent(out) :: err real(psb_spk_), Optional, Intent(out) :: err
! = local data ! = local data
complex(psb_spk_), allocatable :: alpha(:), h(:,:) complex(psb_spk_), allocatable :: alpha(:), h(:,:)
type(psb_c_vect_type), allocatable :: z(:), c(:), c_scale(:) type(psb_c_vect_type), allocatable :: z(:), c(:), c_scale(:)
type(psb_c_vect_type) :: r type(psb_c_vect_type) :: r
real(psb_dpk_) :: r_norm, b_norm, a_norm, x_norm, derr real(psb_dpk_) :: r_norm, b_norm, a_norm, derr
integer(psb_ipk_) :: n_col, mglob, naux, err_act integer(psb_ipk_) :: n_col, mglob, naux, err_act
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: np, me, ictxt integer(psb_ipk_) :: np, me, ictxt
@ -140,7 +140,7 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_itconv_type) :: stopdat type(psb_itconv_type) :: stopdat
character(len=*), parameter :: methdname='CGR' character(len=*), parameter :: methdname='CGR'
integer(psb_ipk_) ::int_err(5) integer(psb_ipk_) ::int_err(5)
info = psb_success_ info = psb_success_
name = 'psb_ccgr' name = 'psb_ccgr'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
@ -169,10 +169,10 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,&
istop_ = 2 istop_ = 2
endif endif
! !
! ISTOP_ = 1: Normwise backward error, infinity norm ! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm ! ISTOP_ = 2: ||r||/||b||, 2-norm
! !
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=psb_err_invalid_istop_ info=psb_err_invalid_istop_
@ -235,28 +235,18 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999 goto 9999
end if end if
call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v) call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v)
x_norm = psb_normi(x, desc_a, info)
do i =1,nl+1
call psb_geasb(c(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(z(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(c_scale(i), desc_a,info, scratch=.true.,mold=x%v)
end do
do i =1,nl+1 itx = 0
call psb_geasb(c(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(z(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(c_scale(i), desc_a,info, scratch=.true.,mold=x%v)
end do
itx = 0
if (istop_ == 2) then nrst = -1
b_norm = psb_norm2(b, desc_a, info) call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info)
else if (istop_ == 1) then
a_norm = psb_spnrmi(a,desc_a,info)
b_norm = psb_normi(b, desc_a, info)
endif
nrst = -1
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info)
restart: do restart: do
if (itx>= itmax_) exit restart if (itx>= itmax_) exit restart
h = czero h = czero
@ -285,54 +275,54 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,&
! goto 9999 ! goto 9999
! End If ! End If
nrst = nrst + 1 nrst = nrst + 1
iteration: do iteration: do
itx = itx + 1 itx = itx + 1
it = it + 1 it = it + 1
j = it j = it
!Apply preconditioner !Apply preconditioner
call prec%apply(r,z(j),desc_a,info,work=aux) call prec%apply(r,z(j),desc_a,info,work=aux)
call psb_spmm(cone,a,z(j),czero,c(1),desc_a,info,work=aux) call psb_spmm(cone,a,z(j),czero,c(1),desc_a,info,work=aux)
do i =1, j - 1 do i =1, j - 1
h(i,j) = psb_gedot(c_scale(i), c(i), desc_a, info) h(i,j) = psb_gedot(c_scale(i), c(i), desc_a, info)
call psb_geaxpby(cone, c(i), czero, c(i+1), desc_a, info) call psb_geaxpby(cone, c(i), czero, c(i+1), desc_a, info)
call psb_geaxpby(-h(i,j), c_scale(i), cone, c(i+1), desc_a, info) call psb_geaxpby(-h(i,j), c_scale(i), cone, c(i+1), desc_a, info)
end do end do
h(j,j) = psb_norm2(c(j), desc_a, info) h(j,j) = psb_norm2(c(j), desc_a, info)
hjj = cone/h(j,j) hjj = cone/h(j,j)
call psb_geaxpby(hjj, c(j), czero, c_scale(j), desc_a, info) call psb_geaxpby(hjj, c(j), czero, c_scale(j), desc_a, info)
alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) alpha(j) = psb_gedot(c_scale(j), r, desc_a, info)
!Update residual !Update residual
call psb_geaxpby(cone, r, czero, r, desc_a, info) call psb_geaxpby(cone, r, czero, r, desc_a, info)
call psb_geaxpby(-alpha(j), c_scale(j), cone, r, desc_a, info) call psb_geaxpby(-alpha(j), c_scale(j), cone, r, desc_a, info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (j >= irst) exit iteration if (j >= irst) exit iteration
end do iteration end do iteration
m = j m = j
!Compute solution !Compute solution
call ctrsm('l','u','n','n',m,1,cone,h,size(h,1),alpha,size(alpha,1)) call ctrsm('l','u','n','n',m,1,cone,h,size(h,1),alpha,size(alpha,1))
if (nrst == 0 ) then if (nrst == 0 ) then
call x%set(czero) call x%set(czero)
endif endif
do i=1,m do i=1,m
call psb_geaxpby(alpha(i), z(i), cone, x, desc_a, info) call psb_geaxpby(alpha(i), z(i), cone, x, desc_a, info)
enddo enddo

98
krylov/psb_dcgr.f90
Unescape Escape View File

@ -124,12 +124,12 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop
integer(psb_ipk_), Optional, Intent(out) :: iter integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_dpk_), Optional, Intent(out) :: err real(psb_dpk_), Optional, Intent(out) :: err
! = local data ! = local data
real(psb_dpk_), allocatable :: alpha(:), h(:,:) real(psb_dpk_), allocatable :: alpha(:), h(:,:)
type(psb_d_vect_type), allocatable :: z(:), c(:), c_scale(:) type(psb_d_vect_type), allocatable :: z(:), c(:), c_scale(:)
type(psb_d_vect_type) :: r type(psb_d_vect_type) :: r
real(psb_dpk_) :: r_norm, b_norm, a_norm, x_norm, derr real(psb_dpk_) :: r_norm, b_norm, a_norm, derr
integer(psb_ipk_) :: n_col, mglob, naux, err_act integer(psb_ipk_) :: n_col, mglob, naux, err_act
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: np, me, ictxt integer(psb_ipk_) :: np, me, ictxt
@ -140,7 +140,7 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_itconv_type) :: stopdat type(psb_itconv_type) :: stopdat
character(len=*), parameter :: methdname='CGR' character(len=*), parameter :: methdname='CGR'
integer(psb_ipk_) ::int_err(5) integer(psb_ipk_) ::int_err(5)
info = psb_success_ info = psb_success_
name = 'psb_dcgr' name = 'psb_dcgr'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
@ -169,10 +169,10 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,&
istop_ = 2 istop_ = 2
endif endif
! !
! ISTOP_ = 1: Normwise backward error, infinity norm ! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm ! ISTOP_ = 2: ||r||/||b||, 2-norm
! !
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=psb_err_invalid_istop_ info=psb_err_invalid_istop_
@ -235,28 +235,18 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999 goto 9999
end if end if
call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v) call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v)
x_norm = psb_normi(x, desc_a, info)
do i =1,nl+1
call psb_geasb(c(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(z(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(c_scale(i), desc_a,info, scratch=.true.,mold=x%v)
end do
do i =1,nl+1 itx = 0
call psb_geasb(c(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(z(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(c_scale(i), desc_a,info, scratch=.true.,mold=x%v)
end do
itx = 0
if (istop_ == 2) then nrst = -1
b_norm = psb_norm2(b, desc_a, info) call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info)
else if (istop_ == 1) then
a_norm = psb_spnrmi(a,desc_a,info)
b_norm = psb_normi(b, desc_a, info)
endif
nrst = -1
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info)
restart: do restart: do
if (itx>= itmax_) exit restart if (itx>= itmax_) exit restart
h = dzero h = dzero
@ -285,54 +275,54 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,&
! goto 9999 ! goto 9999
! End If ! End If
nrst = nrst + 1 nrst = nrst + 1
iteration: do iteration: do
itx = itx + 1 itx = itx + 1
it = it + 1 it = it + 1
j = it j = it
!Apply preconditioner !Apply preconditioner
call prec%apply(r,z(j),desc_a,info,work=aux) call prec%apply(r,z(j),desc_a,info,work=aux)
call psb_spmm(done,a,z(j),dzero,c(1),desc_a,info,work=aux) call psb_spmm(done,a,z(j),dzero,c(1),desc_a,info,work=aux)
do i =1, j - 1 do i =1, j - 1
h(i,j) = psb_gedot(c_scale(i), c(i), desc_a, info) h(i,j) = psb_gedot(c_scale(i), c(i), desc_a, info)
call psb_geaxpby(done, c(i), dzero, c(i+1), desc_a, info) call psb_geaxpby(done, c(i), dzero, c(i+1), desc_a, info)
call psb_geaxpby(-h(i,j), c_scale(i), done, c(i+1), desc_a, info) call psb_geaxpby(-h(i,j), c_scale(i), done, c(i+1), desc_a, info)
end do end do
h(j,j) = psb_norm2(c(j), desc_a, info) h(j,j) = psb_norm2(c(j), desc_a, info)
hjj = done/h(j,j) hjj = done/h(j,j)
call psb_geaxpby(hjj, c(j), dzero, c_scale(j), desc_a, info) call psb_geaxpby(hjj, c(j), dzero, c_scale(j), desc_a, info)
alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) alpha(j) = psb_gedot(c_scale(j), r, desc_a, info)
!Update residual !Update residual
call psb_geaxpby(done, r, dzero, r, desc_a, info) call psb_geaxpby(done, r, dzero, r, desc_a, info)
call psb_geaxpby(-alpha(j), c_scale(j), done, r, desc_a, info) call psb_geaxpby(-alpha(j), c_scale(j), done, r, desc_a, info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (j >= irst) exit iteration if (j >= irst) exit iteration
end do iteration end do iteration
m = j m = j
!Compute solution !Compute solution
call dtrsm('l','u','n','n',m,1,done,h,size(h,1),alpha,size(alpha,1)) call dtrsm('l','u','n','n',m,1,done,h,size(h,1),alpha,size(alpha,1))
if (nrst == 0 ) then if (nrst == 0 ) then
call x%set(dzero) call x%set(dzero)
endif endif
do i=1,m do i=1,m
call psb_geaxpby(alpha(i), z(i), done, x, desc_a, info) call psb_geaxpby(alpha(i), z(i), done, x, desc_a, info)
enddo enddo

98
krylov/psb_scgr.f90
Unescape Escape View File

@ -124,12 +124,12 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop
integer(psb_ipk_), Optional, Intent(out) :: iter integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_spk_), Optional, Intent(out) :: err real(psb_spk_), Optional, Intent(out) :: err
! = local data ! = local data
real(psb_spk_), allocatable :: alpha(:), h(:,:) real(psb_spk_), allocatable :: alpha(:), h(:,:)
type(psb_s_vect_type), allocatable :: z(:), c(:), c_scale(:) type(psb_s_vect_type), allocatable :: z(:), c(:), c_scale(:)
type(psb_s_vect_type) :: r type(psb_s_vect_type) :: r
real(psb_dpk_) :: r_norm, b_norm, a_norm, x_norm, derr real(psb_dpk_) :: r_norm, b_norm, a_norm, derr
integer(psb_ipk_) :: n_col, mglob, naux, err_act integer(psb_ipk_) :: n_col, mglob, naux, err_act
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: np, me, ictxt integer(psb_ipk_) :: np, me, ictxt
@ -140,7 +140,7 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_itconv_type) :: stopdat type(psb_itconv_type) :: stopdat
character(len=*), parameter :: methdname='CGR' character(len=*), parameter :: methdname='CGR'
integer(psb_ipk_) ::int_err(5) integer(psb_ipk_) ::int_err(5)
info = psb_success_ info = psb_success_
name = 'psb_scgr' name = 'psb_scgr'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
@ -169,10 +169,10 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,&
istop_ = 2 istop_ = 2
endif endif
! !
! ISTOP_ = 1: Normwise backward error, infinity norm ! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm ! ISTOP_ = 2: ||r||/||b||, 2-norm
! !
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=psb_err_invalid_istop_ info=psb_err_invalid_istop_
@ -235,28 +235,18 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999 goto 9999
end if end if
call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v) call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v)
x_norm = psb_normi(x, desc_a, info)
do i =1,nl+1
call psb_geasb(c(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(z(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(c_scale(i), desc_a,info, scratch=.true.,mold=x%v)
end do
do i =1,nl+1 itx = 0
call psb_geasb(c(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(z(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(c_scale(i), desc_a,info, scratch=.true.,mold=x%v)
end do
itx = 0
if (istop_ == 2) then nrst = -1
b_norm = psb_norm2(b, desc_a, info) call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info)
else if (istop_ == 1) then
a_norm = psb_spnrmi(a,desc_a,info)
b_norm = psb_normi(b, desc_a, info)
endif
nrst = -1
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info)
restart: do restart: do
if (itx>= itmax_) exit restart if (itx>= itmax_) exit restart
h = szero h = szero
@ -285,54 +275,54 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,&
! goto 9999 ! goto 9999
! End If ! End If
nrst = nrst + 1 nrst = nrst + 1
iteration: do iteration: do
itx = itx + 1 itx = itx + 1
it = it + 1 it = it + 1
j = it j = it
!Apply preconditioner !Apply preconditioner
call prec%apply(r,z(j),desc_a,info,work=aux) call prec%apply(r,z(j),desc_a,info,work=aux)
call psb_spmm(sone,a,z(j),szero,c(1),desc_a,info,work=aux) call psb_spmm(sone,a,z(j),szero,c(1),desc_a,info,work=aux)
do i =1, j - 1 do i =1, j - 1
h(i,j) = psb_gedot(c_scale(i), c(i), desc_a, info) h(i,j) = psb_gedot(c_scale(i), c(i), desc_a, info)
call psb_geaxpby(sone, c(i), szero, c(i+1), desc_a, info) call psb_geaxpby(sone, c(i), szero, c(i+1), desc_a, info)
call psb_geaxpby(-h(i,j), c_scale(i), sone, c(i+1), desc_a, info) call psb_geaxpby(-h(i,j), c_scale(i), sone, c(i+1), desc_a, info)
end do end do
h(j,j) = psb_norm2(c(j), desc_a, info) h(j,j) = psb_norm2(c(j), desc_a, info)
hjj = sone/h(j,j) hjj = sone/h(j,j)
call psb_geaxpby(hjj, c(j), szero, c_scale(j), desc_a, info) call psb_geaxpby(hjj, c(j), szero, c_scale(j), desc_a, info)
alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) alpha(j) = psb_gedot(c_scale(j), r, desc_a, info)
!Update residual !Update residual
call psb_geaxpby(sone, r, szero, r, desc_a, info) call psb_geaxpby(sone, r, szero, r, desc_a, info)
call psb_geaxpby(-alpha(j), c_scale(j), sone, r, desc_a, info) call psb_geaxpby(-alpha(j), c_scale(j), sone, r, desc_a, info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (j >= irst) exit iteration if (j >= irst) exit iteration
end do iteration end do iteration
m = j m = j
!Compute solution !Compute solution
call strsm('l','u','n','n',m,1,sone,h,size(h,1),alpha,size(alpha,1)) call strsm('l','u','n','n',m,1,sone,h,size(h,1),alpha,size(alpha,1))
if (nrst == 0 ) then if (nrst == 0 ) then
call x%set(szero) call x%set(szero)
endif endif
do i=1,m do i=1,m
call psb_geaxpby(alpha(i), z(i), sone, x, desc_a, info) call psb_geaxpby(alpha(i), z(i), sone, x, desc_a, info)
enddo enddo

98
krylov/psb_zcgr.f90
Unescape Escape View File

@ -124,12 +124,12 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop
integer(psb_ipk_), Optional, Intent(out) :: iter integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_dpk_), Optional, Intent(out) :: err real(psb_dpk_), Optional, Intent(out) :: err
! = local data ! = local data
complex(psb_dpk_), allocatable :: alpha(:), h(:,:) complex(psb_dpk_), allocatable :: alpha(:), h(:,:)
type(psb_z_vect_type), allocatable :: z(:), c(:), c_scale(:) type(psb_z_vect_type), allocatable :: z(:), c(:), c_scale(:)
type(psb_z_vect_type) :: r type(psb_z_vect_type) :: r
real(psb_dpk_) :: r_norm, b_norm, a_norm, x_norm, derr real(psb_dpk_) :: r_norm, b_norm, a_norm, derr
integer(psb_ipk_) :: n_col, mglob, naux, err_act integer(psb_ipk_) :: n_col, mglob, naux, err_act
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: np, me, ictxt integer(psb_ipk_) :: np, me, ictxt
@ -140,7 +140,7 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_itconv_type) :: stopdat type(psb_itconv_type) :: stopdat
character(len=*), parameter :: methdname='CGR' character(len=*), parameter :: methdname='CGR'
integer(psb_ipk_) ::int_err(5) integer(psb_ipk_) ::int_err(5)
info = psb_success_ info = psb_success_
name = 'psb_zcgr' name = 'psb_zcgr'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
@ -169,10 +169,10 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,&
istop_ = 2 istop_ = 2
endif endif
! !
! ISTOP_ = 1: Normwise backward error, infinity norm ! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm ! ISTOP_ = 2: ||r||/||b||, 2-norm
! !
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=psb_err_invalid_istop_ info=psb_err_invalid_istop_
@ -235,28 +235,18 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999 goto 9999
end if end if
call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v) call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v)
x_norm = psb_normi(x, desc_a, info)
do i =1,nl+1
call psb_geasb(c(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(z(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(c_scale(i), desc_a,info, scratch=.true.,mold=x%v)
end do
do i =1,nl+1 itx = 0
call psb_geasb(c(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(z(i), desc_a,info, scratch=.true.,mold=x%v)
call psb_geasb(c_scale(i), desc_a,info, scratch=.true.,mold=x%v)
end do
itx = 0
if (istop_ == 2) then nrst = -1
b_norm = psb_norm2(b, desc_a, info) call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info)
else if (istop_ == 1) then
a_norm = psb_spnrmi(a,desc_a,info)
b_norm = psb_normi(b, desc_a, info)
endif
nrst = -1
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info)
restart: do restart: do
if (itx>= itmax_) exit restart if (itx>= itmax_) exit restart
h = zzero h = zzero
@ -285,54 +275,54 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,&
! goto 9999 ! goto 9999
! End If ! End If
nrst = nrst + 1 nrst = nrst + 1
iteration: do iteration: do
itx = itx + 1 itx = itx + 1
it = it + 1 it = it + 1
j = it j = it
!Apply preconditioner !Apply preconditioner
call prec%apply(r,z(j),desc_a,info,work=aux) call prec%apply(r,z(j),desc_a,info,work=aux)
call psb_spmm(zone,a,z(j),zzero,c(1),desc_a,info,work=aux) call psb_spmm(zone,a,z(j),zzero,c(1),desc_a,info,work=aux)
do i =1, j - 1 do i =1, j - 1
h(i,j) = psb_gedot(c_scale(i), c(i), desc_a, info) h(i,j) = psb_gedot(c_scale(i), c(i), desc_a, info)
call psb_geaxpby(zone, c(i), zzero, c(i+1), desc_a, info) call psb_geaxpby(zone, c(i), zzero, c(i+1), desc_a, info)
call psb_geaxpby(-h(i,j), c_scale(i), zone, c(i+1), desc_a, info) call psb_geaxpby(-h(i,j), c_scale(i), zone, c(i+1), desc_a, info)
end do end do
h(j,j) = psb_norm2(c(j), desc_a, info) h(j,j) = psb_norm2(c(j), desc_a, info)
hjj = zone/h(j,j) hjj = zone/h(j,j)
call psb_geaxpby(hjj, c(j), zzero, c_scale(j), desc_a, info) call psb_geaxpby(hjj, c(j), zzero, c_scale(j), desc_a, info)
alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) alpha(j) = psb_gedot(c_scale(j), r, desc_a, info)
!Update residual !Update residual
call psb_geaxpby(zone, r, zzero, r, desc_a, info) call psb_geaxpby(zone, r, zzero, r, desc_a, info)
call psb_geaxpby(-alpha(j), c_scale(j), zone, r, desc_a, info) call psb_geaxpby(-alpha(j), c_scale(j), zone, r, desc_a, info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (j >= irst) exit iteration if (j >= irst) exit iteration
end do iteration end do iteration
m = j m = j
!Compute solution !Compute solution
call ztrsm('l','u','n','n',m,1,zone,h,size(h,1),alpha,size(alpha,1)) call ztrsm('l','u','n','n',m,1,zone,h,size(h,1),alpha,size(alpha,1))
if (nrst == 0 ) then if (nrst == 0 ) then
call x%set(zzero) call x%set(zzero)
endif endif
do i=1,m do i=1,m
call psb_geaxpby(alpha(i), z(i), zone, x, desc_a, info) call psb_geaxpby(alpha(i), z(i), zone, x, desc_a, info)
enddo enddo

Write Preview
Loading…
Cancel
Save