From c83985fe1661f0cde23b12622ca702c8cba007f0 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 23 May 2016 15:03:44 +0000 Subject: [PATCH] psblas: krylov/psb_ccgr.f90 krylov/psb_dcgr.f90 krylov/psb_scgr.f90 krylov/psb_zcgr.f90 Cosmetic changes --- krylov/psb_ccgr.f90 | 208 +++++++++++++++++++++----------------------- krylov/psb_dcgr.f90 | 208 +++++++++++++++++++++----------------------- krylov/psb_scgr.f90 | 208 +++++++++++++++++++++----------------------- krylov/psb_zcgr.f90 | 208 +++++++++++++++++++++----------------------- 4 files changed, 396 insertions(+), 436 deletions(-) diff --git a/krylov/psb_ccgr.f90 b/krylov/psb_ccgr.f90 index 28067bab..9a5531c8 100644 --- a/krylov/psb_ccgr.f90 +++ b/krylov/psb_ccgr.f90 @@ -112,8 +112,8 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,& use psb_c_krylov_conv_mod use psb_krylov_mod implicit none - - + + type(psb_cspmat_type), intent(in) :: a Type(psb_desc_type), Intent(in) :: desc_a class(psb_cprec_type), intent(inout) :: prec @@ -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(out) :: iter real(psb_spk_), Optional, Intent(out) :: err -! = local data + ! = local data complex(psb_spk_), allocatable :: alpha(:), h(:,:) type(psb_c_vect_type), allocatable :: z(:), c(:), c_scale(:) 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_) :: debug_level, debug_unit integer(psb_ipk_) :: np, me, ictxt @@ -140,14 +140,14 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,& type(psb_itconv_type) :: stopdat character(len=*), parameter :: methdname='CGR' integer(psb_ipk_) ::int_err(5) - info = psb_success_ + info = psb_success_ name = 'psb_ccgr' call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - + ictxt = desc_a%get_context() - + call psb_info(ictxt, me, np) if (.not.allocated(b%v)) then info = psb_err_invalid_vect_state_ @@ -159,21 +159,21 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 endif - + mglob = desc_a%get_global_rows() n_col = desc_a%get_local_cols() - + if (present(istop)) then istop_ = istop else istop_ = 2 endif - -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - + + ! + ! ISTOP_ = 1: Normwise backward error, infinity norm + ! ISTOP_ = 2: ||r||/||b||, 2-norm + ! + if ((istop_ < 1 ).or.(istop_ > 2 ) ) then info=psb_err_invalid_istop_ int_err(1)=istop_ @@ -181,8 +181,8 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,i_err=int_err) goto 9999 endif - - + + call psb_chkvect(mglob,ione,x%get_nrows(),ione,ione,desc_a,info) if (info == psb_success_)& & call psb_chkvect(mglob,ione,b%get_nrows(),ione,ione,desc_a,info) @@ -191,19 +191,19 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,a_err='psb_chkvect on X/B') goto 9999 end if - + if (present(itmax)) then itmax_ = itmax else itmax_ = 1000 endif - + if (present(itrace)) then itrace_ = itrace else itrace_ = 0 end if - + if (present(irst)) then nl = irst if (debug_level >= psb_debug_ext_) & @@ -215,7 +215,7 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,& & write(debug_unit,*) me,' ',trim(name),& & ' not present: irst: ',irst,nl endif - + if (nl <=0 ) then info=psb_err_invalid_istop_ int_err(1)=nl @@ -223,54 +223,44 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,i_err=int_err) goto 9999 endif - + naux=4*n_col allocate(aux(naux),h(nl+1,nl+1),& &c_scale(nl+1),c(nl+1),z(nl+1), alpha(nl+1), stat=info) - + h = czero if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - - 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 - - itx = 0 - - if (istop_ == 2) then - b_norm = psb_norm2(b, desc_a, 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) + + call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v) + + 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 + + itx = 0 + + nrst = -1 + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) restart: do if (itx>= itmax_) exit restart h = czero - + it = 0 ! compute r0 = b-ax0 - + if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - - + + call psb_geaxpby(cone, b, czero, r, desc_a, info) call psb_spmm(-cone,a,x,cone,r,desc_a,info,work=aux) if (info /= psb_success_) then @@ -278,65 +268,65 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 end if - + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) ! if (info /= psb_success_) Then ! call psb_errpush(psb_err_from_subroutine_non_,name) ! goto 9999 ! End If - - nrst = nrst + 1 - + + nrst = nrst + 1 + iteration: do - - itx = itx + 1 - it = it + 1 - j = it - !Apply preconditioner - 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) - do i =1, j - 1 - + + itx = itx + 1 + it = it + 1 + j = it + !Apply preconditioner + 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) + do i =1, j - 1 + 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(-h(i,j), c_scale(i), cone, c(i+1), desc_a, info) - end do - - h(j,j) = psb_norm2(c(j), desc_a, info) - hjj = cone/h(j,j) - call psb_geaxpby(hjj, c(j), czero, c_scale(j), desc_a, info) - - alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) - - !Update residual - call psb_geaxpby(cone, r, czero, 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 (j >= irst) exit iteration - - - end do iteration - - m = j - - !Compute solution - - call ctrsm('l','u','n','n',m,1,cone,h,size(h,1),alpha,size(alpha,1)) - - if (nrst == 0 ) then - call x%set(czero) - endif - do i=1,m - call psb_geaxpby(alpha(i), z(i), cone, x, desc_a, info) - enddo - - - - + end do + + h(j,j) = psb_norm2(c(j), desc_a, info) + hjj = cone/h(j,j) + call psb_geaxpby(hjj, c(j), czero, c_scale(j), desc_a, info) + + alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) + + !Update residual + call psb_geaxpby(cone, r, czero, 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 (j >= irst) exit iteration + + + end do iteration + + m = j + + !Compute solution + + call ctrsm('l','u','n','n',m,1,cone,h,size(h,1),alpha,size(alpha,1)) + + if (nrst == 0 ) then + call x%set(czero) + endif + do i=1,m + call psb_geaxpby(alpha(i), z(i), cone, x, desc_a, info) + enddo + + + + end do restart m = j !Compute solution @@ -347,38 +337,38 @@ subroutine psb_ccgr_vect(a,prec,b,x,eps,desc_a,info,& enddo iter = j - + call psb_end_conv(methdname,itx,desc_a,stopdat,info,derr,iter) if (present(err)) err = derr - + if (info == psb_success_) call psb_gefree(r,desc_a,info) - + do j = 1,m if (info == psb_success_) call psb_gefree(z(j),desc_a,info) if (info == psb_success_) call psb_gefree(c_scale(j),desc_a,info) enddo - + do i =1,nl+1 if (info == psb_success_) call psb_gefree(c(i),desc_a,info) end do - + if (info == psb_success_) deallocate(aux,h,c_scale,z,c,alpha,stat=info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 continue call psb_erractionrestore(err_act) if (err_act.eq.psb_act_abort_) then call psb_error() return end if - + return end subroutine psb_ccgr_vect diff --git a/krylov/psb_dcgr.f90 b/krylov/psb_dcgr.f90 index e2377e54..3b60685f 100644 --- a/krylov/psb_dcgr.f90 +++ b/krylov/psb_dcgr.f90 @@ -112,8 +112,8 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,& use psb_d_krylov_conv_mod use psb_krylov_mod implicit none - - + + type(psb_dspmat_type), intent(in) :: a Type(psb_desc_type), Intent(in) :: desc_a class(psb_dprec_type), intent(inout) :: prec @@ -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(out) :: iter real(psb_dpk_), Optional, Intent(out) :: err -! = local data + ! = local data real(psb_dpk_), allocatable :: alpha(:), h(:,:) type(psb_d_vect_type), allocatable :: z(:), c(:), c_scale(:) 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_) :: debug_level, debug_unit integer(psb_ipk_) :: np, me, ictxt @@ -140,14 +140,14 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,& type(psb_itconv_type) :: stopdat character(len=*), parameter :: methdname='CGR' integer(psb_ipk_) ::int_err(5) - info = psb_success_ + info = psb_success_ name = 'psb_dcgr' call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - + ictxt = desc_a%get_context() - + call psb_info(ictxt, me, np) if (.not.allocated(b%v)) then info = psb_err_invalid_vect_state_ @@ -159,21 +159,21 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 endif - + mglob = desc_a%get_global_rows() n_col = desc_a%get_local_cols() - + if (present(istop)) then istop_ = istop else istop_ = 2 endif - -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - + + ! + ! ISTOP_ = 1: Normwise backward error, infinity norm + ! ISTOP_ = 2: ||r||/||b||, 2-norm + ! + if ((istop_ < 1 ).or.(istop_ > 2 ) ) then info=psb_err_invalid_istop_ int_err(1)=istop_ @@ -181,8 +181,8 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,i_err=int_err) goto 9999 endif - - + + call psb_chkvect(mglob,ione,x%get_nrows(),ione,ione,desc_a,info) if (info == psb_success_)& & call psb_chkvect(mglob,ione,b%get_nrows(),ione,ione,desc_a,info) @@ -191,19 +191,19 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,a_err='psb_chkvect on X/B') goto 9999 end if - + if (present(itmax)) then itmax_ = itmax else itmax_ = 1000 endif - + if (present(itrace)) then itrace_ = itrace else itrace_ = 0 end if - + if (present(irst)) then nl = irst if (debug_level >= psb_debug_ext_) & @@ -215,7 +215,7 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,& & write(debug_unit,*) me,' ',trim(name),& & ' not present: irst: ',irst,nl endif - + if (nl <=0 ) then info=psb_err_invalid_istop_ int_err(1)=nl @@ -223,54 +223,44 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,i_err=int_err) goto 9999 endif - + naux=4*n_col allocate(aux(naux),h(nl+1,nl+1),& &c_scale(nl+1),c(nl+1),z(nl+1), alpha(nl+1), stat=info) - + h = dzero if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - - 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 - - itx = 0 - - if (istop_ == 2) then - b_norm = psb_norm2(b, desc_a, 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) + + call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v) + + 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 + + itx = 0 + + nrst = -1 + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) restart: do if (itx>= itmax_) exit restart h = dzero - + it = 0 ! compute r0 = b-ax0 - + if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - - + + call psb_geaxpby(done, b, dzero, r, desc_a, info) call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux) if (info /= psb_success_) then @@ -278,65 +268,65 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 end if - + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) ! if (info /= psb_success_) Then ! call psb_errpush(psb_err_from_subroutine_non_,name) ! goto 9999 ! End If - - nrst = nrst + 1 - + + nrst = nrst + 1 + iteration: do - - itx = itx + 1 - it = it + 1 - j = it - !Apply preconditioner - 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) - do i =1, j - 1 - + + itx = itx + 1 + it = it + 1 + j = it + !Apply preconditioner + 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) + do i =1, j - 1 + 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(-h(i,j), c_scale(i), done, c(i+1), desc_a, info) - end do - - h(j,j) = psb_norm2(c(j), desc_a, info) - hjj = done/h(j,j) - call psb_geaxpby(hjj, c(j), dzero, c_scale(j), desc_a, info) - - alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) - - !Update residual - call psb_geaxpby(done, r, dzero, 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 (j >= irst) exit iteration - - - end do iteration - - m = j - - !Compute solution - - call dtrsm('l','u','n','n',m,1,done,h,size(h,1),alpha,size(alpha,1)) - - if (nrst == 0 ) then - call x%set(dzero) - endif - do i=1,m - call psb_geaxpby(alpha(i), z(i), done, x, desc_a, info) - enddo - - - - + end do + + h(j,j) = psb_norm2(c(j), desc_a, info) + hjj = done/h(j,j) + call psb_geaxpby(hjj, c(j), dzero, c_scale(j), desc_a, info) + + alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) + + !Update residual + call psb_geaxpby(done, r, dzero, 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 (j >= irst) exit iteration + + + end do iteration + + m = j + + !Compute solution + + call dtrsm('l','u','n','n',m,1,done,h,size(h,1),alpha,size(alpha,1)) + + if (nrst == 0 ) then + call x%set(dzero) + endif + do i=1,m + call psb_geaxpby(alpha(i), z(i), done, x, desc_a, info) + enddo + + + + end do restart m = j !Compute solution @@ -347,38 +337,38 @@ subroutine psb_dcgr_vect(a,prec,b,x,eps,desc_a,info,& enddo iter = j - + call psb_end_conv(methdname,itx,desc_a,stopdat,info,derr,iter) if (present(err)) err = derr - + if (info == psb_success_) call psb_gefree(r,desc_a,info) - + do j = 1,m if (info == psb_success_) call psb_gefree(z(j),desc_a,info) if (info == psb_success_) call psb_gefree(c_scale(j),desc_a,info) enddo - + do i =1,nl+1 if (info == psb_success_) call psb_gefree(c(i),desc_a,info) end do - + if (info == psb_success_) deallocate(aux,h,c_scale,z,c,alpha,stat=info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 continue call psb_erractionrestore(err_act) if (err_act.eq.psb_act_abort_) then call psb_error() return end if - + return end subroutine psb_dcgr_vect diff --git a/krylov/psb_scgr.f90 b/krylov/psb_scgr.f90 index 965343b8..5887b332 100644 --- a/krylov/psb_scgr.f90 +++ b/krylov/psb_scgr.f90 @@ -112,8 +112,8 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,& use psb_s_krylov_conv_mod use psb_krylov_mod implicit none - - + + type(psb_sspmat_type), intent(in) :: a Type(psb_desc_type), Intent(in) :: desc_a class(psb_sprec_type), intent(inout) :: prec @@ -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(out) :: iter real(psb_spk_), Optional, Intent(out) :: err -! = local data + ! = local data real(psb_spk_), allocatable :: alpha(:), h(:,:) type(psb_s_vect_type), allocatable :: z(:), c(:), c_scale(:) 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_) :: debug_level, debug_unit integer(psb_ipk_) :: np, me, ictxt @@ -140,14 +140,14 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,& type(psb_itconv_type) :: stopdat character(len=*), parameter :: methdname='CGR' integer(psb_ipk_) ::int_err(5) - info = psb_success_ + info = psb_success_ name = 'psb_scgr' call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - + ictxt = desc_a%get_context() - + call psb_info(ictxt, me, np) if (.not.allocated(b%v)) then info = psb_err_invalid_vect_state_ @@ -159,21 +159,21 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 endif - + mglob = desc_a%get_global_rows() n_col = desc_a%get_local_cols() - + if (present(istop)) then istop_ = istop else istop_ = 2 endif - -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - + + ! + ! ISTOP_ = 1: Normwise backward error, infinity norm + ! ISTOP_ = 2: ||r||/||b||, 2-norm + ! + if ((istop_ < 1 ).or.(istop_ > 2 ) ) then info=psb_err_invalid_istop_ int_err(1)=istop_ @@ -181,8 +181,8 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,i_err=int_err) goto 9999 endif - - + + call psb_chkvect(mglob,ione,x%get_nrows(),ione,ione,desc_a,info) if (info == psb_success_)& & call psb_chkvect(mglob,ione,b%get_nrows(),ione,ione,desc_a,info) @@ -191,19 +191,19 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,a_err='psb_chkvect on X/B') goto 9999 end if - + if (present(itmax)) then itmax_ = itmax else itmax_ = 1000 endif - + if (present(itrace)) then itrace_ = itrace else itrace_ = 0 end if - + if (present(irst)) then nl = irst if (debug_level >= psb_debug_ext_) & @@ -215,7 +215,7 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,& & write(debug_unit,*) me,' ',trim(name),& & ' not present: irst: ',irst,nl endif - + if (nl <=0 ) then info=psb_err_invalid_istop_ int_err(1)=nl @@ -223,54 +223,44 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,i_err=int_err) goto 9999 endif - + naux=4*n_col allocate(aux(naux),h(nl+1,nl+1),& &c_scale(nl+1),c(nl+1),z(nl+1), alpha(nl+1), stat=info) - + h = szero if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - - 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 - - itx = 0 - - if (istop_ == 2) then - b_norm = psb_norm2(b, desc_a, 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) + + call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v) + + 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 + + itx = 0 + + nrst = -1 + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) restart: do if (itx>= itmax_) exit restart h = szero - + it = 0 ! compute r0 = b-ax0 - + if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - - + + call psb_geaxpby(sone, b, szero, r, desc_a, info) call psb_spmm(-sone,a,x,sone,r,desc_a,info,work=aux) if (info /= psb_success_) then @@ -278,65 +268,65 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 end if - + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) ! if (info /= psb_success_) Then ! call psb_errpush(psb_err_from_subroutine_non_,name) ! goto 9999 ! End If - - nrst = nrst + 1 - + + nrst = nrst + 1 + iteration: do - - itx = itx + 1 - it = it + 1 - j = it - !Apply preconditioner - 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) - do i =1, j - 1 - + + itx = itx + 1 + it = it + 1 + j = it + !Apply preconditioner + 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) + do i =1, j - 1 + 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(-h(i,j), c_scale(i), sone, c(i+1), desc_a, info) - end do - - h(j,j) = psb_norm2(c(j), desc_a, info) - hjj = sone/h(j,j) - call psb_geaxpby(hjj, c(j), szero, c_scale(j), desc_a, info) - - alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) - - !Update residual - call psb_geaxpby(sone, r, szero, 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 (j >= irst) exit iteration - - - end do iteration - - m = j - - !Compute solution - - call strsm('l','u','n','n',m,1,sone,h,size(h,1),alpha,size(alpha,1)) - - if (nrst == 0 ) then - call x%set(szero) - endif - do i=1,m - call psb_geaxpby(alpha(i), z(i), sone, x, desc_a, info) - enddo - - - - + end do + + h(j,j) = psb_norm2(c(j), desc_a, info) + hjj = sone/h(j,j) + call psb_geaxpby(hjj, c(j), szero, c_scale(j), desc_a, info) + + alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) + + !Update residual + call psb_geaxpby(sone, r, szero, 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 (j >= irst) exit iteration + + + end do iteration + + m = j + + !Compute solution + + call strsm('l','u','n','n',m,1,sone,h,size(h,1),alpha,size(alpha,1)) + + if (nrst == 0 ) then + call x%set(szero) + endif + do i=1,m + call psb_geaxpby(alpha(i), z(i), sone, x, desc_a, info) + enddo + + + + end do restart m = j !Compute solution @@ -347,38 +337,38 @@ subroutine psb_scgr_vect(a,prec,b,x,eps,desc_a,info,& enddo iter = j - + call psb_end_conv(methdname,itx,desc_a,stopdat,info,derr,iter) if (present(err)) err = derr - + if (info == psb_success_) call psb_gefree(r,desc_a,info) - + do j = 1,m if (info == psb_success_) call psb_gefree(z(j),desc_a,info) if (info == psb_success_) call psb_gefree(c_scale(j),desc_a,info) enddo - + do i =1,nl+1 if (info == psb_success_) call psb_gefree(c(i),desc_a,info) end do - + if (info == psb_success_) deallocate(aux,h,c_scale,z,c,alpha,stat=info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 continue call psb_erractionrestore(err_act) if (err_act.eq.psb_act_abort_) then call psb_error() return end if - + return end subroutine psb_scgr_vect diff --git a/krylov/psb_zcgr.f90 b/krylov/psb_zcgr.f90 index c62d68eb..7d0f60ce 100644 --- a/krylov/psb_zcgr.f90 +++ b/krylov/psb_zcgr.f90 @@ -112,8 +112,8 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,& use psb_z_krylov_conv_mod use psb_krylov_mod implicit none - - + + type(psb_zspmat_type), intent(in) :: a Type(psb_desc_type), Intent(in) :: desc_a class(psb_zprec_type), intent(inout) :: prec @@ -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(out) :: iter real(psb_dpk_), Optional, Intent(out) :: err -! = local data + ! = local data complex(psb_dpk_), allocatable :: alpha(:), h(:,:) type(psb_z_vect_type), allocatable :: z(:), c(:), c_scale(:) 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_) :: debug_level, debug_unit integer(psb_ipk_) :: np, me, ictxt @@ -140,14 +140,14 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,& type(psb_itconv_type) :: stopdat character(len=*), parameter :: methdname='CGR' integer(psb_ipk_) ::int_err(5) - info = psb_success_ + info = psb_success_ name = 'psb_zcgr' call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - + ictxt = desc_a%get_context() - + call psb_info(ictxt, me, np) if (.not.allocated(b%v)) then info = psb_err_invalid_vect_state_ @@ -159,21 +159,21 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 endif - + mglob = desc_a%get_global_rows() n_col = desc_a%get_local_cols() - + if (present(istop)) then istop_ = istop else istop_ = 2 endif - -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - + + ! + ! ISTOP_ = 1: Normwise backward error, infinity norm + ! ISTOP_ = 2: ||r||/||b||, 2-norm + ! + if ((istop_ < 1 ).or.(istop_ > 2 ) ) then info=psb_err_invalid_istop_ int_err(1)=istop_ @@ -181,8 +181,8 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,i_err=int_err) goto 9999 endif - - + + call psb_chkvect(mglob,ione,x%get_nrows(),ione,ione,desc_a,info) if (info == psb_success_)& & call psb_chkvect(mglob,ione,b%get_nrows(),ione,ione,desc_a,info) @@ -191,19 +191,19 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,a_err='psb_chkvect on X/B') goto 9999 end if - + if (present(itmax)) then itmax_ = itmax else itmax_ = 1000 endif - + if (present(itrace)) then itrace_ = itrace else itrace_ = 0 end if - + if (present(irst)) then nl = irst if (debug_level >= psb_debug_ext_) & @@ -215,7 +215,7 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,& & write(debug_unit,*) me,' ',trim(name),& & ' not present: irst: ',irst,nl endif - + if (nl <=0 ) then info=psb_err_invalid_istop_ int_err(1)=nl @@ -223,54 +223,44 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name,i_err=int_err) goto 9999 endif - + naux=4*n_col allocate(aux(naux),h(nl+1,nl+1),& &c_scale(nl+1),c(nl+1),z(nl+1), alpha(nl+1), stat=info) - + h = zzero if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - - 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 - - itx = 0 - - if (istop_ == 2) then - b_norm = psb_norm2(b, desc_a, 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) + + call psb_geasb(r, desc_a,info, scratch=.true.,mold=x%v) + + 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 + + itx = 0 + + nrst = -1 + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) restart: do if (itx>= itmax_) exit restart h = zzero - + it = 0 ! compute r0 = b-ax0 - + if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - - + + call psb_geaxpby(zone, b, zzero, r, desc_a, info) call psb_spmm(-zone,a,x,zone,r,desc_a,info,work=aux) if (info /= psb_success_) then @@ -278,65 +268,65 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,& call psb_errpush(info,name) goto 9999 end if - + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,b,eps,desc_a,stopdat,info) ! if (info /= psb_success_) Then ! call psb_errpush(psb_err_from_subroutine_non_,name) ! goto 9999 ! End If - - nrst = nrst + 1 - + + nrst = nrst + 1 + iteration: do - - itx = itx + 1 - it = it + 1 - j = it - !Apply preconditioner - 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) - do i =1, j - 1 - + + itx = itx + 1 + it = it + 1 + j = it + !Apply preconditioner + 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) + do i =1, j - 1 + 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(-h(i,j), c_scale(i), zone, c(i+1), desc_a, info) - end do - - h(j,j) = psb_norm2(c(j), desc_a, info) - hjj = zone/h(j,j) - call psb_geaxpby(hjj, c(j), zzero, c_scale(j), desc_a, info) - - alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) - - !Update residual - call psb_geaxpby(zone, r, zzero, 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 (j >= irst) exit iteration - - - end do iteration - - m = j - - !Compute solution - - call ztrsm('l','u','n','n',m,1,zone,h,size(h,1),alpha,size(alpha,1)) - - if (nrst == 0 ) then - call x%set(zzero) - endif - do i=1,m - call psb_geaxpby(alpha(i), z(i), zone, x, desc_a, info) - enddo - - - - + end do + + h(j,j) = psb_norm2(c(j), desc_a, info) + hjj = zone/h(j,j) + call psb_geaxpby(hjj, c(j), zzero, c_scale(j), desc_a, info) + + alpha(j) = psb_gedot(c_scale(j), r, desc_a, info) + + !Update residual + call psb_geaxpby(zone, r, zzero, 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 (j >= irst) exit iteration + + + end do iteration + + m = j + + !Compute solution + + call ztrsm('l','u','n','n',m,1,zone,h,size(h,1),alpha,size(alpha,1)) + + if (nrst == 0 ) then + call x%set(zzero) + endif + do i=1,m + call psb_geaxpby(alpha(i), z(i), zone, x, desc_a, info) + enddo + + + + end do restart m = j !Compute solution @@ -347,38 +337,38 @@ subroutine psb_zcgr_vect(a,prec,b,x,eps,desc_a,info,& enddo iter = j - + call psb_end_conv(methdname,itx,desc_a,stopdat,info,derr,iter) if (present(err)) err = derr - + if (info == psb_success_) call psb_gefree(r,desc_a,info) - + do j = 1,m if (info == psb_success_) call psb_gefree(z(j),desc_a,info) if (info == psb_success_) call psb_gefree(c_scale(j),desc_a,info) enddo - + do i =1,nl+1 if (info == psb_success_) call psb_gefree(c(i),desc_a,info) end do - + if (info == psb_success_) deallocate(aux,h,c_scale,z,c,alpha,stat=info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 continue call psb_erractionrestore(err_act) if (err_act.eq.psb_act_abort_) then call psb_error() return end if - + return end subroutine psb_zcgr_vect