From eee0efe3b9732503d02aef2ec5269631f01b12c1 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Thu, 14 Jun 2018 21:37:36 +0100 Subject: [PATCH] New code for FCG method. Header update for GMRES. --- krylov/psb_cfcg.F90 | 151 +++++++++++++++++++++-------------------- krylov/psb_crgmres.f90 | 8 +-- krylov/psb_dfcg.F90 | 151 +++++++++++++++++++++-------------------- krylov/psb_drgmres.f90 | 8 +-- krylov/psb_sfcg.F90 | 151 +++++++++++++++++++++-------------------- krylov/psb_srgmres.f90 | 8 +-- krylov/psb_zfcg.F90 | 151 +++++++++++++++++++++-------------------- krylov/psb_zrgmres.f90 | 8 +-- 8 files changed, 320 insertions(+), 316 deletions(-) diff --git a/krylov/psb_cfcg.F90 b/krylov/psb_cfcg.F90 index 0e4c0271..c042571b 100644 --- a/krylov/psb_cfcg.F90 +++ b/krylov/psb_cfcg.F90 @@ -121,15 +121,15 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_spk_), Optional, Intent(out) :: err,cond ! = Local data - type(psb_c_vect_type) :: v, w - type(psb_c_vect_type), dimension(0:1) :: d - complex(psb_spk_) :: alpha, tau, tau1, beta, delta + type(psb_c_vect_type) :: v, w, d , q, r + complex(psb_spk_) :: alpha, beta, delta, gamma, theta real(psb_dpk_) :: derr integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_ integer(psb_ipk_) :: n_col, mglob, naux, err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: np, me, ictxt complex(psb_spk_), allocatable, target :: aux(:) + complex(psb_spk_) :: vres(3) character(len=20) :: name type(psb_itconv_type) :: stopdat character(len=*), parameter :: methdname='FCG' @@ -191,112 +191,113 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,& end if - !Assemble w, v - - 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,& + !Assemble w, v, d, q, r, u + call psb_geasb(w, desc_a,info,& & scratch=.true.,mold=x%v) - call psb_geasb(d(1),& - & desc_a,info,& + call psb_geasb(v, 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) - - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) - itx=0 + call psb_init_conv(methdname,istop_,itrace_,itmax_,& + & a,x,b,eps,desc_a,stopdat,info) + itx = 0 restart: do if (itx>= itmax_) exit restart - ! w=b - call psb_geaxpby(cone,b,czero,w,& - & desc_a,info) + ! r=b -Ax + call psb_geaxpby(cone,b,czero,r, 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 + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart - !Compute v = Ax - 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 call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') + & a_err='Error during residue') goto 9999 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 - call psb_errpush(psb_err_from_subroutine_non_,name) - goto 9999 - End If + itx = itx + 1 - 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 - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residual'); goto 9999 end if - tau = psb_gedot(d(idx), v, desc_a, info) - - alpha = delta/tau - !Update solution x - call psb_geaxpby(alpha, d(idx), cone, x, desc_a, info) - !Update residual w - call psb_geaxpby(-alpha, 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.) + vres(3) = psb_gedot(q, v, desc_a, info, global = .false.) - itx = itx + 1 - idx=mod(itx ,2) + call psb_sum(ictxt, vres(1:3)) - call d(idx)%set(czero) - call prec%apply(w,d(idx),desc_a,info,work=aux) + alpha = vres(1) + beta = vres(2) + gamma = vres(3) - tau1= psb_gedot(d(idx), v, desc_a, info) - beta=tau1/tau + ! Compute d = v-(gamma/delta)*d + ! q = w-(gamma/delta)*q + theta= gamma/delta + call psb_geaxpby(cone, v, -theta, d, desc_a, info) + call psb_geaxpby(cone, w, -theta, q , desc_a, info) - if (idx == 1) then - call psb_geaxpby(-beta, d(idx - 1), cone, d(idx), desc_a, info) - else - call psb_geaxpby(-beta, d(idx + 1), cone, d(idx), desc_a, info) - endif - - 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 - if (info /= psb_success_) Then - call psb_errpush(psb_err_from_subroutine_non_,name) - goto 9999 - End If + ! update u and r + ! u = u + (alpha/delta)*d + ! r = r - (alpha/delta)*q + theta= alpha/delta + 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 restart diff --git a/krylov/psb_crgmres.f90 b/krylov/psb_crgmres.f90 index d7e75251..f68dc74b 100644 --- a/krylov/psb_crgmres.f90 +++ b/krylov/psb_crgmres.f90 @@ -29,10 +29,10 @@ ! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN ! 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 -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! C C ! C References: C diff --git a/krylov/psb_dfcg.F90 b/krylov/psb_dfcg.F90 index 2821829a..90e5b116 100644 --- a/krylov/psb_dfcg.F90 +++ b/krylov/psb_dfcg.F90 @@ -121,15 +121,15 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_dpk_), Optional, Intent(out) :: err,cond ! = Local data - type(psb_d_vect_type) :: v, w - type(psb_d_vect_type), dimension(0:1) :: d - real(psb_dpk_) :: alpha, tau, tau1, beta, delta + type(psb_d_vect_type) :: v, w, d , q, r + real(psb_dpk_) :: alpha, beta, delta, gamma, theta real(psb_dpk_) :: derr integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_ integer(psb_ipk_) :: n_col, mglob, naux, err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: np, me, ictxt real(psb_dpk_), allocatable, target :: aux(:) + real(psb_dpk_) :: vres(3) character(len=20) :: name type(psb_itconv_type) :: stopdat character(len=*), parameter :: methdname='FCG' @@ -191,112 +191,113 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,& end if - !Assemble w, v - - 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,& + !Assemble w, v, d, q, r, u + call psb_geasb(w, desc_a,info,& & scratch=.true.,mold=x%v) - call psb_geasb(d(1),& - & desc_a,info,& + call psb_geasb(v, 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) - - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) - itx=0 + call psb_init_conv(methdname,istop_,itrace_,itmax_,& + & a,x,b,eps,desc_a,stopdat,info) + itx = 0 restart: do if (itx>= itmax_) exit restart - ! w=b - call psb_geaxpby(done,b,dzero,w,& - & desc_a,info) + ! r=b -Ax + call psb_geaxpby(done,b,dzero,r, 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 + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart - !Compute v = Ax - 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 call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') + & a_err='Error during residue') goto 9999 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 - call psb_errpush(psb_err_from_subroutine_non_,name) - goto 9999 - End If + itx = itx + 1 - 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 - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residual'); goto 9999 end if - tau = psb_gedot(d(idx), v, desc_a, info) - - alpha = delta/tau - !Update solution x - call psb_geaxpby(alpha, d(idx), done, x, desc_a, info) - !Update residual w - call psb_geaxpby(-alpha, 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.) + vres(3) = psb_gedot(q, v, desc_a, info, global = .false.) - itx = itx + 1 - idx=mod(itx ,2) + call psb_sum(ictxt, vres(1:3)) - call d(idx)%set(dzero) - call prec%apply(w,d(idx),desc_a,info,work=aux) + alpha = vres(1) + beta = vres(2) + gamma = vres(3) - tau1= psb_gedot(d(idx), v, desc_a, info) - beta=tau1/tau + ! Compute d = v-(gamma/delta)*d + ! q = w-(gamma/delta)*q + theta= gamma/delta + call psb_geaxpby(done, v, -theta, d, desc_a, info) + call psb_geaxpby(done, w, -theta, q , desc_a, info) - if (idx == 1) then - call psb_geaxpby(-beta, d(idx - 1), done, d(idx), desc_a, info) - else - call psb_geaxpby(-beta, d(idx + 1), done, d(idx), desc_a, info) - endif - - 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 - if (info /= psb_success_) Then - call psb_errpush(psb_err_from_subroutine_non_,name) - goto 9999 - End If + ! update u and r + ! u = u + (alpha/delta)*d + ! r = r - (alpha/delta)*q + theta= alpha/delta + 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 restart diff --git a/krylov/psb_drgmres.f90 b/krylov/psb_drgmres.f90 index 368a6482..9b577903 100644 --- a/krylov/psb_drgmres.f90 +++ b/krylov/psb_drgmres.f90 @@ -29,10 +29,10 @@ ! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN ! 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 -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! C C ! C References: C diff --git a/krylov/psb_sfcg.F90 b/krylov/psb_sfcg.F90 index 2fae24b4..1f074371 100644 --- a/krylov/psb_sfcg.F90 +++ b/krylov/psb_sfcg.F90 @@ -121,15 +121,15 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_spk_), Optional, Intent(out) :: err,cond ! = Local data - type(psb_s_vect_type) :: v, w - type(psb_s_vect_type), dimension(0:1) :: d - real(psb_spk_) :: alpha, tau, tau1, beta, delta + type(psb_s_vect_type) :: v, w, d , q, r + real(psb_spk_) :: alpha, beta, delta, gamma, theta real(psb_dpk_) :: derr integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_ integer(psb_ipk_) :: n_col, mglob, naux, err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: np, me, ictxt real(psb_spk_), allocatable, target :: aux(:) + real(psb_spk_) :: vres(3) character(len=20) :: name type(psb_itconv_type) :: stopdat character(len=*), parameter :: methdname='FCG' @@ -191,112 +191,113 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,& end if - !Assemble w, v - - 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,& + !Assemble w, v, d, q, r, u + call psb_geasb(w, desc_a,info,& & scratch=.true.,mold=x%v) - call psb_geasb(d(1),& - & desc_a,info,& + call psb_geasb(v, 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) - - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) - itx=0 + call psb_init_conv(methdname,istop_,itrace_,itmax_,& + & a,x,b,eps,desc_a,stopdat,info) + itx = 0 restart: do if (itx>= itmax_) exit restart - ! w=b - call psb_geaxpby(sone,b,szero,w,& - & desc_a,info) + ! r=b -Ax + call psb_geaxpby(sone,b,szero,r, 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 + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart - !Compute v = Ax - 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 call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') + & a_err='Error during residue') goto 9999 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 - call psb_errpush(psb_err_from_subroutine_non_,name) - goto 9999 - End If + itx = itx + 1 - 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 - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residual'); goto 9999 end if - tau = psb_gedot(d(idx), v, desc_a, info) - - alpha = delta/tau - !Update solution x - call psb_geaxpby(alpha, d(idx), sone, x, desc_a, info) - !Update residual w - call psb_geaxpby(-alpha, 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.) + vres(3) = psb_gedot(q, v, desc_a, info, global = .false.) - itx = itx + 1 - idx=mod(itx ,2) + call psb_sum(ictxt, vres(1:3)) - call d(idx)%set(szero) - call prec%apply(w,d(idx),desc_a,info,work=aux) + alpha = vres(1) + beta = vres(2) + gamma = vres(3) - tau1= psb_gedot(d(idx), v, desc_a, info) - beta=tau1/tau + ! Compute d = v-(gamma/delta)*d + ! q = w-(gamma/delta)*q + theta= gamma/delta + call psb_geaxpby(sone, v, -theta, d, desc_a, info) + call psb_geaxpby(sone, w, -theta, q , desc_a, info) - if (idx == 1) then - call psb_geaxpby(-beta, d(idx - 1), sone, d(idx), desc_a, info) - else - call psb_geaxpby(-beta, d(idx + 1), sone, d(idx), desc_a, info) - endif - - 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 - if (info /= psb_success_) Then - call psb_errpush(psb_err_from_subroutine_non_,name) - goto 9999 - End If + ! update u and r + ! u = u + (alpha/delta)*d + ! r = r - (alpha/delta)*q + theta= alpha/delta + 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 restart diff --git a/krylov/psb_srgmres.f90 b/krylov/psb_srgmres.f90 index 311bc0fb..2e458d82 100644 --- a/krylov/psb_srgmres.f90 +++ b/krylov/psb_srgmres.f90 @@ -29,10 +29,10 @@ ! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN ! 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 -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! C C ! C References: C diff --git a/krylov/psb_zfcg.F90 b/krylov/psb_zfcg.F90 index 0e2ad2da..1e180685 100644 --- a/krylov/psb_zfcg.F90 +++ b/krylov/psb_zfcg.F90 @@ -121,15 +121,15 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,& integer(psb_ipk_), Optional, Intent(out) :: iter real(psb_dpk_), Optional, Intent(out) :: err,cond ! = Local data - type(psb_z_vect_type) :: v, w - type(psb_z_vect_type), dimension(0:1) :: d - complex(psb_dpk_) :: alpha, tau, tau1, beta, delta + type(psb_z_vect_type) :: v, w, d , q, r + complex(psb_dpk_) :: alpha, beta, delta, gamma, theta real(psb_dpk_) :: derr integer(psb_ipk_) :: i, idx, nc2l, it, itx, istop_, itmax_, itrace_ integer(psb_ipk_) :: n_col, mglob, naux, err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: np, me, ictxt complex(psb_dpk_), allocatable, target :: aux(:) + complex(psb_dpk_) :: vres(3) character(len=20) :: name type(psb_itconv_type) :: stopdat character(len=*), parameter :: methdname='FCG' @@ -191,112 +191,113 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,& end if - !Assemble w, v - - 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,& + !Assemble w, v, d, q, r, u + call psb_geasb(w, desc_a,info,& & scratch=.true.,mold=x%v) - call psb_geasb(d(1),& - & desc_a,info,& + call psb_geasb(v, 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) - - call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) - itx=0 + call psb_init_conv(methdname,istop_,itrace_,itmax_,& + & a,x,b,eps,desc_a,stopdat,info) + itx = 0 restart: do if (itx>= itmax_) exit restart - ! w=b - call psb_geaxpby(zone,b,zzero,w,& - & desc_a,info) + ! r=b -Ax + call psb_geaxpby(zone,b,zzero,r, 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 + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart - !Compute v = Ax - 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 call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') + & a_err='Error during residue') goto 9999 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 - call psb_errpush(psb_err_from_subroutine_non_,name) - goto 9999 - End If + itx = itx + 1 - 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 - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residual'); goto 9999 end if - tau = psb_gedot(d(idx), v, desc_a, info) - - alpha = delta/tau - !Update solution x - call psb_geaxpby(alpha, d(idx), zone, x, desc_a, info) - !Update residual w - call psb_geaxpby(-alpha, 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.) + vres(3) = psb_gedot(q, v, desc_a, info, global = .false.) - itx = itx + 1 - idx=mod(itx ,2) + call psb_sum(ictxt, vres(1:3)) - call d(idx)%set(zzero) - call prec%apply(w,d(idx),desc_a,info,work=aux) + alpha = vres(1) + beta = vres(2) + gamma = vres(3) - tau1= psb_gedot(d(idx), v, desc_a, info) - beta=tau1/tau + ! Compute d = v-(gamma/delta)*d + ! q = w-(gamma/delta)*q + theta= gamma/delta + call psb_geaxpby(zone, v, -theta, d, desc_a, info) + call psb_geaxpby(zone, w, -theta, q , desc_a, info) - if (idx == 1) then - call psb_geaxpby(-beta, d(idx - 1), zone, d(idx), desc_a, info) - else - call psb_geaxpby(-beta, d(idx + 1), zone, d(idx), desc_a, info) - endif - - 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 - if (info /= psb_success_) Then - call psb_errpush(psb_err_from_subroutine_non_,name) - goto 9999 - End If + ! update u and r + ! u = u + (alpha/delta)*d + ! r = r - (alpha/delta)*q + theta= alpha/delta + 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 restart diff --git a/krylov/psb_zrgmres.f90 b/krylov/psb_zrgmres.f90 index deb4111f..10738699 100644 --- a/krylov/psb_zrgmres.f90 +++ b/krylov/psb_zrgmres.f90 @@ -29,10 +29,10 @@ ! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN ! 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 -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! C C ! C References: C