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

New code for FCG method.

Browse Source 
Header update for GMRES.
waitall
Salvatore Filippone 7 years ago
parent bf71bc6a04
commit eee0efe3b9
8 changed files with 320 additions and 316 deletions
Show Stats Download Patch File Download Diff File

151
krylov/psb_cfcg.F90
Unescape Escape View File

@ -121,15 +121,15 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(out) :: iter integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_spk_), Optional, Intent(out) :: err,cond real(psb_spk_), Optional, Intent(out) :: err,cond
! = Local data ! = Local data
type(psb_c_vect_type) :: v, w type(psb_c_vect_type) :: v, w, d , q, r
type(psb_c_vect_type), dimension(0:1) :: d complex(psb_spk_) :: alpha, beta, delta, gamma, theta
complex(psb_spk_) :: alpha, tau, tau1, beta, delta
real(psb_dpk_) :: derr real(psb_dpk_) :: derr
integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_ integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_
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
complex(psb_spk_), allocatable, target :: aux(:) complex(psb_spk_), allocatable, target :: aux(:)
complex(psb_spk_) :: vres(3)
character(len=20) :: name character(len=20) :: name
type(psb_itconv_type) :: stopdat type(psb_itconv_type) :: stopdat
character(len=*), parameter :: methdname='FCG' character(len=*), parameter :: methdname='FCG'
@ -191,112 +191,113 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
end if end if
!Assemble w, v !Assemble w, v, d, q, r, u
call psb_geasb(w, desc_a,info,&
call psb_geasb(w,&
& desc_a,info,&
& scratch=.true.,mold=b%v)
call psb_geasb(v,&
& desc_a,info,&
& scratch=.true.,mold=b%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& desc_a,info,&
& scratch=.true.,mold=x%v) & scratch=.true.,mold=x%v)
call psb_geasb(d(1),& call psb_geasb(v, desc_a,info,&
& desc_a,info,& & scratch=.true.,mold=x%v)
call psb_geasb(d, desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(q, desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(r, desc_a,info,&
& scratch=.true.,mold=x%v) & scratch=.true.,mold=x%v)
call psb_init_conv(methdname,istop_,itrace_,itmax_,&
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) & a,x,b,eps,desc_a,stopdat,info)
itx=0 itx = 0
restart: do restart: do
if (itx>= itmax_) exit restart if (itx>= itmax_) exit restart
! w=b ! r=b -Ax
call psb_geaxpby(cone,b,czero,w,& call psb_geaxpby(cone,b,czero,r, desc_a,info)
& desc_a,info) if (info == psb_success_) call psb_spmm(-cone,a,x,cone,r,desc_a,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (psb_errstatus_fatal()) then
nc2l = desc_a%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
!Compute v = Ax if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
call psb_spmm(cone,a,x,czero,v,desc_a,info) ! Apply the preconditioner v=Pr
! Compute w = Av
call prec%apply(r,v,desc_a,info,work=aux)
if (info == psb_success_) call psb_spmm(cone,a,v,czero,w,desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
!Compute w = -Ax + b
call psb_geaxpby(-cone, v, cone, w, desc_a, info) vres(1) = psb_gedot(r, v, desc_a, info, global = .false.)
vres(2) = psb_gedot(w, v, desc_a, info, global = .false.)
!Apply the preconditioner
idx=0 call psb_sum(ictxt, vres(1:2))
call prec%apply(w,d(idx),desc_a,info,work=aux) alpha = vres(1)
beta = vres(2)
delta = psb_gedot(d(idx), w, desc_a, info) ! d = v
call psb_geaxpby(cone, v, czero, d, desc_a, info)
! q = w
call psb_geaxpby(cone, w, czero, q, desc_a, info)
! compute delta=beta
! then
! x = x + (alpha/delta)*d
! r = r - (alpha/delta)*q
!Loop delta = beta
theta = alpha/delta
call psb_geaxpby(theta, d, cone, x, desc_a, info)
call psb_geaxpby(-theta, q, cone, r, desc_a, info)
if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart iteration: do
if (info /= psb_success_) Then itx = itx + 1
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
End If
iteration: do if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
call psb_spmm(cone,a,d(idx),czero,v,desc_a,info) ! Apply the preconditioner v = Pr
! Compute w = Av
call prec%apply(r,v,desc_a,info,work=aux)
if (info == psb_success_) call psb_spmm(cone,a,v,czero,w,desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residual'); goto 9999
goto 9999
end if end if
tau = psb_gedot(d(idx), v, desc_a, info)
alpha = delta/tau vres(1) = psb_gedot(r, v, desc_a, info, global = .false.)
!Update solution x vres(2) = psb_gedot(w, v, desc_a, info, global = .false.)
call psb_geaxpby(alpha, d(idx), cone, x, desc_a, info) vres(3) = psb_gedot(q, v, desc_a, info, global = .false.)
!Update residual w
call psb_geaxpby(-alpha, v, cone, w, desc_a, info)
itx = itx + 1
idx=mod(itx ,2)
call d(idx)%set(czero) call psb_sum(ictxt, vres(1:3))
call prec%apply(w,d(idx),desc_a,info,work=aux)
tau1= psb_gedot(d(idx), v, desc_a, info) alpha = vres(1)
beta=tau1/tau beta = vres(2)
gamma = vres(3)
if (idx == 1) then ! Compute d = v-(gamma/delta)*d
call psb_geaxpby(-beta, d(idx - 1), cone, d(idx), desc_a, info) ! q = w-(gamma/delta)*q
else theta= gamma/delta
call psb_geaxpby(-beta, d(idx + 1), cone, d(idx), desc_a, info) call psb_geaxpby(cone, v, -theta, d, desc_a, info)
endif call psb_geaxpby(cone, w, -theta, q , desc_a, info)
delta = psb_gedot(w, d(idx), desc_a, info) ! update delta
delta = beta - (gamma*gamma)/delta
if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart ! update u and r
if (info /= psb_success_) Then ! u = u + (alpha/delta)*d
call psb_errpush(psb_err_from_subroutine_non_,name) ! r = r - (alpha/delta)*q
goto 9999 theta= alpha/delta
End If call psb_geaxpby(theta, d, cone, x, desc_a, info)
call psb_geaxpby(-theta, q, cone, r, desc_a, info)
end do iteration end do iteration
end do restart end do restart

8
krylov/psb_crgmres.f90
Unescape Escape View File

@ -29,10 +29,10 @@
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
!!$ !
!!$ !
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C C ! C C
! C References: C ! C References: C

151
krylov/psb_dfcg.F90
Unescape Escape View File

@ -121,15 +121,15 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(out) :: iter integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_dpk_), Optional, Intent(out) :: err,cond real(psb_dpk_), Optional, Intent(out) :: err,cond
! = Local data ! = Local data
type(psb_d_vect_type) :: v, w type(psb_d_vect_type) :: v, w, d , q, r
type(psb_d_vect_type), dimension(0:1) :: d real(psb_dpk_) :: alpha, beta, delta, gamma, theta
real(psb_dpk_) :: alpha, tau, tau1, beta, delta
real(psb_dpk_) :: derr real(psb_dpk_) :: derr
integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_ integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_
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
real(psb_dpk_), allocatable, target :: aux(:) real(psb_dpk_), allocatable, target :: aux(:)
real(psb_dpk_) :: vres(3)
character(len=20) :: name character(len=20) :: name
type(psb_itconv_type) :: stopdat type(psb_itconv_type) :: stopdat
character(len=*), parameter :: methdname='FCG' character(len=*), parameter :: methdname='FCG'
@ -191,112 +191,113 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
end if end if
!Assemble w, v !Assemble w, v, d, q, r, u
call psb_geasb(w, desc_a,info,&
call psb_geasb(w,&
& desc_a,info,&
& scratch=.true.,mold=b%v)
call psb_geasb(v,&
& desc_a,info,&
& scratch=.true.,mold=b%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& desc_a,info,&
& scratch=.true.,mold=x%v) & scratch=.true.,mold=x%v)
call psb_geasb(d(1),& call psb_geasb(v, desc_a,info,&
& desc_a,info,& & scratch=.true.,mold=x%v)
call psb_geasb(d, desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(q, desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(r, desc_a,info,&
& scratch=.true.,mold=x%v) & scratch=.true.,mold=x%v)
call psb_init_conv(methdname,istop_,itrace_,itmax_,&
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) & a,x,b,eps,desc_a,stopdat,info)
itx=0 itx = 0
restart: do restart: do
if (itx>= itmax_) exit restart if (itx>= itmax_) exit restart
! w=b ! r=b -Ax
call psb_geaxpby(done,b,dzero,w,& call psb_geaxpby(done,b,dzero,r, desc_a,info)
& desc_a,info) if (info == psb_success_) call psb_spmm(-done,a,x,done,r,desc_a,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (psb_errstatus_fatal()) then
nc2l = desc_a%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
!Compute v = Ax if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
call psb_spmm(done,a,x,dzero,v,desc_a,info) ! Apply the preconditioner v=Pr
! Compute w = Av
call prec%apply(r,v,desc_a,info,work=aux)
if (info == psb_success_) call psb_spmm(done,a,v,dzero,w,desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
!Compute w = -Ax + b
call psb_geaxpby(-done, v, done, w, desc_a, info) vres(1) = psb_gedot(r, v, desc_a, info, global = .false.)
vres(2) = psb_gedot(w, v, desc_a, info, global = .false.)
!Apply the preconditioner
idx=0 call psb_sum(ictxt, vres(1:2))
call prec%apply(w,d(idx),desc_a,info,work=aux) alpha = vres(1)
beta = vres(2)
delta = psb_gedot(d(idx), w, desc_a, info) ! d = v
call psb_geaxpby(done, v, dzero, d, desc_a, info)
! q = w
call psb_geaxpby(done, w, dzero, q, desc_a, info)
! compute delta=beta
! then
! x = x + (alpha/delta)*d
! r = r - (alpha/delta)*q
!Loop delta = beta
theta = alpha/delta
call psb_geaxpby(theta, d, done, x, desc_a, info)
call psb_geaxpby(-theta, q, done, r, desc_a, info)
if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart iteration: do
if (info /= psb_success_) Then itx = itx + 1
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
End If
iteration: do if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
call psb_spmm(done,a,d(idx),dzero,v,desc_a,info) ! Apply the preconditioner v = Pr
! Compute w = Av
call prec%apply(r,v,desc_a,info,work=aux)
if (info == psb_success_) call psb_spmm(done,a,v,dzero,w,desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residual'); goto 9999
goto 9999
end if end if
tau = psb_gedot(d(idx), v, desc_a, info)
alpha = delta/tau vres(1) = psb_gedot(r, v, desc_a, info, global = .false.)
!Update solution x vres(2) = psb_gedot(w, v, desc_a, info, global = .false.)
call psb_geaxpby(alpha, d(idx), done, x, desc_a, info) vres(3) = psb_gedot(q, v, desc_a, info, global = .false.)
!Update residual w
call psb_geaxpby(-alpha, v, done, w, desc_a, info)
itx = itx + 1
idx=mod(itx ,2)
call d(idx)%set(dzero) call psb_sum(ictxt, vres(1:3))
call prec%apply(w,d(idx),desc_a,info,work=aux)
tau1= psb_gedot(d(idx), v, desc_a, info) alpha = vres(1)
beta=tau1/tau beta = vres(2)
gamma = vres(3)
if (idx == 1) then ! Compute d = v-(gamma/delta)*d
call psb_geaxpby(-beta, d(idx - 1), done, d(idx), desc_a, info) ! q = w-(gamma/delta)*q
else theta= gamma/delta
call psb_geaxpby(-beta, d(idx + 1), done, d(idx), desc_a, info) call psb_geaxpby(done, v, -theta, d, desc_a, info)
endif call psb_geaxpby(done, w, -theta, q , desc_a, info)
delta = psb_gedot(w, d(idx), desc_a, info) ! update delta
delta = beta - (gamma*gamma)/delta
if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart ! update u and r
if (info /= psb_success_) Then ! u = u + (alpha/delta)*d
call psb_errpush(psb_err_from_subroutine_non_,name) ! r = r - (alpha/delta)*q
goto 9999 theta= alpha/delta
End If call psb_geaxpby(theta, d, done, x, desc_a, info)
call psb_geaxpby(-theta, q, done, r, desc_a, info)
end do iteration end do iteration
end do restart end do restart

8
krylov/psb_drgmres.f90
Unescape Escape View File

@ -29,10 +29,10 @@
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
!!$ !
!!$ !
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C C ! C C
! C References: C ! C References: C

151
krylov/psb_sfcg.F90
Unescape Escape View File

@ -121,15 +121,15 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(out) :: iter integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_spk_), Optional, Intent(out) :: err,cond real(psb_spk_), Optional, Intent(out) :: err,cond
! = Local data ! = Local data
type(psb_s_vect_type) :: v, w type(psb_s_vect_type) :: v, w, d , q, r
type(psb_s_vect_type), dimension(0:1) :: d real(psb_spk_) :: alpha, beta, delta, gamma, theta
real(psb_spk_) :: alpha, tau, tau1, beta, delta
real(psb_dpk_) :: derr real(psb_dpk_) :: derr
integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_ integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_
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
real(psb_spk_), allocatable, target :: aux(:) real(psb_spk_), allocatable, target :: aux(:)
real(psb_spk_) :: vres(3)
character(len=20) :: name character(len=20) :: name
type(psb_itconv_type) :: stopdat type(psb_itconv_type) :: stopdat
character(len=*), parameter :: methdname='FCG' character(len=*), parameter :: methdname='FCG'
@ -191,112 +191,113 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
end if end if
!Assemble w, v !Assemble w, v, d, q, r, u
call psb_geasb(w, desc_a,info,&
call psb_geasb(w,&
& desc_a,info,&
& scratch=.true.,mold=b%v)
call psb_geasb(v,&
& desc_a,info,&
& scratch=.true.,mold=b%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& desc_a,info,&
& scratch=.true.,mold=x%v) & scratch=.true.,mold=x%v)
call psb_geasb(d(1),& call psb_geasb(v, desc_a,info,&
& desc_a,info,& & scratch=.true.,mold=x%v)
call psb_geasb(d, desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(q, desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(r, desc_a,info,&
& scratch=.true.,mold=x%v) & scratch=.true.,mold=x%v)
call psb_init_conv(methdname,istop_,itrace_,itmax_,&
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) & a,x,b,eps,desc_a,stopdat,info)
itx=0 itx = 0
restart: do restart: do
if (itx>= itmax_) exit restart if (itx>= itmax_) exit restart
! w=b ! r=b -Ax
call psb_geaxpby(sone,b,szero,w,& call psb_geaxpby(sone,b,szero,r, desc_a,info)
& desc_a,info) if (info == psb_success_) call psb_spmm(-sone,a,x,sone,r,desc_a,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (psb_errstatus_fatal()) then
nc2l = desc_a%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
!Compute v = Ax if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
call psb_spmm(sone,a,x,szero,v,desc_a,info) ! Apply the preconditioner v=Pr
! Compute w = Av
call prec%apply(r,v,desc_a,info,work=aux)
if (info == psb_success_) call psb_spmm(sone,a,v,szero,w,desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
!Compute w = -Ax + b
call psb_geaxpby(-sone, v, sone, w, desc_a, info) vres(1) = psb_gedot(r, v, desc_a, info, global = .false.)
vres(2) = psb_gedot(w, v, desc_a, info, global = .false.)
!Apply the preconditioner
idx=0 call psb_sum(ictxt, vres(1:2))
call prec%apply(w,d(idx),desc_a,info,work=aux) alpha = vres(1)
beta = vres(2)
delta = psb_gedot(d(idx), w, desc_a, info) ! d = v
call psb_geaxpby(sone, v, szero, d, desc_a, info)
! q = w
call psb_geaxpby(sone, w, szero, q, desc_a, info)
! compute delta=beta
! then
! x = x + (alpha/delta)*d
! r = r - (alpha/delta)*q
!Loop delta = beta
theta = alpha/delta
call psb_geaxpby(theta, d, sone, x, desc_a, info)
call psb_geaxpby(-theta, q, sone, r, desc_a, info)
if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart iteration: do
if (info /= psb_success_) Then itx = itx + 1
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
End If
iteration: do if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
call psb_spmm(sone,a,d(idx),szero,v,desc_a,info) ! Apply the preconditioner v = Pr
! Compute w = Av
call prec%apply(r,v,desc_a,info,work=aux)
if (info == psb_success_) call psb_spmm(sone,a,v,szero,w,desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residual'); goto 9999
goto 9999
end if end if
tau = psb_gedot(d(idx), v, desc_a, info)
alpha = delta/tau vres(1) = psb_gedot(r, v, desc_a, info, global = .false.)
!Update solution x vres(2) = psb_gedot(w, v, desc_a, info, global = .false.)
call psb_geaxpby(alpha, d(idx), sone, x, desc_a, info) vres(3) = psb_gedot(q, v, desc_a, info, global = .false.)
!Update residual w
call psb_geaxpby(-alpha, v, sone, w, desc_a, info)
itx = itx + 1
idx=mod(itx ,2)
call d(idx)%set(szero) call psb_sum(ictxt, vres(1:3))
call prec%apply(w,d(idx),desc_a,info,work=aux)
tau1= psb_gedot(d(idx), v, desc_a, info) alpha = vres(1)
beta=tau1/tau beta = vres(2)
gamma = vres(3)
if (idx == 1) then ! Compute d = v-(gamma/delta)*d
call psb_geaxpby(-beta, d(idx - 1), sone, d(idx), desc_a, info) ! q = w-(gamma/delta)*q
else theta= gamma/delta
call psb_geaxpby(-beta, d(idx + 1), sone, d(idx), desc_a, info) call psb_geaxpby(sone, v, -theta, d, desc_a, info)
endif call psb_geaxpby(sone, w, -theta, q , desc_a, info)
delta = psb_gedot(w, d(idx), desc_a, info) ! update delta
delta = beta - (gamma*gamma)/delta
if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart ! update u and r
if (info /= psb_success_) Then ! u = u + (alpha/delta)*d
call psb_errpush(psb_err_from_subroutine_non_,name) ! r = r - (alpha/delta)*q
goto 9999 theta= alpha/delta
End If call psb_geaxpby(theta, d, sone, x, desc_a, info)
call psb_geaxpby(-theta, q, sone, r, desc_a, info)
end do iteration end do iteration
end do restart end do restart

8
krylov/psb_srgmres.f90
Unescape Escape View File

@ -29,10 +29,10 @@
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
!!$ !
!!$ !
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C C ! C C
! C References: C ! C References: C

151
krylov/psb_zfcg.F90
Unescape Escape View File

@ -121,15 +121,15 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(out) :: iter integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_dpk_), Optional, Intent(out) :: err,cond real(psb_dpk_), Optional, Intent(out) :: err,cond
! = Local data ! = Local data
type(psb_z_vect_type) :: v, w type(psb_z_vect_type) :: v, w, d , q, r
type(psb_z_vect_type), dimension(0:1) :: d complex(psb_dpk_) :: alpha, beta, delta, gamma, theta
complex(psb_dpk_) :: alpha, tau, tau1, beta, delta
real(psb_dpk_) :: derr real(psb_dpk_) :: derr
integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_ integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_
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
complex(psb_dpk_), allocatable, target :: aux(:) complex(psb_dpk_), allocatable, target :: aux(:)
complex(psb_dpk_) :: vres(3)
character(len=20) :: name character(len=20) :: name
type(psb_itconv_type) :: stopdat type(psb_itconv_type) :: stopdat
character(len=*), parameter :: methdname='FCG' character(len=*), parameter :: methdname='FCG'
@ -191,112 +191,113 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
end if end if
!Assemble w, v !Assemble w, v, d, q, r, u
call psb_geasb(w, desc_a,info,&
call psb_geasb(w,&
& desc_a,info,&
& scratch=.true.,mold=b%v)
call psb_geasb(v,&
& desc_a,info,&
& scratch=.true.,mold=b%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& desc_a,info,&
& scratch=.true.,mold=x%v) & scratch=.true.,mold=x%v)
call psb_geasb(d(1),& call psb_geasb(v, desc_a,info,&
& desc_a,info,& & scratch=.true.,mold=x%v)
call psb_geasb(d, desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(q, desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(r, desc_a,info,&
& scratch=.true.,mold=x%v) & scratch=.true.,mold=x%v)
call psb_init_conv(methdname,istop_,itrace_,itmax_,&
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) & a,x,b,eps,desc_a,stopdat,info)
itx=0 itx = 0
restart: do restart: do
if (itx>= itmax_) exit restart if (itx>= itmax_) exit restart
! w=b ! r=b -Ax
call psb_geaxpby(zone,b,zzero,w,& call psb_geaxpby(zone,b,zzero,r, desc_a,info)
& desc_a,info) if (info == psb_success_) call psb_spmm(-zone,a,x,zone,r,desc_a,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (psb_errstatus_fatal()) then
nc2l = desc_a%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
!Compute v = Ax if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
call psb_spmm(zone,a,x,zzero,v,desc_a,info) ! Apply the preconditioner v=Pr
! Compute w = Av
call prec%apply(r,v,desc_a,info,work=aux)
if (info == psb_success_) call psb_spmm(zone,a,v,zzero,w,desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
!Compute w = -Ax + b
call psb_geaxpby(-zone, v, zone, w, desc_a, info) vres(1) = psb_gedot(r, v, desc_a, info, global = .false.)
vres(2) = psb_gedot(w, v, desc_a, info, global = .false.)
!Apply the preconditioner
idx=0 call psb_sum(ictxt, vres(1:2))
call prec%apply(w,d(idx),desc_a,info,work=aux) alpha = vres(1)
beta = vres(2)
delta = psb_gedot(d(idx), w, desc_a, info) ! d = v
call psb_geaxpby(zone, v, zzero, d, desc_a, info)
! q = w
call psb_geaxpby(zone, w, zzero, q, desc_a, info)
! compute delta=beta
! then
! x = x + (alpha/delta)*d
! r = r - (alpha/delta)*q
!Loop delta = beta
theta = alpha/delta
call psb_geaxpby(theta, d, zone, x, desc_a, info)
call psb_geaxpby(-theta, q, zone, r, desc_a, info)
if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart iteration: do
if (info /= psb_success_) Then itx = itx + 1
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
End If
iteration: do if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
call psb_spmm(zone,a,d(idx),zzero,v,desc_a,info) ! Apply the preconditioner v = Pr
! Compute w = Av
call prec%apply(r,v,desc_a,info,work=aux)
if (info == psb_success_) call psb_spmm(zone,a,v,zzero,w,desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residual'); goto 9999
goto 9999
end if end if
tau = psb_gedot(d(idx), v, desc_a, info)
alpha = delta/tau vres(1) = psb_gedot(r, v, desc_a, info, global = .false.)
!Update solution x vres(2) = psb_gedot(w, v, desc_a, info, global = .false.)
call psb_geaxpby(alpha, d(idx), zone, x, desc_a, info) vres(3) = psb_gedot(q, v, desc_a, info, global = .false.)
!Update residual w
call psb_geaxpby(-alpha, v, zone, w, desc_a, info)
itx = itx + 1
idx=mod(itx ,2)
call d(idx)%set(zzero) call psb_sum(ictxt, vres(1:3))
call prec%apply(w,d(idx),desc_a,info,work=aux)
tau1= psb_gedot(d(idx), v, desc_a, info) alpha = vres(1)
beta=tau1/tau beta = vres(2)
gamma = vres(3)
if (idx == 1) then ! Compute d = v-(gamma/delta)*d
call psb_geaxpby(-beta, d(idx - 1), zone, d(idx), desc_a, info) ! q = w-(gamma/delta)*q
else theta= gamma/delta
call psb_geaxpby(-beta, d(idx + 1), zone, d(idx), desc_a, info) call psb_geaxpby(zone, v, -theta, d, desc_a, info)
endif call psb_geaxpby(zone, w, -theta, q , desc_a, info)
delta = psb_gedot(w, d(idx), desc_a, info) ! update delta
delta = beta - (gamma*gamma)/delta
if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart ! update u and r
if (info /= psb_success_) Then ! u = u + (alpha/delta)*d
call psb_errpush(psb_err_from_subroutine_non_,name) ! r = r - (alpha/delta)*q
goto 9999 theta= alpha/delta
End If call psb_geaxpby(theta, d, zone, x, desc_a, info)
call psb_geaxpby(-theta, q, zone, r, desc_a, info)
end do iteration end do iteration
end do restart end do restart

8
krylov/psb_zrgmres.f90
Unescape Escape View File

@ -29,10 +29,10 @@
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
!!$ !
!!$ !
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C C ! C C
! C References: C ! C References: C

Write Preview
Loading…
Cancel
Save