diff --git a/krylov/Makefile b/krylov/Makefile index f71dbb1a..c08d6740 100644 --- a/krylov/Makefile +++ b/krylov/Makefile @@ -11,6 +11,7 @@ MODOBJS= psb_base_krylov_conv_mod.o \ psb_d_krylov_conv_mod.o psb_z_krylov_conv_mod.o \ psb_krylov_mod.o F90OBJS=psb_dkrylov.o psb_skrylov.o psb_ckrylov.o psb_zkrylov.o \ + psb_drichardson.o psb_srichardson.o psb_crichardson.o psb_zrichardson.o \ psb_dcgstab.o psb_dcg.o psb_dfcg.o psb_dgcr.o psb_dcgs.o \ psb_dbicg.o psb_dcgstabl.o psb_drgmres.o\ psb_scgstab.o psb_scg.o psb_sfcg.o psb_sgcr.o psb_scgs.o \ diff --git a/krylov/psb_ckrylov.f90 b/krylov/psb_ckrylov.f90 index 01228234..3afa525e 100644 --- a/krylov/psb_ckrylov.f90 +++ b/krylov/psb_ckrylov.f90 @@ -42,6 +42,7 @@ ! ! methd - character The specific method; can take the values: ! CG +! FCG ! CGS ! BICG ! BICGSTAB diff --git a/krylov/psb_crichardson.f90 b/krylov/psb_crichardson.f90 new file mode 100644 index 00000000..08678653 --- /dev/null +++ b/krylov/psb_crichardson.f90 @@ -0,0 +1,216 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! 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. +! +! +! +! File: psb_richardson_mod.f90 +! Interfaces for Richardson iterative methods. +! +! +! Subroutine: psb_crichardson +! +! Front-end for the Richardson iterations, complexversion +! +! Arguments: +! +! a - type(psb_cspmat_type) Input: sparse matrix containing A. +! prec - class(psb_cprec_type) Input: preconditioner +! b - complex,dimension(:) Input: vector containing the +! right hand side B +! x - complex,dimension(:) Input/Output: vector containing the +! initial guess and final solution X. +! eps - real Input: Stopping tolerance; the iteration is +! stopped when the error +! estimate |err| <= eps +! +! desc_a - type(psb_desc_type). Input: The communication descriptor. +! info - integer. Output: Return code +! +! itmax - integer(optional) Input: maximum number of iterations to be +! performed. +! iter - integer(optional) Output: how many iterations have been +! performed. +! err - real (optional) Output: error estimate on exit +! itrace - integer(optional) Input: print an informational message +! with the error estimate every itrace +! iterations +! istop - integer(optional) Input: stopping criterion, or how +! to estimate the error. +! 1: err = |r|/(|a||x|+|b|) +! 2: err = |r|/|b| +! where r is the (preconditioned, recursive +! estimate of) residual +! +Subroutine psb_crichardson_vect(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,istop) + + use psb_base_mod + use psb_prec_mod + use psb_c_krylov_conv_mod + use psb_krylov_mod, psb_protect_name => psb_crichardson_vect + + Type(psb_cspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_cprec_type), intent(inout) :: prec + type(psb_c_vect_type), Intent(inout) :: b + type(psb_c_vect_type), Intent(inout) :: x + Real(psb_spk_), Intent(in) :: eps + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop + integer(psb_ipk_), Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err + + + logical :: do_alloc_wrk + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: me,np,err_act + complex(psb_spk_), allocatable, target :: aux(:) + type(psb_c_vect_type), allocatable, target :: wwrk(:) + type(psb_c_vect_type), pointer :: q, p, r, z, w + real(psb_dpk_) :: derr + integer(psb_ipk_) :: itmax_, istop_, naux, it, itx, itrace_,& + & n_col, n_row,ieg,nspl, istebz + integer(psb_lpk_) :: mglob + integer(psb_ipk_) :: debug_level, debug_unit + type(psb_itconv_type) :: stopdat + character(len=20) :: name + character(len=*), parameter :: methdname='RICHARDSON' + + info = psb_success_ + name = 'psb_crichardson' + call psb_erractionsave(err_act) + + ctxt=desc_a%get_context() + + call psb_info(ctxt, me, np) + + if (present(itrace)) then + itrace_ = itrace + else + itrace_ = -1 + end if + + if (present(istop)) then + istop_ = istop + else + istop_ = 2 + endif + if (present(itmax)) then + itmax_ = itmax + else + itmax_ = 1000 + endif + + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a) + + if (.not.allocated(b%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + mglob = desc_a%get_global_rows() + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) + if (info == psb_success_)& + & call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_chkvect on X/B') + goto 9999 + end if + + naux=4*n_col + allocate(aux(naux), stat=info) + if (info == psb_success_) call psb_geall(wwrk,desc_a,info,n=5_psb_ipk_) + if (info == psb_success_) call psb_geasb(wwrk,desc_a,info,mold=x%v,scratch=.true.) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + p => wwrk(1) + q => wwrk(2) + r => wwrk(3) + z => wwrk(4) + w => wwrk(5) + + 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,work=aux) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + if (info /= psb_success_) Then + call psb_errpush(psb_err_from_subroutine_non_,name) + goto 9999 + End If + + loop: do itx=1,itmax_ + call prec%apply(r,z,desc_a,info,work=aux) + call psb_geaxpby(cone,z,cone,x,desc_a,info) + 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,work=aux) + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit loop + end do loop + call psb_end_conv(methdname,itx,desc_a,stopdat,info,derr,iter) + if (present(err)) err = derr + + if (info == psb_success_) call psb_gefree(wwrk,desc_a,info) + if (info == psb_success_) deallocate(aux,stat=info) + if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) + + if(info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info,name,a_err=trim(methdname)) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + +end subroutine psb_crichardson_vect + diff --git a/krylov/psb_dkrylov.f90 b/krylov/psb_dkrylov.f90 index d5d40eaf..2bc24d6a 100644 --- a/krylov/psb_dkrylov.f90 +++ b/krylov/psb_dkrylov.f90 @@ -42,6 +42,7 @@ ! ! methd - character The specific method; can take the values: ! CG +! FCG ! CGS ! BICG ! BICGSTAB diff --git a/krylov/psb_drichardson.f90 b/krylov/psb_drichardson.f90 new file mode 100644 index 00000000..2c057320 --- /dev/null +++ b/krylov/psb_drichardson.f90 @@ -0,0 +1,216 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! 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. +! +! +! +! File: psb_richardson_mod.f90 +! Interfaces for Richardson iterative methods. +! +! +! Subroutine: psb_drichardson +! +! Front-end for the Richardson iterations, realversion +! +! Arguments: +! +! a - type(psb_dspmat_type) Input: sparse matrix containing A. +! prec - class(psb_dprec_type) Input: preconditioner +! b - real,dimension(:) Input: vector containing the +! right hand side B +! x - real,dimension(:) Input/Output: vector containing the +! initial guess and final solution X. +! eps - real Input: Stopping tolerance; the iteration is +! stopped when the error +! estimate |err| <= eps +! +! desc_a - type(psb_desc_type). Input: The communication descriptor. +! info - integer. Output: Return code +! +! itmax - integer(optional) Input: maximum number of iterations to be +! performed. +! iter - integer(optional) Output: how many iterations have been +! performed. +! err - real (optional) Output: error estimate on exit +! itrace - integer(optional) Input: print an informational message +! with the error estimate every itrace +! iterations +! istop - integer(optional) Input: stopping criterion, or how +! to estimate the error. +! 1: err = |r|/(|a||x|+|b|) +! 2: err = |r|/|b| +! where r is the (preconditioned, recursive +! estimate of) residual +! +Subroutine psb_drichardson_vect(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,istop) + + use psb_base_mod + use psb_prec_mod + use psb_d_krylov_conv_mod + use psb_krylov_mod, psb_protect_name => psb_drichardson_vect + + Type(psb_dspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_dprec_type), intent(inout) :: prec + type(psb_d_vect_type), Intent(inout) :: b + type(psb_d_vect_type), Intent(inout) :: x + Real(psb_dpk_), Intent(in) :: eps + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop + integer(psb_ipk_), Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err + + + logical :: do_alloc_wrk + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: me,np,err_act + real(psb_dpk_), allocatable, target :: aux(:) + type(psb_d_vect_type), allocatable, target :: wwrk(:) + type(psb_d_vect_type), pointer :: q, p, r, z, w + real(psb_dpk_) :: derr + integer(psb_ipk_) :: itmax_, istop_, naux, it, itx, itrace_,& + & n_col, n_row,ieg,nspl, istebz + integer(psb_lpk_) :: mglob + integer(psb_ipk_) :: debug_level, debug_unit + type(psb_itconv_type) :: stopdat + character(len=20) :: name + character(len=*), parameter :: methdname='RICHARDSON' + + info = psb_success_ + name = 'psb_drichardson' + call psb_erractionsave(err_act) + + ctxt=desc_a%get_context() + + call psb_info(ctxt, me, np) + + if (present(itrace)) then + itrace_ = itrace + else + itrace_ = -1 + end if + + if (present(istop)) then + istop_ = istop + else + istop_ = 2 + endif + if (present(itmax)) then + itmax_ = itmax + else + itmax_ = 1000 + endif + + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a) + + if (.not.allocated(b%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + mglob = desc_a%get_global_rows() + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) + if (info == psb_success_)& + & call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_chkvect on X/B') + goto 9999 + end if + + naux=4*n_col + allocate(aux(naux), stat=info) + if (info == psb_success_) call psb_geall(wwrk,desc_a,info,n=5_psb_ipk_) + if (info == psb_success_) call psb_geasb(wwrk,desc_a,info,mold=x%v,scratch=.true.) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + p => wwrk(1) + q => wwrk(2) + r => wwrk(3) + z => wwrk(4) + w => wwrk(5) + + 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,work=aux) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + if (info /= psb_success_) Then + call psb_errpush(psb_err_from_subroutine_non_,name) + goto 9999 + End If + + loop: do itx=1,itmax_ + call prec%apply(r,z,desc_a,info,work=aux) + call psb_geaxpby(done,z,done,x,desc_a,info) + 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,work=aux) + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit loop + end do loop + call psb_end_conv(methdname,itx,desc_a,stopdat,info,derr,iter) + if (present(err)) err = derr + + if (info == psb_success_) call psb_gefree(wwrk,desc_a,info) + if (info == psb_success_) deallocate(aux,stat=info) + if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) + + if(info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info,name,a_err=trim(methdname)) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + +end subroutine psb_drichardson_vect + diff --git a/krylov/psb_krylov_mod.f90 b/krylov/psb_krylov_mod.f90 index d8d4d904..e9a94e18 100644 --- a/krylov/psb_krylov_mod.f90 +++ b/krylov/psb_krylov_mod.f90 @@ -127,4 +127,88 @@ Module psb_krylov_mod end interface + interface psb_richardson + + Subroutine psb_srichardson_vect(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,istop) + + use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_sspmat_type, & + & psb_spk_, psb_s_vect_type + use psb_prec_mod, only : psb_sprec_type + + Type(psb_sspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_sprec_type), intent(inout) :: prec + type(psb_s_vect_type), Intent(inout) :: b + type(psb_s_vect_type), Intent(inout) :: x + Real(psb_spk_), Intent(in) :: eps + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop + integer(psb_ipk_), Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err + + end Subroutine psb_srichardson_vect + + Subroutine psb_crichardson_vect(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,istop) + + use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_cspmat_type, & + & psb_spk_, psb_c_vect_type + use psb_prec_mod, only : psb_cprec_type + + Type(psb_cspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_cprec_type), intent(inout) :: prec + type(psb_c_vect_type), Intent(inout) :: b + type(psb_c_vect_type), Intent(inout) :: x + Real(psb_spk_), Intent(in) :: eps + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop + integer(psb_ipk_), Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err + + end Subroutine psb_crichardson_vect + + Subroutine psb_drichardson_vect(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,istop) + + use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_dspmat_type, & + & psb_dpk_, psb_d_vect_type + use psb_prec_mod, only : psb_dprec_type + + Type(psb_dspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_dprec_type), intent(inout) :: prec + type(psb_d_vect_type), Intent(inout) :: b + type(psb_d_vect_type), Intent(inout) :: x + Real(psb_dpk_), Intent(in) :: eps + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop + integer(psb_ipk_), Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err + + end Subroutine psb_drichardson_vect + + Subroutine psb_zrichardson_vect(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,istop) + + use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_zspmat_type, & + & psb_dpk_, psb_z_vect_type + use psb_prec_mod, only : psb_zprec_type + + Type(psb_zspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_zprec_type), intent(inout) :: prec + type(psb_z_vect_type), Intent(inout) :: b + type(psb_z_vect_type), Intent(inout) :: x + Real(psb_dpk_), Intent(in) :: eps + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop + integer(psb_ipk_), Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err + + end Subroutine psb_zrichardson_vect + + end interface + end module psb_krylov_mod diff --git a/krylov/psb_skrylov.f90 b/krylov/psb_skrylov.f90 index 39aecc36..35d2024f 100644 --- a/krylov/psb_skrylov.f90 +++ b/krylov/psb_skrylov.f90 @@ -42,6 +42,7 @@ ! ! methd - character The specific method; can take the values: ! CG +! FCG ! CGS ! BICG ! BICGSTAB diff --git a/krylov/psb_srichardson.f90 b/krylov/psb_srichardson.f90 new file mode 100644 index 00000000..a06f6cb4 --- /dev/null +++ b/krylov/psb_srichardson.f90 @@ -0,0 +1,216 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! 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. +! +! +! +! File: psb_richardson_mod.f90 +! Interfaces for Richardson iterative methods. +! +! +! Subroutine: psb_srichardson +! +! Front-end for the Richardson iterations, realversion +! +! Arguments: +! +! a - type(psb_sspmat_type) Input: sparse matrix containing A. +! prec - class(psb_sprec_type) Input: preconditioner +! b - real,dimension(:) Input: vector containing the +! right hand side B +! x - real,dimension(:) Input/Output: vector containing the +! initial guess and final solution X. +! eps - real Input: Stopping tolerance; the iteration is +! stopped when the error +! estimate |err| <= eps +! +! desc_a - type(psb_desc_type). Input: The communication descriptor. +! info - integer. Output: Return code +! +! itmax - integer(optional) Input: maximum number of iterations to be +! performed. +! iter - integer(optional) Output: how many iterations have been +! performed. +! err - real (optional) Output: error estimate on exit +! itrace - integer(optional) Input: print an informational message +! with the error estimate every itrace +! iterations +! istop - integer(optional) Input: stopping criterion, or how +! to estimate the error. +! 1: err = |r|/(|a||x|+|b|) +! 2: err = |r|/|b| +! where r is the (preconditioned, recursive +! estimate of) residual +! +Subroutine psb_srichardson_vect(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,istop) + + use psb_base_mod + use psb_prec_mod + use psb_s_krylov_conv_mod + use psb_krylov_mod, psb_protect_name => psb_srichardson_vect + + Type(psb_sspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_sprec_type), intent(inout) :: prec + type(psb_s_vect_type), Intent(inout) :: b + type(psb_s_vect_type), Intent(inout) :: x + Real(psb_spk_), Intent(in) :: eps + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop + integer(psb_ipk_), Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err + + + logical :: do_alloc_wrk + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: me,np,err_act + real(psb_spk_), allocatable, target :: aux(:) + type(psb_s_vect_type), allocatable, target :: wwrk(:) + type(psb_s_vect_type), pointer :: q, p, r, z, w + real(psb_dpk_) :: derr + integer(psb_ipk_) :: itmax_, istop_, naux, it, itx, itrace_,& + & n_col, n_row,ieg,nspl, istebz + integer(psb_lpk_) :: mglob + integer(psb_ipk_) :: debug_level, debug_unit + type(psb_itconv_type) :: stopdat + character(len=20) :: name + character(len=*), parameter :: methdname='RICHARDSON' + + info = psb_success_ + name = 'psb_srichardson' + call psb_erractionsave(err_act) + + ctxt=desc_a%get_context() + + call psb_info(ctxt, me, np) + + if (present(itrace)) then + itrace_ = itrace + else + itrace_ = -1 + end if + + if (present(istop)) then + istop_ = istop + else + istop_ = 2 + endif + if (present(itmax)) then + itmax_ = itmax + else + itmax_ = 1000 + endif + + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a) + + if (.not.allocated(b%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + mglob = desc_a%get_global_rows() + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) + if (info == psb_success_)& + & call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_chkvect on X/B') + goto 9999 + end if + + naux=4*n_col + allocate(aux(naux), stat=info) + if (info == psb_success_) call psb_geall(wwrk,desc_a,info,n=5_psb_ipk_) + if (info == psb_success_) call psb_geasb(wwrk,desc_a,info,mold=x%v,scratch=.true.) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + p => wwrk(1) + q => wwrk(2) + r => wwrk(3) + z => wwrk(4) + w => wwrk(5) + + 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,work=aux) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + if (info /= psb_success_) Then + call psb_errpush(psb_err_from_subroutine_non_,name) + goto 9999 + End If + + loop: do itx=1,itmax_ + call prec%apply(r,z,desc_a,info,work=aux) + call psb_geaxpby(sone,z,sone,x,desc_a,info) + 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,work=aux) + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit loop + end do loop + call psb_end_conv(methdname,itx,desc_a,stopdat,info,derr,iter) + if (present(err)) err = derr + + if (info == psb_success_) call psb_gefree(wwrk,desc_a,info) + if (info == psb_success_) deallocate(aux,stat=info) + if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) + + if(info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info,name,a_err=trim(methdname)) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + +end subroutine psb_srichardson_vect + diff --git a/krylov/psb_zkrylov.f90 b/krylov/psb_zkrylov.f90 index a70cc98a..bcfd6806 100644 --- a/krylov/psb_zkrylov.f90 +++ b/krylov/psb_zkrylov.f90 @@ -42,6 +42,7 @@ ! ! methd - character The specific method; can take the values: ! CG +! FCG ! CGS ! BICG ! BICGSTAB diff --git a/krylov/psb_zrichardson.f90 b/krylov/psb_zrichardson.f90 new file mode 100644 index 00000000..19f2f10e --- /dev/null +++ b/krylov/psb_zrichardson.f90 @@ -0,0 +1,216 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! 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. +! +! +! +! File: psb_richardson_mod.f90 +! Interfaces for Richardson iterative methods. +! +! +! Subroutine: psb_zrichardson +! +! Front-end for the Richardson iterations, complexversion +! +! Arguments: +! +! a - type(psb_zspmat_type) Input: sparse matrix containing A. +! prec - class(psb_zprec_type) Input: preconditioner +! b - complex,dimension(:) Input: vector containing the +! right hand side B +! x - complex,dimension(:) Input/Output: vector containing the +! initial guess and final solution X. +! eps - real Input: Stopping tolerance; the iteration is +! stopped when the error +! estimate |err| <= eps +! +! desc_a - type(psb_desc_type). Input: The communication descriptor. +! info - integer. Output: Return code +! +! itmax - integer(optional) Input: maximum number of iterations to be +! performed. +! iter - integer(optional) Output: how many iterations have been +! performed. +! err - real (optional) Output: error estimate on exit +! itrace - integer(optional) Input: print an informational message +! with the error estimate every itrace +! iterations +! istop - integer(optional) Input: stopping criterion, or how +! to estimate the error. +! 1: err = |r|/(|a||x|+|b|) +! 2: err = |r|/|b| +! where r is the (preconditioned, recursive +! estimate of) residual +! +Subroutine psb_zrichardson_vect(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,istop) + + use psb_base_mod + use psb_prec_mod + use psb_z_krylov_conv_mod + use psb_krylov_mod, psb_protect_name => psb_zrichardson_vect + + Type(psb_zspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_zprec_type), intent(inout) :: prec + type(psb_z_vect_type), Intent(inout) :: b + type(psb_z_vect_type), Intent(inout) :: x + Real(psb_dpk_), Intent(in) :: eps + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop + integer(psb_ipk_), Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err + + + logical :: do_alloc_wrk + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: me,np,err_act + complex(psb_dpk_), allocatable, target :: aux(:) + type(psb_z_vect_type), allocatable, target :: wwrk(:) + type(psb_z_vect_type), pointer :: q, p, r, z, w + real(psb_dpk_) :: derr + integer(psb_ipk_) :: itmax_, istop_, naux, it, itx, itrace_,& + & n_col, n_row,ieg,nspl, istebz + integer(psb_lpk_) :: mglob + integer(psb_ipk_) :: debug_level, debug_unit + type(psb_itconv_type) :: stopdat + character(len=20) :: name + character(len=*), parameter :: methdname='RICHARDSON' + + info = psb_success_ + name = 'psb_zrichardson' + call psb_erractionsave(err_act) + + ctxt=desc_a%get_context() + + call psb_info(ctxt, me, np) + + if (present(itrace)) then + itrace_ = itrace + else + itrace_ = -1 + end if + + if (present(istop)) then + istop_ = istop + else + istop_ = 2 + endif + if (present(itmax)) then + itmax_ = itmax + else + itmax_ = 1000 + endif + + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a) + + if (.not.allocated(b%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + mglob = desc_a%get_global_rows() + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) + if (info == psb_success_)& + & call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_chkvect on X/B') + goto 9999 + end if + + naux=4*n_col + allocate(aux(naux), stat=info) + if (info == psb_success_) call psb_geall(wwrk,desc_a,info,n=5_psb_ipk_) + if (info == psb_success_) call psb_geasb(wwrk,desc_a,info,mold=x%v,scratch=.true.) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + p => wwrk(1) + q => wwrk(2) + r => wwrk(3) + z => wwrk(4) + w => wwrk(5) + + 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,work=aux) + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + + call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info) + if (info /= psb_success_) Then + call psb_errpush(psb_err_from_subroutine_non_,name) + goto 9999 + End If + + loop: do itx=1,itmax_ + call prec%apply(r,z,desc_a,info,work=aux) + call psb_geaxpby(zone,z,zone,x,desc_a,info) + 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,work=aux) + if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit loop + end do loop + call psb_end_conv(methdname,itx,desc_a,stopdat,info,derr,iter) + if (present(err)) err = derr + + if (info == psb_success_) call psb_gefree(wwrk,desc_a,info) + if (info == psb_success_) deallocate(aux,stat=info) + if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) + + if(info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info,name,a_err=trim(methdname)) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + +end subroutine psb_zrichardson_vect + diff --git a/test/cudakern/dpdegenmv.F90 b/test/cudakern/dpdegenmv.F90 index 85059e81..a5463f0b 100644 --- a/test/cudakern/dpdegenmv.F90 +++ b/test/cudakern/dpdegenmv.F90 @@ -594,7 +594,7 @@ program pdgenmv ! solver parameters integer(psb_epk_) :: amatsize, precsize, descsize, annz, nbytes - real(psb_dpk_) :: err, eps + real(psb_dpk_) :: err, eps, tnv, tng,tdot, dnrm2,ddot integer, parameter :: ntests=200, ngpu=50, ncnv=20 type(psb_d_coo_sparse_mat), target :: acoo type(psb_d_csr_sparse_mat), target :: acsr @@ -745,7 +745,7 @@ program pdgenmv call psb_geall(x0,desc_a,info) do i=1, nr call desc_a%l2g(i,ig,info) - x0(i) = 1.0 + (1.0*ig)/nrg + x0(i) = 1.0 + (1.0*ig)/(nrg**2) end do call a%cscnv(aux_a,info,mold=acoo) tcnvcsr = 0 @@ -843,6 +843,12 @@ program pdgenmv call bg%sync() x1 = bv%get_vect() x2 = bg%get_vect() + tnv = psb_genrm2(bv,desc_a,info) + tng = psb_genrm2(bg,desc_a,info) + tdot = psb_gedot(bg,bg,desc_a,info) + write(0,*) ' bv ',tnv,' bg ',tng, ' dot ',tdot,eps,& + & dnrm2(desc_a%get_local_rows(),x2,1),& + & ddot(desc_a%get_local_rows(),x1,1,x2,1) call psb_geaxpby(-done,bg,+done,bv,desc_a,info) eps = psb_geamax(bv,desc_a,info) diff --git a/test/cudakern/spdegenmv.F90 b/test/cudakern/spdegenmv.F90 index f953e163..fbce6726 100644 --- a/test/cudakern/spdegenmv.F90 +++ b/test/cudakern/spdegenmv.F90 @@ -580,7 +580,7 @@ program pdgenmv ! solver parameters integer(psb_epk_) :: amatsize, precsize, descsize, annz, nbytes - real(psb_spk_) :: err, eps + real(psb_spk_) :: err, eps, tnv, tng,tdot, snrm2,sdot integer, parameter :: ntests=200, ngpu=50, ncnv=20 type(psb_s_coo_sparse_mat), target :: acoo type(psb_s_csr_sparse_mat), target :: acsr @@ -728,7 +728,7 @@ program pdgenmv call psb_geall(x0,desc_a,info) do i=1, nr call desc_a%l2g(i,ig,info) - x0(i) = 1.0 + (1.0*ig)/nrg + x0(i) = 1.0 + (1.0*ig)/(nrg**2) end do call a%cscnv(aux_a,info,mold=acoo) tcnvcsr = 0 @@ -826,10 +826,16 @@ program pdgenmv call psb_amx(ctxt,gt2) call bg%sync() x1 = bv%get_vect() - x2 = bg%get_vect() + x2 = bg%get_vect() + tnv = psb_genrm2(bv,desc_a,info) + tng = psb_genrm2(bg,desc_a,info) + tdot = psb_gedot(bg,bg,desc_a,info) + write(0,*) ' bv ',tnv,' bg ',tng, ' dot ',tdot,eps,& + & snrm2(desc_a%get_local_rows(),x2,1),& + & sdot(desc_a%get_local_rows(),x1,1,x2,1) call psb_geaxpby(-sone,bg,+sone,bv,desc_a,info) eps = psb_geamax(bv,desc_a,info) - + call psb_amx(ctxt,t2) eps = maxval(abs(x1(1:nr)-x2(1:nr))) call psb_amx(ctxt,eps)