diff --git a/base/modules/Makefile b/base/modules/Makefile index fedc33e0..8ce37c35 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -5,7 +5,7 @@ UTIL_MODS = psb_string_mod.o \ psb_desc_type.o psb_sort_mod.o psb_serial_mod.o \ psb_base_tools_mod.o psb_s_tools_mod.o psb_d_tools_mod.o\ psb_c_tools_mod.o psb_z_tools_mod.o psb_tools_mod.o \ - psi_comm_buffers_mod.o psi_penv_mod.o psi_bcast_mod.o \ + psb_penv_mod.o psi_comm_buffers_mod.o psi_penv_mod.o psi_bcast_mod.o \ psi_reduce_mod.o psi_p2p_mod.o psb_error_impl.o \ psb_linmap_type_mod.o psb_linmap_mod.o psb_comm_mod.o\ psb_s_psblas_mod.o psb_c_psblas_mod.o \ diff --git a/configure b/configure index cc555ed8..bad37952 100755 --- a/configure +++ b/configure @@ -7063,7 +7063,7 @@ if test "X$psblas_cv_fc" == X"nag" ; then F90COPT="$F90COPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv" EXTRA_OPT="-mismatch_all" F03COPT="${F90COPT}" - F03="nagfor" + F03="${MPIFC}" elif test "X$psblas_cv_fc" == X"xlf" ; then F03="xlf2003_r" F03COPT="-O3 -qarch=auto -qsuffix=f=f03:cpp=F03 $F03COPT" diff --git a/configure.ac b/configure.ac index fca9a9ea..7370ba53 100755 --- a/configure.ac +++ b/configure.ac @@ -428,7 +428,7 @@ if test "X$psblas_cv_fc" == X"nag" ; then F90COPT="$F90COPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv" EXTRA_OPT="-mismatch_all" F03COPT="${F90COPT}" - F03="nagfor" + F03="${MPIFC}" elif test "X$psblas_cv_fc" == X"xlf" ; then F03="xlf2003_r" F03COPT="-O3 -qarch=auto -qsuffix=f=f03:cpp=F03 $F03COPT" diff --git a/krylov/Makefile b/krylov/Makefile index f180e1c0..1873a0f4 100644 --- a/krylov/Makefile +++ b/krylov/Makefile @@ -4,8 +4,11 @@ include ../Make.inc HERE=. LIBDIR=../lib -MODOBJS= psb_krylov_mod.o -F90OBJS= psb_dcgstab.o psb_dcg.o psb_dcgs.o \ +MODOBJS= psb_base_inner_krylov_mod.o \ + psb_s_inner_krylov_mod.o psb_c_inner_krylov_mod.o psb_d_inner_krylov_mod.o psb_z_inner_krylov_mod.o \ + psb_inner_krylov_mod.o psb_krylov_mod.o +F90OBJS=psb_dkrylov.o psb_skrylov.o psb_ckrylov.o psb_zkrylov.o \ + psb_dcgstab.o psb_dcg.o psb_dcgs.o \ psb_dbicg.o psb_dcgstabl.o psb_drgmres.o\ psb_scgstab.o psb_scg.o psb_scgs.o \ psb_sbicg.o psb_scgstabl.o psb_srgmres.o\ @@ -16,7 +19,7 @@ F90OBJS= psb_dcgstab.o psb_dcg.o psb_dcgs.o \ OBJS=$(F90OBJS) $(MODOBJS) LIBMOD=psb_krylov_mod$(.mod) -LOCAL_MODS=$(LIBMOD) +LOCAL_MODS=$(MODOBJS:.o=$(.mod)) LIBNAME=$(METHDLIBNAME) FINCLUDES=$(FMFLAG)$(LIBDIR) $(FMFLAG). @@ -27,6 +30,8 @@ lib: $(OBJS) /bin/cp -p $(HERE)/$(LIBNAME) $(LIBDIR) /bin/cp -p $(LIBMOD) $(LIBDIR) +psb_s_inner_krylov_mod.o psb_c_inner_krylov_mod.o psb_d_inner_krylov_mod.o psb_z_inner_krylov_mod.o: psb_base_inner_krylov_mod.o +psb_inner_krylov_mod.o: psb_s_inner_krylov_mod.o psb_c_inner_krylov_mod.o psb_d_inner_krylov_mod.o psb_z_inner_krylov_mod.o $(F90OBJS): $(MODOBJS) $(OBJS): $(LIBDIR)/psb_prec_mod$(.mod) $(LIBDIR)/psb_sparse_mod$(.mod) veryclean: clean diff --git a/krylov/psb_base_inner_krylov_mod.f90 b/krylov/psb_base_inner_krylov_mod.f90 new file mode 100644 index 00000000..225016f2 --- /dev/null +++ b/krylov/psb_base_inner_krylov_mod.f90 @@ -0,0 +1,161 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.2 +!!$ (C) Copyright 2006/2007/2008 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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_krylov_mod.f90 +! Interfaces for Krylov subspace iterative methods. +! +Module psb_base_inner_krylov_mod + + use psb_const_mod + + interface psb_end_conv + module procedure psb_d_end_conv + end interface + + integer, parameter :: psb_ik_bni_=1, psb_ik_rni_=2, psb_ik_ani_=3 + integer, parameter :: psb_ik_xni_=4, psb_ik_bn2_=5, psb_ik_xn2_=6 + integer, parameter :: psb_ik_errnum_=7, psb_ik_errden_=8, psb_ik_eps_=9, psb_ik_rn2_=10 + integer, parameter :: psb_ik_stopc_=1, psb_ik_trace_=2, psb_ik_itmax_=3 + integer, parameter :: psb_ik_ivsz_=16 + type psb_itconv_type + integer :: controls(psb_ik_ivsz_) + real(psb_dpk_) :: values(psb_ik_ivsz_) + end type psb_itconv_type + +contains + + subroutine log_header(methdname) + !use psb_base_mod + implicit none + character(len=*), intent(in) :: methdname + character(len=*), parameter :: fmt='(a18,1x,a4,3(2x,a15))' + integer, parameter :: outlen=18 + character(len=len(methdname)) :: mname + character(len=outlen) :: outname + + mname = adjustl(trim(methdname)) + write(outname,'(a)') mname(1:min(len_trim(mname),outlen-1))//':' + write(*,fmt) adjustl(outname),'Iteration','Error Estimate','Tolerance' + + end subroutine log_header + + + subroutine log_conv(methdname,me,itx,itrace,errnum,errden,eps) + !use psb_base_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: me, itx, itrace + real(psb_dpk_), intent(in) :: errnum, errden, eps + character(len=*), parameter :: fmt='(a18,1x,i4,3(2x,es15.9))' + integer, parameter :: outlen=18 + character(len=len(methdname)) :: mname + character(len=outlen) :: outname + + if ((mod(itx,itrace) == 0).and.(me == 0)) then + mname = adjustl(trim(methdname)) + write(outname,'(a)') mname(1:min(len_trim(mname),outlen-1))//':' + if (errden > dzero ) then + write(*,fmt) adjustl(outname),itx,errnum/errden,eps + else + write(*,fmt) adjustl(outname),itx,errnum,eps + end if + endif + + end subroutine log_conv + + subroutine log_end(methdname,me,it,errnum,errden,eps,err,iter) + !use psb_base_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: me, it + real(psb_dpk_), intent(in) :: errnum, errden, eps + real(psb_dpk_), optional, intent(out) :: err + integer, optional, intent(out) :: iter + + character(len=*), parameter :: fmt='(a,2x,es15.9,1x,a,1x,i4,1x,a)' + character(len=*), parameter :: fmt1='(a,3(2x,es15.9))' + + if (errden == dzero) then + if (errnum > eps) then + if (me == 0) then + write(*,fmt) trim(methdname)//' failed to converge to ',eps,& + & ' in ',it,' iterations. ' + write(*,fmt1) 'Last iteration error estimate: ',& + & errnum + end if + end if + if (present(err)) err=errnum + else + if (errnum/errden > eps) then + if (me == 0) then + write(*,fmt) trim(methdname)//' failed to converge to ',eps,& + & ' in ',it,' iterations. ' + write(*,fmt1) 'Last iteration error estimate: ',& + & errnum/errden + end if + endif + if (present(err)) err=errnum/errden + end if + if (present(iter)) iter = it + + end subroutine log_end + + + subroutine psb_d_end_conv(methdname,it,desc_a,stopdat,info,err,iter) + use psb_sparse_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: it + type(psb_desc_type), intent(in) :: desc_a + type(psb_itconv_type) :: stopdat + integer, intent(out) :: info + real(psb_dpk_), optional, intent(out) :: err + integer, optional, intent(out) :: iter + + integer :: ictxt, me, np, err_act + real(psb_dpk_) :: errnum, errden, eps + character(len=20) :: name + + info = psb_success_ + name = 'psb_end_conv' + + ictxt = psb_cd_get_context(desc_a) + call psb_info(ictxt,me,np) + + errnum = stopdat%values(psb_ik_errnum_) + errden = stopdat%values(psb_ik_errden_) + eps = stopdat%values(psb_ik_eps_) + call log_end(methdname,me,it,errnum,errden,eps,err,iter) + + end subroutine psb_d_end_conv + +end module psb_base_inner_krylov_mod diff --git a/krylov/psb_c_inner_krylov_mod.f90 b/krylov/psb_c_inner_krylov_mod.f90 new file mode 100644 index 00000000..cc29742d --- /dev/null +++ b/krylov/psb_c_inner_krylov_mod.f90 @@ -0,0 +1,191 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.2 +!!$ (C) Copyright 2006/2007/2008 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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_krylov_mod.f90 +! Interfaces for Krylov subspace iterative methods. +! +Module psb_c_inner_krylov_mod + + use psb_base_inner_krylov_mod + + interface psb_init_conv + module procedure psb_c_init_conv + end interface + + interface psb_check_conv + module procedure psb_c_check_conv + end interface + + +contains + subroutine psb_c_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + use psb_sparse_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: stopc, trace, itmax + type(psb_c_sparse_mat), intent(in) :: a + complex(psb_spk_), intent(in) :: b(:) + real(psb_spk_), intent(in) :: eps + type(psb_desc_type), intent(in) :: desc_a + type(psb_itconv_type) :: stopdat + integer, intent(out) :: info + + integer :: ictxt, me, np, err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_init_conv' + call psb_erractionsave(err_act) + + + ictxt=psb_cd_get_context(desc_a) + + call psb_info(ictxt, me, np) + + stopdat%controls(:) = 0 + stopdat%values(:) = 0.0d0 + + stopdat%controls(psb_ik_stopc_) = stopc + stopdat%controls(psb_ik_trace_) = trace + stopdat%controls(psb_ik_itmax_) = itmax + + select case(stopdat%controls(psb_ik_stopc_)) + case (1) + stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) + if (info == psb_success_) stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) + + case (2) + stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + + case default + info=psb_err_invalid_istop_ + call psb_errpush(info,name,i_err=(/stopc,0,0,0,0/)) + goto 9999 + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err="Init conv check data") + goto 9999 + end if + + stopdat%values(psb_ik_eps_) = eps + stopdat%values(psb_ik_errnum_) = dzero + stopdat%values(psb_ik_errden_) = done + + if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))& + & call log_header(methdname) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + + end subroutine psb_c_init_conv + + + function psb_c_check_conv(methdname,it,x,r,desc_a,stopdat,info) + use psb_sparse_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: it + complex(psb_spk_), intent(in) :: x(:), r(:) + type(psb_desc_type), intent(in) :: desc_a + type(psb_itconv_type) :: stopdat + logical :: psb_c_check_conv + integer, intent(out) :: info + + integer :: ictxt, me, np, err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_check_conv' + call psb_erractionsave(err_act) + + ictxt = psb_cd_get_context(desc_a) + call psb_info(ictxt,me,np) + psb_c_check_conv = .false. + + select case(stopdat%controls(psb_ik_stopc_)) + case(1) + stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) + if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) + stopdat%values(psb_ik_errden_) = & + & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)+stopdat%values(psb_ik_bni_)) + case(2) + stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) + stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + + case default + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err="Control data in stopdat messed up!") + goto 9999 + end select + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + if (stopdat%values(psb_ik_errden_) == dzero) then + psb_c_check_conv = (stopdat%values(psb_ik_errnum_) <= stopdat%values(psb_ik_eps_)) + else + psb_c_check_conv = & + & (stopdat%values(psb_ik_errnum_) <= stopdat%values(psb_ik_eps_)*stopdat%values(psb_ik_errden_)) + end if + + psb_c_check_conv = (psb_c_check_conv.or.(stopdat%controls(psb_ik_itmax_) <= it)) + + if ( (stopdat%controls(psb_ik_trace_) > 0).and.& + & ((mod(it,stopdat%controls(psb_ik_trace_)) == 0).or.psb_c_check_conv)) then + call log_conv(methdname,me,it,1,stopdat%values(psb_ik_errnum_),& + & stopdat%values(psb_ik_errden_),stopdat%values(psb_ik_eps_)) + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + + end function psb_c_check_conv + +end module psb_c_inner_krylov_mod diff --git a/krylov/psb_cbicg.f90 b/krylov/psb_cbicg.f90 index 27dba3b9..57ab7731 100644 --- a/krylov/psb_cbicg.f90 +++ b/krylov/psb_cbicg.f90 @@ -96,7 +96,8 @@ subroutine psb_cbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_cbicg + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ parameters diff --git a/krylov/psb_ccg.f90 b/krylov/psb_ccg.f90 index 6af51fba..0f2426b3 100644 --- a/krylov/psb_ccg.f90 +++ b/krylov/psb_ccg.f90 @@ -98,7 +98,8 @@ subroutine psb_ccg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_ccg + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ Parameters diff --git a/krylov/psb_ccgs.f90 b/krylov/psb_ccgs.f90 index 25890d8b..7252a4b6 100644 --- a/krylov/psb_ccgs.f90 +++ b/krylov/psb_ccgs.f90 @@ -95,7 +95,8 @@ Subroutine psb_ccgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_ccgs + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ parameters diff --git a/krylov/psb_ccgstab.f90 b/krylov/psb_ccgstab.f90 index 3e7dc901..6f9ba0a9 100644 --- a/krylov/psb_ccgstab.f90 +++ b/krylov/psb_ccgstab.f90 @@ -96,7 +96,8 @@ subroutine psb_ccgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_ccgstab + use psb_inner_krylov_mod + use psb_krylov_mod Implicit None !!$ parameters Type(psb_c_sparse_mat), Intent(in) :: a diff --git a/krylov/psb_ccgstabl.f90 b/krylov/psb_ccgstabl.f90 index 3fe024e7..11ae8e2e 100644 --- a/krylov/psb_ccgstabl.f90 +++ b/krylov/psb_ccgstabl.f90 @@ -106,7 +106,8 @@ Subroutine psb_ccgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_ccgstabl + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ parameters diff --git a/krylov/psb_ckrylov.f90 b/krylov/psb_ckrylov.f90 new file mode 100644 index 00000000..01caa4ec --- /dev/null +++ b/krylov/psb_ckrylov.f90 @@ -0,0 +1,246 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.2 +!!$ (C) Copyright 2006/2007/2008 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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_krylov_mod.f90 +! Interfaces for Krylov subspace iterative methods. +! + + ! + ! Subroutine: psb_ckrylov + ! + ! Front-end for the Krylov subspace iterations, complexversion + ! + ! Arguments: + ! + ! methd - character The specific method; can take the values: + ! CG + ! CGS + ! BICG + ! BICGSTAB + ! BICGSTABL + ! RGMRES + ! + ! a - type(psb_c_sparse_mat) 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 + ! irst - integer(optional) Input: restart parameter for RGMRES and + ! BICGSTAB(L) methods + ! istop - integer(optional) Input: stopping criterion, or how + ! to estimate the error. + ! 1: err = |r|/|b| + ! 2: err = |r|/(|a||x|+|b|) + ! where r is the (preconditioned, recursive + ! estimate of) residual + ! +Subroutine psb_ckrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) + use psb_sparse_mod + use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type + use psb_krylov_mod, psb_protect_name => psb_ckrylov + character(len=*) :: method + Type(psb_c_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_cprec_type), intent(in) :: prec + complex(psb_spk_), Intent(in) :: b(:) + complex(psb_spk_), Intent(inout) :: x(:) + Real(psb_spk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err + interface + subroutine psb_ccg(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_cprec_type + type(psb_c_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + complex(psb_spk_), intent(in) :: b(:) + complex(psb_spk_), intent(inout) :: x(:) + real(psb_spk_), intent(in) :: eps + class(psb_cprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_spk_), optional, intent(out) :: err + end subroutine psb_ccg + subroutine psb_cbicg(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_cprec_type + type(psb_c_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + complex(psb_spk_), intent(in) :: b(:) + complex(psb_spk_), intent(inout) :: x(:) + real(psb_spk_), intent(in) :: eps + class(psb_cprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_spk_), optional, intent(out) :: err + end subroutine psb_cbicg + subroutine psb_ccgstab(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_cprec_type + type(psb_c_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + complex(psb_spk_), intent(in) :: b(:) + complex(psb_spk_), intent(inout) :: x(:) + real(psb_spk_), intent(in) :: eps + class(psb_cprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_spk_), optional, intent(out) :: err + end subroutine psb_ccgstab + Subroutine psb_ccgstabl(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_cprec_type + Type(psb_c_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_cprec_type), intent(in) :: prec + complex(psb_spk_), Intent(in) :: b(:) + complex(psb_spk_), Intent(inout) :: x(:) + Real(psb_spk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err + end subroutine psb_ccgstabl + Subroutine psb_crgmres(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_cprec_type + Type(psb_c_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_cprec_type), intent(in) :: prec + complex(psb_spk_), Intent(in) :: b(:) + complex(psb_spk_), Intent(inout) :: x(:) + Real(psb_spk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err + end subroutine psb_crgmres + subroutine psb_ccgs(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_cprec_type + type(psb_c_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + complex(psb_spk_), intent(in) :: b(:) + complex(psb_spk_), intent(inout) :: x(:) + real(psb_spk_), intent(in) :: eps + class(psb_cprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_spk_), optional, intent(out) :: err + end subroutine psb_ccgs + end interface + + + integer :: ictxt,me,np,err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_krylov' + call psb_erractionsave(err_act) + + + ictxt=psb_cd_get_context(desc_a) + + call psb_info(ictxt, me, np) + + + select case(psb_toupper(method)) + case('CG') + call psb_ccg(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('CGS') + call psb_ccgs(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('BICG') + call psb_cbicg(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('BICGSTAB') + call psb_ccgstab(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,istop) + case('RGMRES') + call psb_crgmres(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,irst,istop) + case('BICGSTABL') + call psb_ccgstabl(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + case default + if (me == 0) write(0,*) trim(name),': Warning: Unknown method ',method,& + & ', defaulting to BiCGSTAB' + call psb_ccgstab(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + end select + + if(info /= psb_success_) then + 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 == psb_act_abort_) then + call psb_error(ictxt) + return + end if + +end subroutine psb_ckrylov + diff --git a/krylov/psb_crgmres.f90 b/krylov/psb_crgmres.f90 index 43e03f16..9fe4ee5a 100644 --- a/krylov/psb_crgmres.f90 +++ b/krylov/psb_crgmres.f90 @@ -108,7 +108,8 @@ Subroutine psb_crgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_crgmres + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ Parameters diff --git a/krylov/psb_d_inner_krylov_mod.f90 b/krylov/psb_d_inner_krylov_mod.f90 new file mode 100644 index 00000000..56679b80 --- /dev/null +++ b/krylov/psb_d_inner_krylov_mod.f90 @@ -0,0 +1,193 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.2 +!!$ (C) Copyright 2006/2007/2008 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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_krylov_mod.f90 +! Interfaces for Krylov subspace iterative methods. +! +Module psb_d_inner_krylov_mod + + use psb_base_inner_krylov_mod + + interface psb_init_conv + module procedure psb_d_init_conv + end interface + + interface psb_check_conv + module procedure psb_d_check_conv + end interface + + + +contains + + subroutine psb_d_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + use psb_sparse_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: stopc, trace,itmax + type(psb_d_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: b(:), eps + type(psb_desc_type), intent(in) :: desc_a + type(psb_itconv_type) :: stopdat + integer, intent(out) :: info + + integer :: ictxt, me, np, err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_init_conv' + call psb_erractionsave(err_act) + + + ictxt=psb_cd_get_context(desc_a) + + call psb_info(ictxt, me, np) + + stopdat%controls(:) = 0 + stopdat%values(:) = 0.0d0 + + stopdat%controls(psb_ik_stopc_) = stopc + stopdat%controls(psb_ik_trace_) = trace + stopdat%controls(psb_ik_itmax_) = itmax + + select case(stopdat%controls(psb_ik_stopc_)) + case (1) + stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) + if (info == psb_success_) stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) + + case (2) + stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + + case default + info=psb_err_invalid_istop_ + call psb_errpush(info,name,i_err=(/stopc,0,0,0,0/)) + goto 9999 + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err="Init conv check data") + goto 9999 + end if + + stopdat%values(psb_ik_eps_) = eps + stopdat%values(psb_ik_errnum_) = dzero + stopdat%values(psb_ik_errden_) = done + + if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))& + & call log_header(methdname) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + + end subroutine psb_d_init_conv + + function psb_d_check_conv(methdname,it,x,r,desc_a,stopdat,info) + use psb_sparse_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: it + real(psb_dpk_), intent(in) :: x(:), r(:) + type(psb_desc_type), intent(in) :: desc_a + type(psb_itconv_type) :: stopdat + logical :: psb_d_check_conv + integer, intent(out) :: info + + integer :: ictxt, me, np, err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_check_conv' + call psb_erractionsave(err_act) + + ictxt = psb_cd_get_context(desc_a) + call psb_info(ictxt,me,np) + + psb_d_check_conv = .false. + + select case(stopdat%controls(psb_ik_stopc_)) + case(1) + stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) + if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) + stopdat%values(psb_ik_errden_) = & + & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)+stopdat%values(psb_ik_bni_)) + case(2) + stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) + stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + + case default + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err="Control data in stopdat messed up!") + goto 9999 + end select + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + if (stopdat%values(psb_ik_errden_) == dzero) then + psb_d_check_conv = (stopdat%values(psb_ik_errnum_) <= stopdat%values(psb_ik_eps_)) + else + psb_d_check_conv = & + & (stopdat%values(psb_ik_errnum_) <= stopdat%values(psb_ik_eps_)*stopdat%values(psb_ik_errden_)) + end if + + psb_d_check_conv = (psb_d_check_conv.or.(stopdat%controls(psb_ik_itmax_) <= it)) + + if ( (stopdat%controls(psb_ik_trace_) > 0).and.& + & ((mod(it,stopdat%controls(psb_ik_trace_)) == 0).or.psb_d_check_conv)) then + call log_conv(methdname,me,it,1,stopdat%values(psb_ik_errnum_),& + & stopdat%values(psb_ik_errden_),stopdat%values(psb_ik_eps_)) + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + + end function psb_d_check_conv + + +end module psb_d_inner_krylov_mod diff --git a/krylov/psb_dbicg.f90 b/krylov/psb_dbicg.f90 index 3d0e214b..20373fda 100644 --- a/krylov/psb_dbicg.f90 +++ b/krylov/psb_dbicg.f90 @@ -97,7 +97,8 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_dbicg + use psb_inner_krylov_mod + use psb_krylov_mod implicit none type(psb_d_sparse_mat), intent(in) :: a diff --git a/krylov/psb_dcg.F90 b/krylov/psb_dcg.F90 index 3992ab40..46492a6b 100644 --- a/krylov/psb_dcg.F90 +++ b/krylov/psb_dcg.F90 @@ -98,7 +98,8 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,cond) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_dcg + use psb_inner_krylov_mod + use psb_krylov_mod implicit none type(psb_d_sparse_mat), intent(in) :: a diff --git a/krylov/psb_dcgs.f90 b/krylov/psb_dcgs.f90 index 34481edb..e88b8276 100644 --- a/krylov/psb_dcgs.f90 +++ b/krylov/psb_dcgs.f90 @@ -97,7 +97,8 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_dcgs + use psb_inner_krylov_mod + use psb_krylov_mod implicit none type(psb_d_sparse_mat), intent(in) :: a diff --git a/krylov/psb_dcgstab.F90 b/krylov/psb_dcgstab.F90 index f7352a12..3e0276a3 100644 --- a/krylov/psb_dcgstab.F90 +++ b/krylov/psb_dcgstab.F90 @@ -97,7 +97,8 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_dcgstab + use psb_inner_krylov_mod + use psb_krylov_mod implicit none type(psb_d_sparse_mat), intent(in) :: a diff --git a/krylov/psb_dcgstabl.f90 b/krylov/psb_dcgstabl.f90 index 4d5acbbe..956fd328 100644 --- a/krylov/psb_dcgstabl.f90 +++ b/krylov/psb_dcgstabl.f90 @@ -106,7 +106,8 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_dcgstabl + use psb_inner_krylov_mod + use psb_krylov_mod implicit none type(psb_d_sparse_mat), intent(in) :: a diff --git a/krylov/psb_dkrylov.f90 b/krylov/psb_dkrylov.f90 new file mode 100644 index 00000000..43373d32 --- /dev/null +++ b/krylov/psb_dkrylov.f90 @@ -0,0 +1,244 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.2 +!!$ (C) Copyright 2006/2007/2008 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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_krylov_mod.f90 +! Interfaces for Krylov subspace iterative methods. + + ! Subroutine: psb_dkrylov + ! + ! Front-end for the Krylov subspace iterations, real version + ! + ! Arguments: + ! + ! methd - character The specific method; can take the values: + ! CG + ! CGS + ! BICG + ! BICGSTAB + ! BICGSTABL + ! RGMRES + ! + ! a - type(psb_d_sparse_mat) 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 + ! irst - integer(optional) Input: restart parameter for RGMRES and + ! BICGSTAB(L) methods + ! istop - integer(optional) Input: stopping criterion, or how + ! to estimate the error. + ! 1: err = |r|/|b| + ! 2: err = |r|/(|a||x|+|b|) + ! where r is the (preconditioned, recursive + ! estimate of) residual + ! + +Subroutine psb_dkrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop,cond) + + use psb_sparse_mod + use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type + use psb_krylov_mod, psb_protect_name => psb_dkrylov + + character(len=*) :: method + Type(psb_d_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_dprec_type), intent(in) :: prec + Real(psb_dpk_), Intent(in) :: b(:) + Real(psb_dpk_), Intent(inout) :: x(:) + Real(psb_dpk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err,cond + + interface + subroutine psb_dcg(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop,cond) + use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_dprec_type + type(psb_d_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: b(:) + real(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(in) :: eps + class(psb_dprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_dpk_), optional, intent(out) :: err,cond + end subroutine psb_dcg + subroutine psb_dbicg(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_dprec_type + type(psb_d_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: b(:) + real(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(in) :: eps + class(psb_dprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_dpk_), optional, intent(out) :: err + end subroutine psb_dbicg + subroutine psb_dcgstab(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_dprec_type + type(psb_d_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: b(:) + real(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(in) :: eps + class(psb_dprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_dpk_), optional, intent(out) :: err + end subroutine psb_dcgstab + Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err, itrace,irst,istop) + use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_dprec_type + Type(psb_d_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_dprec_type), intent(in) :: prec + Real(psb_dpk_), Intent(in) :: b(:) + Real(psb_dpk_), Intent(inout) :: x(:) + Real(psb_dpk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err + end subroutine psb_dcgstabl + Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_dprec_type + Type(psb_d_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_dprec_type), intent(in) :: prec + Real(psb_dpk_), Intent(in) :: b(:) + Real(psb_dpk_), Intent(inout) :: x(:) + Real(psb_dpk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err + end subroutine psb_drgmres + subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_dprec_type + type(psb_d_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + class(psb_dprec_type), intent(in) :: prec + real(psb_dpk_), intent(in) :: b(:) + real(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(in) :: eps + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_dpk_), optional, intent(out) :: err + end subroutine psb_dcgs + end interface + integer :: ictxt,me,np,err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_krylov' + call psb_erractionsave(err_act) + + + ictxt=psb_cd_get_context(desc_a) + + call psb_info(ictxt, me, np) + + select case(psb_toupper(method)) + case('CG') + call psb_dcg(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop,cond) + case('CGS') + call psb_dcgs(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('BICG') + call psb_dbicg(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('BICGSTAB') + call psb_dcgstab(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('RGMRES') + call psb_drgmres(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + case('BICGSTABL') + call psb_dcgstabl(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + case default + if (me == 0) write(0,*) trim(name),': Warning: Unknown method ',method,& + & ', defaulting to BiCGSTAB' + call psb_dcgstab(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + end select + + if(info /= psb_success_) then + 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 == psb_act_abort_) then + call psb_error(ictxt) + return + end if + +end subroutine psb_dkrylov + diff --git a/krylov/psb_drgmres.f90 b/krylov/psb_drgmres.f90 index 5ae68ae4..be20ae23 100644 --- a/krylov/psb_drgmres.f90 +++ b/krylov/psb_drgmres.f90 @@ -109,7 +109,8 @@ subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_drgmres + use psb_inner_krylov_mod + use psb_krylov_mod implicit none type(psb_d_sparse_mat), intent(in) :: a diff --git a/krylov/psb_inner_krylov_mod.f90 b/krylov/psb_inner_krylov_mod.f90 new file mode 100644 index 00000000..dc76279c --- /dev/null +++ b/krylov/psb_inner_krylov_mod.f90 @@ -0,0 +1,41 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.2 +!!$ (C) Copyright 2006/2007/2008 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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_krylov_mod.f90 +! Interfaces for Krylov subspace iterative methods. +! +module psb_inner_krylov_mod + use psb_s_inner_krylov_mod + use psb_d_inner_krylov_mod + use psb_c_inner_krylov_mod + use psb_z_inner_krylov_mod +end module psb_inner_krylov_mod diff --git a/krylov/psb_krylov_mod.f90 b/krylov/psb_krylov_mod.f90 index 43250990..2f12d6f7 100644 --- a/krylov/psb_krylov_mod.f90 +++ b/krylov/psb_krylov_mod.f90 @@ -34,1554 +34,78 @@ ! Interfaces for Krylov subspace iterative methods. ! Module psb_krylov_mod - + use psb_const_mod public interface psb_krylov - module procedure psb_skrylov, psb_dkrylov, psb_ckrylov, psb_zkrylov - end interface - - interface psb_cg - subroutine psb_scg(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop,cond) - use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_sprec_type - type(psb_s_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(psb_spk_), intent(in) :: b(:) - real(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(in) :: eps - class(psb_sprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_spk_), optional, intent(out) :: err,cond - end subroutine psb_scg - subroutine psb_dcg(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop,cond) - use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_dprec_type - type(psb_d_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(psb_dpk_), intent(in) :: b(:) - real(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(in) :: eps - class(psb_dprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_dpk_), optional, intent(out) :: err,cond - end subroutine psb_dcg - subroutine psb_ccg(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_cprec_type - type(psb_c_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - complex(psb_spk_), intent(in) :: b(:) - complex(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(in) :: eps - class(psb_cprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_spk_), optional, intent(out) :: err - end subroutine psb_ccg - subroutine psb_zcg(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_zprec_type - type(psb_z_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - complex(psb_dpk_), intent(in) :: b(:) - complex(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(in) :: eps - class(psb_zprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_dpk_), optional, intent(out) :: err - end subroutine psb_zcg - end interface - - interface psb_bicg - subroutine psb_sbicg(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_sprec_type - type(psb_s_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(psb_spk_), intent(in) :: b(:) - real(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(in) :: eps - class(psb_sprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_spk_), optional, intent(out) :: err - end subroutine psb_sbicg - subroutine psb_dbicg(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_dprec_type - type(psb_d_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(psb_dpk_), intent(in) :: b(:) - real(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(in) :: eps - class(psb_dprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_dpk_), optional, intent(out) :: err - end subroutine psb_dbicg - subroutine psb_cbicg(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_cprec_type - type(psb_c_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - complex(psb_spk_), intent(in) :: b(:) - complex(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(in) :: eps - class(psb_cprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_spk_), optional, intent(out) :: err - end subroutine psb_cbicg - subroutine psb_zbicg(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_zprec_type - type(psb_z_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - complex(psb_dpk_), intent(in) :: b(:) - complex(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(in) :: eps - class(psb_zprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_dpk_), optional, intent(out) :: err - end subroutine psb_zbicg - end interface - - interface psb_bicgstab - subroutine psb_scgstab(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_sprec_type - type(psb_s_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(psb_spk_), intent(in) :: b(:) - real(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(in) :: eps - class(psb_sprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_spk_), optional, intent(out) :: err - end subroutine psb_scgstab - subroutine psb_dcgstab(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_dprec_type - type(psb_d_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(psb_dpk_), intent(in) :: b(:) - real(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(in) :: eps - class(psb_dprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_dpk_), optional, intent(out) :: err - end subroutine psb_dcgstab - subroutine psb_ccgstab(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_cprec_type - type(psb_c_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - complex(psb_spk_), intent(in) :: b(:) - complex(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(in) :: eps - class(psb_cprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_spk_), optional, intent(out) :: err - end subroutine psb_ccgstab - subroutine psb_zcgstab(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_zprec_type - type(psb_z_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - complex(psb_dpk_), intent(in) :: b(:) - complex(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(in) :: eps - class(psb_zprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_dpk_), optional, intent(out) :: err - end subroutine psb_zcgstab - end interface - - interface psb_bicgstabl - Subroutine psb_scgstabl(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,irst,istop) + + Subroutine psb_skrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop,cond) use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_sprec_type + use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type + + character(len=*) :: method Type(psb_s_sparse_mat), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(psb_sprec_type), intent(in) :: prec - Real(psb_spk_), Intent(in) :: b(:) - Real(psb_spk_), Intent(inout) :: x(:) - Real(psb_spk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(psb_spk_), Optional, Intent(out) :: err - end subroutine psb_scgstabl - Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,irst,istop) - use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_dprec_type - type(psb_d_sparse_mat), intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(psb_dprec_type), intent(in) :: prec - Real(psb_dpk_), Intent(in) :: b(:) - Real(psb_dpk_), Intent(inout) :: x(:) - Real(psb_dpk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(psb_dpk_), Optional, Intent(out) :: err - end subroutine psb_dcgstabl - Subroutine psb_ccgstabl(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_sprec_type), intent(in) :: prec + Real(psb_spk_), Intent(in) :: b(:) + Real(psb_spk_), Intent(inout) :: x(:) + Real(psb_spk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err,cond + end Subroutine psb_skrylov + Subroutine psb_ckrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_cprec_type + use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type + character(len=*) :: method Type(psb_c_sparse_mat), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a + Type(psb_desc_type), Intent(in) :: desc_a class(psb_cprec_type), intent(in) :: prec - complex(psb_spk_), Intent(in) :: b(:) - complex(psb_spk_), Intent(inout) :: x(:) - Real(psb_spk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(psb_spk_), Optional, Intent(out) :: err - end subroutine psb_ccgstabl - Subroutine psb_zcgstabl(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_zprec_type - Type(psb_z_sparse_mat), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(psb_zprec_type), intent(in) :: prec - complex(psb_dpk_), Intent(in) :: b(:) - complex(psb_dpk_), Intent(inout) :: x(:) - Real(psb_dpk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(psb_dpk_), Optional, Intent(out) :: err - end subroutine psb_zcgstabl - end interface - - interface psb_rgmres - Subroutine psb_srgmres(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_sprec_type - Type(psb_s_sparse_mat), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(psb_sprec_type), intent(in) :: prec - Real(psb_spk_), Intent(in) :: b(:) - Real(psb_spk_), Intent(inout) :: x(:) - Real(psb_spk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter + complex(psb_spk_), Intent(in) :: b(:) + complex(psb_spk_), Intent(inout) :: x(:) + Real(psb_spk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter Real(psb_spk_), Optional, Intent(out) :: err - end subroutine psb_srgmres - Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) + end Subroutine psb_ckrylov + Subroutine psb_dkrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop,cond) + use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_dprec_type - type(psb_d_sparse_mat), intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a + use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type + + character(len=*) :: method + Type(psb_d_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a class(psb_dprec_type), intent(in) :: prec - Real(psb_dpk_), Intent(in) :: b(:) - Real(psb_dpk_), Intent(inout) :: x(:) - Real(psb_dpk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(psb_dpk_), Optional, Intent(out) :: err - end subroutine psb_drgmres - Subroutine psb_crgmres(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_cprec_type - Type(psb_c_sparse_mat), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(psb_cprec_type), intent(in) :: prec - complex(psb_spk_), Intent(in) :: b(:) - complex(psb_spk_), Intent(inout) :: x(:) - Real(psb_spk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(psb_spk_), Optional, Intent(out) :: err - end subroutine psb_crgmres - Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) + Real(psb_dpk_), Intent(in) :: b(:) + Real(psb_dpk_), Intent(inout) :: x(:) + Real(psb_dpk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err,cond + + end Subroutine psb_dkrylov + Subroutine psb_zkrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_zprec_type + use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type + character(len=*) :: method Type(psb_z_sparse_mat), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a + Type(psb_desc_type), Intent(in) :: desc_a class(psb_zprec_type), intent(in) :: prec - complex(psb_dpk_), Intent(in) :: b(:) - complex(psb_dpk_), Intent(inout) :: x(:) - Real(psb_dpk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter + complex(psb_dpk_), Intent(in) :: b(:) + complex(psb_dpk_), Intent(inout) :: x(:) + Real(psb_dpk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter Real(psb_dpk_), Optional, Intent(out) :: err - end subroutine psb_zrgmres - end interface - - interface psb_cgs - subroutine psb_scgs(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_sprec_type - type(psb_s_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - class(psb_sprec_type), intent(in) :: prec - real(psb_spk_), intent(in) :: b(:) - real(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(in) :: eps - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_spk_), optional, intent(out) :: err - end subroutine psb_scgs - subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_dprec_type - type(psb_d_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - class(psb_dprec_type), intent(in) :: prec - real(psb_dpk_), intent(in) :: b(:) - real(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(in) :: eps - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_dpk_), optional, intent(out) :: err - end subroutine psb_dcgs - subroutine psb_ccgs(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_c_sparse_mat, psb_spk_ - use psb_prec_mod, only : psb_cprec_type - type(psb_c_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - complex(psb_spk_), intent(in) :: b(:) - complex(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(in) :: eps - class(psb_cprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_spk_), optional, intent(out) :: err - end subroutine psb_ccgs - subroutine psb_zcgs(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ - use psb_prec_mod, only : psb_zprec_type - type(psb_z_sparse_mat), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - complex(psb_dpk_), intent(in) :: b(:) - complex(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(in) :: eps - class(psb_zprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(psb_dpk_), optional, intent(out) :: err - end subroutine psb_zcgs - end interface - + end Subroutine psb_zkrylov - - interface psb_init_conv - module procedure psb_s_init_conv, psb_d_init_conv,& - & psb_c_init_conv, psb_z_init_conv end interface - interface psb_check_conv - module procedure psb_s_check_conv, psb_d_check_conv,& - & psb_c_check_conv, psb_z_check_conv - end interface - - interface psb_end_conv - module procedure psb_d_end_conv - end interface - - integer, parameter, private :: bni_=1, rni_=2, ani_=3, xni_=4, bn2_=5, xn2_=6 - integer, parameter, private :: errnum_=7, errden_=8, eps_=9, rn2_=10 - integer, parameter, private :: stopc_=1, trace_=2, itmax_=3 - integer, parameter, private :: ivsz_=16 - type psb_itconv_type - private - integer :: controls(ivsz_) - real(psb_dpk_) :: values(ivsz_) - end type psb_itconv_type - -contains - - ! Subroutine: psb_skrylov - ! - ! Front-end for the Krylov subspace iterations, real version - ! - ! Arguments: - ! - ! methd - character The specific method; can take the values: - ! CG - ! CGS - ! BICG - ! BICGSTAB - ! BICGSTABL - ! RGMRES - ! - ! a - type(psb_s_sparse_mat) 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 - ! irst - integer(optional) Input: restart parameter for RGMRES and - ! BICGSTAB(L) methods - ! istop - integer(optional) Input: stopping criterion, or how - ! to estimate the error. - ! 1: err = |r|/|b| - ! 2: err = |r|/(|a||x|+|b|) - ! where r is the (preconditioned, recursive - ! estimate of) residual - ! - - Subroutine psb_skrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop,cond) - - use psb_sparse_mod - use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type - - character(len=*) :: method - Type(psb_s_sparse_mat), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(psb_sprec_type), intent(in) :: prec - Real(psb_spk_), Intent(in) :: b(:) - Real(psb_spk_), Intent(inout) :: x(:) - Real(psb_spk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(psb_spk_), Optional, Intent(out) :: err,cond - - integer :: ictxt,me,np,err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_krylov' - call psb_erractionsave(err_act) - - - ictxt=psb_cd_get_context(desc_a) - - call psb_info(ictxt, me, np) - - select case(psb_toupper(method)) - case('CG') - call psb_cg(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop,cond) - case('CGS') - call psb_cgs(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('BICG') - call psb_bicg(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('BICGSTAB') - call psb_bicgstab(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('RGMRES') - call psb_rgmres(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - case('BICGSTABL') - call psb_bicgstabl(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - case default - if (me == 0) write(0,*) trim(name),': Warning: Unknown method ',method,& - & ', defaulting to BiCGSTAB' - call psb_bicgstab(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - end select - - if(info /= psb_success_) then - 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 == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end subroutine psb_skrylov - - - ! Subroutine: psb_dkrylov - ! - ! Front-end for the Krylov subspace iterations, real version - ! - ! Arguments: - ! - ! methd - character The specific method; can take the values: - ! CG - ! CGS - ! BICG - ! BICGSTAB - ! BICGSTABL - ! RGMRES - ! - ! a - type(psb_d_sparse_mat) 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 - ! irst - integer(optional) Input: restart parameter for RGMRES and - ! BICGSTAB(L) methods - ! istop - integer(optional) Input: stopping criterion, or how - ! to estimate the error. - ! 1: err = |r|/|b| - ! 2: err = |r|/(|a||x|+|b|) - ! where r is the (preconditioned, recursive - ! estimate of) residual - ! - - Subroutine psb_dkrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop,cond) - - use psb_sparse_mod - use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type - - type(psb_d_sparse_mat), intent(in) :: a - - character(len=*) :: method - Type(psb_desc_type), Intent(in) :: desc_a - class(psb_dprec_type), intent(in) :: prec - Real(psb_dpk_), Intent(in) :: b(:) - Real(psb_dpk_), Intent(inout) :: x(:) - Real(psb_dpk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(psb_dpk_), Optional, Intent(out) :: err,cond - - integer :: ictxt,me,np,err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_krylov' - call psb_erractionsave(err_act) - - - ictxt=psb_cd_get_context(desc_a) - - call psb_info(ictxt, me, np) - - select case(psb_toupper(method)) - case('CG') - call psb_cg(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop,cond) - case('CGS') - call psb_cgs(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('BICG') - call psb_bicg(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('BICGSTAB') - call psb_bicgstab(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('RGMRES') - call psb_rgmres(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - case('BICGSTABL') - call psb_bicgstabl(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - case default - if (me == 0) write(0,*) trim(name),': Warning: Unknown method ',method,& - & ', defaulting to BiCGSTAB' - call psb_bicgstab(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - end select - - if(info /= psb_success_) then - 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 == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end subroutine psb_dkrylov - - - ! - ! Subroutine: psb_ckrylov - ! - ! Front-end for the Krylov subspace iterations, complexversion - ! - ! Arguments: - ! - ! methd - character The specific method; can take the values: - ! CG - ! CGS - ! BICG - ! BICGSTAB - ! BICGSTABL - ! RGMRES - ! - ! a - type(psb_c_sparse_mat) 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 - ! irst - integer(optional) Input: restart parameter for RGMRES and - ! BICGSTAB(L) methods - ! istop - integer(optional) Input: stopping criterion, or how - ! to estimate the error. - ! 1: err = |r|/|b| - ! 2: err = |r|/(|a||x|+|b|) - ! where r is the (preconditioned, recursive - ! estimate of) residual - ! - Subroutine psb_ckrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) - use psb_sparse_mod - use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type - character(len=*) :: method - Type(psb_c_sparse_mat), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(psb_cprec_type), intent(in) :: prec - complex(psb_spk_), Intent(in) :: b(:) - complex(psb_spk_), Intent(inout) :: x(:) - Real(psb_spk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(psb_spk_), Optional, Intent(out) :: err - - integer :: ictxt,me,np,err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_krylov' - call psb_erractionsave(err_act) - - - ictxt=psb_cd_get_context(desc_a) - - call psb_info(ictxt, me, np) - - - select case(psb_toupper(method)) - case('CG') - call psb_cg(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('CGS') - call psb_cgs(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('BICG') - call psb_bicg(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('BICGSTAB') - call psb_bicgstab(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop) - case('RGMRES') - call psb_rgmres(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop) - case('BICGSTABL') - call psb_bicgstabl(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - case default - if (me == 0) write(0,*) trim(name),': Warning: Unknown method ',method,& - & ', defaulting to BiCGSTAB' - call psb_bicgstab(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - end select - - if(info /= psb_success_) then - 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 == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end subroutine psb_ckrylov - - ! - ! Subroutine: psb_zkrylov - ! - ! Front-end for the Krylov subspace iterations, complexversion - ! - ! Arguments: - ! - ! methd - character The specific method; can take the values: - ! CG - ! CGS - ! BICG - ! BICGSTAB - ! BICGSTABL - ! RGMRES - ! - ! a - type(psb_z_sparse_mat) 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 - ! irst - integer(optional) Input: restart parameter for RGMRES and - ! BICGSTAB(L) methods - ! istop - integer(optional) Input: stopping criterion, or how - ! to estimate the error. - ! 1: err = |r|/|b| - ! 2: err = |r|/(|a||x|+|b|) - ! where r is the (preconditioned, recursive - ! estimate of) residual - ! - Subroutine psb_zkrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) - use psb_sparse_mod - use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type - character(len=*) :: method - Type(psb_z_sparse_mat), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(psb_zprec_type), intent(in) :: prec - complex(psb_dpk_), Intent(in) :: b(:) - complex(psb_dpk_), Intent(inout) :: x(:) - Real(psb_dpk_), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(psb_dpk_), Optional, Intent(out) :: err - - integer :: ictxt,me,np,err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_krylov' - call psb_erractionsave(err_act) - - - ictxt=psb_cd_get_context(desc_a) - - call psb_info(ictxt, me, np) - - - select case(psb_toupper(method)) - case('CG') - call psb_cg(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('CGS') - call psb_cgs(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('BICG') - call psb_bicg(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - case('BICGSTAB') - call psb_bicgstab(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,istop) - case('RGMRES') - call psb_rgmres(a,prec,b,x,eps,desc_a,info,& - & itmax,iter,err,itrace,irst,istop) - case('BICGSTABL') - call psb_bicgstabl(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - case default - if (me == 0) write(0,*) trim(name),': Warning: Unknown method ',method,& - & ', defaulting to BiCGSTAB' - call psb_bicgstab(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - end select - - if(info /= psb_success_) then - 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 == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end subroutine psb_zkrylov - - subroutine log_header(methdname) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - character(len=*), parameter :: fmt='(a18,1x,a4,3(2x,a15))' - integer, parameter :: outlen=18 - character(len=len(methdname)) :: mname - character(len=outlen) :: outname - - mname = adjustl(trim(methdname)) - write(outname,'(a)') mname(1:min(len_trim(mname),outlen-1))//':' - write(*,fmt) adjustl(outname),'Iteration','Error Estimate','Tolerance' - - end subroutine log_header - - - subroutine log_conv(methdname,me,itx,itrace,errnum,errden,eps) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - integer, intent(in) :: me, itx, itrace - real(psb_dpk_), intent(in) :: errnum, errden, eps - character(len=*), parameter :: fmt='(a18,1x,i4,3(2x,es15.9))' - integer, parameter :: outlen=18 - character(len=len(methdname)) :: mname - character(len=outlen) :: outname - - if ((mod(itx,itrace) == 0).and.(me == 0)) then - mname = adjustl(trim(methdname)) - write(outname,'(a)') mname(1:min(len_trim(mname),outlen-1))//':' - if (errden > dzero ) then - write(*,fmt) adjustl(outname),itx,errnum/errden,eps - else - write(*,fmt) adjustl(outname),itx,errnum,eps - end if - endif - - end subroutine log_conv - - subroutine log_end(methdname,me,it,errnum,errden,eps,err,iter) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - integer, intent(in) :: me, it - real(psb_dpk_), intent(in) :: errnum, errden, eps - real(psb_dpk_), optional, intent(out) :: err - integer, optional, intent(out) :: iter - - character(len=*), parameter :: fmt='(a,2x,es15.9,1x,a,1x,i4,1x,a)' - character(len=*), parameter :: fmt1='(a,3(2x,es15.9))' - - if (errden == dzero) then - if (errnum > eps) then - if (me == 0) then - write(*,fmt) trim(methdname)//' failed to converge to ',eps,& - & ' in ',it,' iterations. ' - write(*,fmt1) 'Last iteration error estimate: ',& - & errnum - end if - end if - if (present(err)) err=errnum - else - if (errnum/errden > eps) then - if (me == 0) then - write(*,fmt) trim(methdname)//' failed to converge to ',eps,& - & ' in ',it,' iterations. ' - write(*,fmt1) 'Last iteration error estimate: ',& - & errnum/errden - end if - endif - if (present(err)) err=errnum/errden - end if - if (present(iter)) iter = it - - end subroutine log_end - - subroutine psb_s_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - integer, intent(in) :: stopc, trace,itmax - type(psb_s_sparse_mat), intent(in) :: a - real(psb_spk_), intent(in) :: b(:), eps - type(psb_desc_type), intent(in) :: desc_a - type(psb_itconv_type) :: stopdat - integer, intent(out) :: info - - integer :: ictxt, me, np, err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_init_conv' - call psb_erractionsave(err_act) - - - ictxt=psb_cd_get_context(desc_a) - - call psb_info(ictxt, me, np) - - stopdat%controls(:) = 0 - stopdat%values(:) = 0.0d0 - - stopdat%controls(stopc_) = stopc - stopdat%controls(trace_) = trace - stopdat%controls(itmax_) = itmax - - select case(stopdat%controls(stopc_)) - case (1) - stopdat%values(ani_) = psb_spnrmi(a,desc_a,info) - if (info == psb_success_) stopdat%values(bni_) = psb_geamax(b,desc_a,info) - - case (2) - stopdat%values(bn2_) = psb_genrm2(b,desc_a,info) - - case default - info=psb_err_invalid_istop_ - call psb_errpush(info,name,i_err=(/stopc,0,0,0,0/)) - goto 9999 - end select - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err="Init conv check data") - goto 9999 - end if - - stopdat%values(eps_) = eps - stopdat%values(errnum_) = dzero - stopdat%values(errden_) = done - - if ((stopdat%controls(trace_) > 0).and. (me == 0))& - & call log_header(methdname) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end subroutine psb_s_init_conv - - subroutine psb_d_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) - use psb_sparse_mod - implicit none - type(psb_d_sparse_mat), intent(in) :: a - character(len=*), intent(in) :: methdname - integer, intent(in) :: stopc, trace,itmax - real(psb_dpk_), intent(in) :: b(:), eps - type(psb_desc_type), intent(in) :: desc_a - type(psb_itconv_type) :: stopdat - integer, intent(out) :: info - - integer :: ictxt, me, np, err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_init_conv' - call psb_erractionsave(err_act) - - - ictxt=psb_cd_get_context(desc_a) - - call psb_info(ictxt, me, np) - - stopdat%controls(:) = 0 - stopdat%values(:) = 0.0d0 - - stopdat%controls(stopc_) = stopc - stopdat%controls(trace_) = trace - stopdat%controls(itmax_) = itmax - - select case(stopdat%controls(stopc_)) - case (1) - stopdat%values(ani_) = psb_spnrmi(a,desc_a,info) - if (info == psb_success_) stopdat%values(bni_) = psb_geamax(b,desc_a,info) - - case (2) - stopdat%values(bn2_) = psb_genrm2(b,desc_a,info) - - case default - info=psb_err_invalid_istop_ - call psb_errpush(info,name,i_err=(/stopc,0,0,0,0/)) - goto 9999 - end select - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err="Init conv check data") - goto 9999 - end if - - stopdat%values(eps_) = eps - stopdat%values(errnum_) = dzero - stopdat%values(errden_) = done - - if ((stopdat%controls(trace_) > 0).and. (me == 0))& - & call log_header(methdname) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end subroutine psb_d_init_conv - - subroutine psb_c_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - integer, intent(in) :: stopc, trace, itmax - type(psb_c_sparse_mat), intent(in) :: a - complex(psb_spk_), intent(in) :: b(:) - real(psb_spk_), intent(in) :: eps - type(psb_desc_type), intent(in) :: desc_a - type(psb_itconv_type) :: stopdat - integer, intent(out) :: info - - integer :: ictxt, me, np, err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_init_conv' - call psb_erractionsave(err_act) - - - ictxt=psb_cd_get_context(desc_a) - - call psb_info(ictxt, me, np) - - stopdat%controls(:) = 0 - stopdat%values(:) = 0.0d0 - - stopdat%controls(stopc_) = stopc - stopdat%controls(trace_) = trace - stopdat%controls(itmax_) = itmax - - select case(stopdat%controls(stopc_)) - case (1) - stopdat%values(ani_) = psb_spnrmi(a,desc_a,info) - if (info == psb_success_) stopdat%values(bni_) = psb_geamax(b,desc_a,info) - - case (2) - stopdat%values(bn2_) = psb_genrm2(b,desc_a,info) - - case default - info=psb_err_invalid_istop_ - call psb_errpush(info,name,i_err=(/stopc,0,0,0,0/)) - goto 9999 - end select - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err="Init conv check data") - goto 9999 - end if - - stopdat%values(eps_) = eps - stopdat%values(errnum_) = dzero - stopdat%values(errden_) = done - - if ((stopdat%controls(trace_) > 0).and. (me == 0))& - & call log_header(methdname) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end subroutine psb_c_init_conv - - subroutine psb_z_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - integer, intent(in) :: stopc, trace, itmax - type(psb_z_sparse_mat), intent(in) :: a - complex(psb_dpk_), intent(in) :: b(:) - real(psb_dpk_), intent(in) :: eps - type(psb_desc_type), intent(in) :: desc_a - type(psb_itconv_type) :: stopdat - integer, intent(out) :: info - - integer :: ictxt, me, np, err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_init_conv' - call psb_erractionsave(err_act) - - - ictxt=psb_cd_get_context(desc_a) - - call psb_info(ictxt, me, np) - - stopdat%controls(:) = 0 - stopdat%values(:) = 0.0d0 - - stopdat%controls(stopc_) = stopc - stopdat%controls(trace_) = trace - stopdat%controls(itmax_) = itmax - - select case(stopdat%controls(stopc_)) - case (1) - stopdat%values(ani_) = psb_spnrmi(a,desc_a,info) - if (info == psb_success_) stopdat%values(bni_) = psb_geamax(b,desc_a,info) - - case (2) - stopdat%values(bn2_) = psb_genrm2(b,desc_a,info) - - case default - info=psb_err_invalid_istop_ - call psb_errpush(info,name,i_err=(/stopc,0,0,0,0/)) - goto 9999 - end select - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err="Init conv check data") - goto 9999 - end if - - stopdat%values(eps_) = eps - stopdat%values(errnum_) = dzero - stopdat%values(errden_) = done - - if ((stopdat%controls(trace_) > 0).and. (me == 0))& - & call log_header(methdname) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end subroutine psb_z_init_conv - - - function psb_s_check_conv(methdname,it,x,r,desc_a,stopdat,info) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - integer, intent(in) :: it - real(psb_spk_), intent(in) :: x(:), r(:) - type(psb_desc_type), intent(in) :: desc_a - type(psb_itconv_type) :: stopdat - logical :: psb_s_check_conv - integer, intent(out) :: info - - integer :: ictxt, me, np, err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_check_conv' - call psb_erractionsave(err_act) - - ictxt = psb_cd_get_context(desc_a) - call psb_info(ictxt,me,np) - - psb_s_check_conv = .false. - - select case(stopdat%controls(stopc_)) - case(1) - stopdat%values(rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(xni_) = psb_geamax(x,desc_a,info) - stopdat%values(errnum_) = stopdat%values(rni_) - stopdat%values(errden_) = & - & (stopdat%values(ani_)*stopdat%values(xni_)+stopdat%values(bni_)) - case(2) - stopdat%values(rn2_) = psb_genrm2(r,desc_a,info) - stopdat%values(errnum_) = stopdat%values(rn2_) - stopdat%values(errden_) = stopdat%values(bn2_) - - case default - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err="Control data in stopdat messed up!") - goto 9999 - end select - if (info /= psb_success_) then - info=psb_err_from_subroutine_non_ - call psb_errpush(info,name) - goto 9999 - end if - - if (stopdat%values(errden_) == dzero) then - psb_s_check_conv = (stopdat%values(errnum_) <= stopdat%values(eps_)) - else - psb_s_check_conv = & - & (stopdat%values(errnum_) <= stopdat%values(eps_)*stopdat%values(errden_)) - end if - - psb_s_check_conv = (psb_s_check_conv.or.(stopdat%controls(itmax_) <= it)) - - if ( (stopdat%controls(trace_) > 0).and.& - & ((mod(it,stopdat%controls(trace_)) == 0).or.psb_s_check_conv)) then - call log_conv(methdname,me,it,1,stopdat%values(errnum_),& - & stopdat%values(errden_),stopdat%values(eps_)) - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end function psb_s_check_conv - - function psb_d_check_conv(methdname,it,x,r,desc_a,stopdat,info) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - integer, intent(in) :: it - real(psb_dpk_), intent(in) :: x(:), r(:) - type(psb_desc_type), intent(in) :: desc_a - type(psb_itconv_type) :: stopdat - logical :: psb_d_check_conv - integer, intent(out) :: info - - integer :: ictxt, me, np, err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_check_conv' - call psb_erractionsave(err_act) - - ictxt = psb_cd_get_context(desc_a) - call psb_info(ictxt,me,np) - - psb_d_check_conv = .false. - - select case(stopdat%controls(stopc_)) - case(1) - stopdat%values(rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(xni_) = psb_geamax(x,desc_a,info) - stopdat%values(errnum_) = stopdat%values(rni_) - stopdat%values(errden_) = & - & (stopdat%values(ani_)*stopdat%values(xni_)+stopdat%values(bni_)) - case(2) - stopdat%values(rn2_) = psb_genrm2(r,desc_a,info) - stopdat%values(errnum_) = stopdat%values(rn2_) - stopdat%values(errden_) = stopdat%values(bn2_) - - case default - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err="Control data in stopdat messed up!") - goto 9999 - end select - if (info /= psb_success_) then - info=psb_err_from_subroutine_non_ - call psb_errpush(info,name) - goto 9999 - end if - - if (stopdat%values(errden_) == dzero) then - psb_d_check_conv = (stopdat%values(errnum_) <= stopdat%values(eps_)) - else - psb_d_check_conv = & - & (stopdat%values(errnum_) <= stopdat%values(eps_)*stopdat%values(errden_)) - end if - - psb_d_check_conv = (psb_d_check_conv.or.(stopdat%controls(itmax_) <= it)) - - if ( (stopdat%controls(trace_) > 0).and.& - & ((mod(it,stopdat%controls(trace_)) == 0).or.psb_d_check_conv)) then - call log_conv(methdname,me,it,1,stopdat%values(errnum_),& - & stopdat%values(errden_),stopdat%values(eps_)) - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end function psb_d_check_conv - - - function psb_c_check_conv(methdname,it,x,r,desc_a,stopdat,info) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - integer, intent(in) :: it - complex(psb_spk_), intent(in) :: x(:), r(:) - type(psb_desc_type), intent(in) :: desc_a - type(psb_itconv_type) :: stopdat - logical :: psb_c_check_conv - integer, intent(out) :: info - - integer :: ictxt, me, np, err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_check_conv' - call psb_erractionsave(err_act) - - ictxt = psb_cd_get_context(desc_a) - call psb_info(ictxt,me,np) - psb_c_check_conv = .false. - - select case(stopdat%controls(stopc_)) - case(1) - stopdat%values(rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(xni_) = psb_geamax(x,desc_a,info) - stopdat%values(errnum_) = stopdat%values(rni_) - stopdat%values(errden_) = & - & (stopdat%values(ani_)*stopdat%values(xni_)+stopdat%values(bni_)) - case(2) - stopdat%values(rn2_) = psb_genrm2(r,desc_a,info) - stopdat%values(errnum_) = stopdat%values(rn2_) - stopdat%values(errden_) = stopdat%values(bn2_) - - case default - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err="Control data in stopdat messed up!") - goto 9999 - end select - if (info /= psb_success_) then - info=psb_err_from_subroutine_non_ - call psb_errpush(info,name) - goto 9999 - end if - - if (stopdat%values(errden_) == dzero) then - psb_c_check_conv = (stopdat%values(errnum_) <= stopdat%values(eps_)) - else - psb_c_check_conv = & - & (stopdat%values(errnum_) <= stopdat%values(eps_)*stopdat%values(errden_)) - end if - - psb_c_check_conv = (psb_c_check_conv.or.(stopdat%controls(itmax_) <= it)) - - if ( (stopdat%controls(trace_) > 0).and.& - & ((mod(it,stopdat%controls(trace_)) == 0).or.psb_c_check_conv)) then - call log_conv(methdname,me,it,1,stopdat%values(errnum_),& - & stopdat%values(errden_),stopdat%values(eps_)) - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end function psb_c_check_conv - - function psb_z_check_conv(methdname,it,x,r,desc_a,stopdat,info) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - integer, intent(in) :: it - complex(psb_dpk_), intent(in) :: x(:), r(:) - type(psb_desc_type), intent(in) :: desc_a - type(psb_itconv_type) :: stopdat - logical :: psb_z_check_conv - integer, intent(out) :: info - - integer :: ictxt, me, np, err_act - character(len=20) :: name - - info = psb_success_ - name = 'psb_check_conv' - call psb_erractionsave(err_act) - - ictxt = psb_cd_get_context(desc_a) - call psb_info(ictxt,me,np) - psb_z_check_conv = .false. - - select case(stopdat%controls(stopc_)) - case(1) - stopdat%values(rni_) = psb_geamax(r,desc_a,info) - if (info == psb_success_) stopdat%values(xni_) = psb_geamax(x,desc_a,info) - stopdat%values(errnum_) = stopdat%values(rni_) - stopdat%values(errden_) = & - & (stopdat%values(ani_)*stopdat%values(xni_)+stopdat%values(bni_)) - case(2) - stopdat%values(rn2_) = psb_genrm2(r,desc_a,info) - stopdat%values(errnum_) = stopdat%values(rn2_) - stopdat%values(errden_) = stopdat%values(bn2_) - - case default - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err="Control data in stopdat messed up!") - goto 9999 - end select - if (info /= psb_success_) then - info=psb_err_from_subroutine_non_ - call psb_errpush(info,name) - goto 9999 - end if - - if (stopdat%values(errden_) == dzero) then - psb_z_check_conv = (stopdat%values(errnum_) <= stopdat%values(eps_)) - else - psb_z_check_conv = & - & (stopdat%values(errnum_) <= stopdat%values(eps_)*stopdat%values(errden_)) - end if - - psb_z_check_conv = (psb_z_check_conv.or.(stopdat%controls(itmax_) <= it)) - - if ( (stopdat%controls(trace_) > 0).and.& - & ((mod(it,stopdat%controls(trace_)) == 0).or.psb_z_check_conv)) then - call log_conv(methdname,me,it,1,stopdat%values(errnum_),& - & stopdat%values(errden_),stopdat%values(eps_)) - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end function psb_z_check_conv - - subroutine psb_d_end_conv(methdname,it,desc_a,stopdat,info,err,iter) - use psb_sparse_mod - implicit none - character(len=*), intent(in) :: methdname - integer, intent(in) :: it - type(psb_desc_type), intent(in) :: desc_a - type(psb_itconv_type) :: stopdat - integer, intent(out) :: info - real(psb_dpk_), optional, intent(out) :: err - integer, optional, intent(out) :: iter - - integer :: ictxt, me, np, err_act - real(psb_dpk_) :: errnum, errden, eps - character(len=20) :: name - - info = psb_success_ - name = 'psb_end_conv' - - ictxt = psb_cd_get_context(desc_a) - call psb_info(ictxt,me,np) - - errnum = stopdat%values(errnum_) - errden = stopdat%values(errden_) - eps = stopdat%values(eps_) - call log_end(methdname,me,it,errnum,errden,eps,err,iter) - - end subroutine psb_d_end_conv end module psb_krylov_mod - - - diff --git a/krylov/psb_s_inner_krylov_mod.f90 b/krylov/psb_s_inner_krylov_mod.f90 new file mode 100644 index 00000000..48363cca --- /dev/null +++ b/krylov/psb_s_inner_krylov_mod.f90 @@ -0,0 +1,192 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.2 +!!$ (C) Copyright 2006/2007/2008 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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_krylov_mod.f90 +! Interfaces for Krylov subspace iterative methods. +! +Module psb_s_inner_krylov_mod + + use psb_base_inner_krylov_mod + + interface psb_init_conv + module procedure psb_s_init_conv + end interface + + interface psb_check_conv + module procedure psb_s_check_conv + end interface + + + +contains + + subroutine psb_s_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + use psb_sparse_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: stopc, trace,itmax + type(psb_s_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: b(:), eps + type(psb_desc_type), intent(in) :: desc_a + type(psb_itconv_type) :: stopdat + integer, intent(out) :: info + + integer :: ictxt, me, np, err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_init_conv' + call psb_erractionsave(err_act) + + + ictxt=psb_cd_get_context(desc_a) + + call psb_info(ictxt, me, np) + + stopdat%controls(:) = 0 + stopdat%values(:) = 0.0d0 + + stopdat%controls(psb_ik_stopc_) = stopc + stopdat%controls(psb_ik_trace_) = trace + stopdat%controls(psb_ik_itmax_) = itmax + + select case(stopdat%controls(psb_ik_stopc_)) + case (1) + stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) + if (info == psb_success_) stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) + + case (2) + stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + + case default + info=psb_err_invalid_istop_ + call psb_errpush(info,name,i_err=(/stopc,0,0,0,0/)) + goto 9999 + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err="Init conv check data") + goto 9999 + end if + + stopdat%values(psb_ik_eps_) = eps + stopdat%values(psb_ik_errnum_) = dzero + stopdat%values(psb_ik_errden_) = done + + if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))& + & call log_header(methdname) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + + end subroutine psb_s_init_conv + + function psb_s_check_conv(methdname,it,x,r,desc_a,stopdat,info) + use psb_sparse_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: it + real(psb_spk_), intent(in) :: x(:), r(:) + type(psb_desc_type), intent(in) :: desc_a + type(psb_itconv_type) :: stopdat + logical :: psb_s_check_conv + integer, intent(out) :: info + + integer :: ictxt, me, np, err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_check_conv' + call psb_erractionsave(err_act) + + ictxt = psb_cd_get_context(desc_a) + call psb_info(ictxt,me,np) + + psb_s_check_conv = .false. + + select case(stopdat%controls(psb_ik_stopc_)) + case(1) + stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) + if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) + stopdat%values(psb_ik_errden_) = & + & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)+stopdat%values(psb_ik_bni_)) + case(2) + stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) + stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + + case default + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err="Control data in stopdat messed up!") + goto 9999 + end select + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + if (stopdat%values(psb_ik_errden_) == dzero) then + psb_s_check_conv = (stopdat%values(psb_ik_errnum_) <= stopdat%values(psb_ik_eps_)) + else + psb_s_check_conv = & + & (stopdat%values(psb_ik_errnum_) <= stopdat%values(psb_ik_eps_)*stopdat%values(psb_ik_errden_)) + end if + + psb_s_check_conv = (psb_s_check_conv.or.(stopdat%controls(psb_ik_itmax_) <= it)) + + if ( (stopdat%controls(psb_ik_trace_) > 0).and.& + & ((mod(it,stopdat%controls(psb_ik_trace_)) == 0).or.psb_s_check_conv)) then + call log_conv(methdname,me,it,1,stopdat%values(psb_ik_errnum_),& + & stopdat%values(psb_ik_errden_),stopdat%values(psb_ik_eps_)) + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + + end function psb_s_check_conv + +end module psb_s_inner_krylov_mod diff --git a/krylov/psb_sbicg.f90 b/krylov/psb_sbicg.f90 index ff9a0141..ec0b924a 100644 --- a/krylov/psb_sbicg.f90 +++ b/krylov/psb_sbicg.f90 @@ -97,7 +97,8 @@ subroutine psb_sbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_sbicg + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ parameters diff --git a/krylov/psb_scg.F90 b/krylov/psb_scg.F90 index 59b01d87..6643490b 100644 --- a/krylov/psb_scg.F90 +++ b/krylov/psb_scg.F90 @@ -98,7 +98,8 @@ subroutine psb_scg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,cond) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_scg + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ Parameters diff --git a/krylov/psb_scgs.f90 b/krylov/psb_scgs.f90 index 55cddad0..31d09977 100644 --- a/krylov/psb_scgs.f90 +++ b/krylov/psb_scgs.f90 @@ -97,7 +97,8 @@ Subroutine psb_scgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_scgs + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ parameters diff --git a/krylov/psb_scgstab.F90 b/krylov/psb_scgstab.F90 index 79edfb44..53b6fc8d 100644 --- a/krylov/psb_scgstab.F90 +++ b/krylov/psb_scgstab.F90 @@ -97,7 +97,8 @@ Subroutine psb_scgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_scgstab + use psb_inner_krylov_mod + use psb_krylov_mod Implicit None !!$ parameters Type(psb_s_sparse_mat), Intent(in) :: a diff --git a/krylov/psb_scgstabl.f90 b/krylov/psb_scgstabl.f90 index 1abf0bf8..7d1e2b42 100644 --- a/krylov/psb_scgstabl.f90 +++ b/krylov/psb_scgstabl.f90 @@ -106,7 +106,8 @@ Subroutine psb_scgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_scgstabl + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ parameters diff --git a/krylov/psb_skrylov.f90 b/krylov/psb_skrylov.f90 new file mode 100644 index 00000000..25ef9112 --- /dev/null +++ b/krylov/psb_skrylov.f90 @@ -0,0 +1,249 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.2 +!!$ (C) Copyright 2006/2007/2008 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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_krylov_mod.f90 +! Interfaces for Krylov subspace iterative methods. + ! Subroutine: psb_skrylov + ! + ! Front-end for the Krylov subspace iterations, real version + ! + ! Arguments: + ! + ! methd - character The specific method; can take the values: + ! CG + ! CGS + ! BICG + ! BICGSTAB + ! BICGSTABL + ! RGMRES + ! + ! a - type(psb_s_sparse_mat) 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 + ! irst - integer(optional) Input: restart parameter for RGMRES and + ! BICGSTAB(L) methods + ! istop - integer(optional) Input: stopping criterion, or how + ! to estimate the error. + ! 1: err = |r|/|b| + ! 2: err = |r|/(|a||x|+|b|) + ! where r is the (preconditioned, recursive + ! estimate of) residual + ! + +Subroutine psb_skrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop,cond) + + use psb_sparse_mod + use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type + use psb_krylov_mod, psb_protect_name => psb_skrylov + character(len=*) :: method + Type(psb_s_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_sprec_type), intent(in) :: prec + Real(psb_spk_), Intent(in) :: b(:) + Real(psb_spk_), Intent(inout) :: x(:) + Real(psb_spk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err,cond + + interface + subroutine psb_scg(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop,cond) + use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_sprec_type + type(psb_s_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_spk_), intent(in) :: b(:) + real(psb_spk_), intent(inout) :: x(:) + real(psb_spk_), intent(in) :: eps + class(psb_sprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_spk_), optional, intent(out) :: err,cond + end subroutine psb_scg + + subroutine psb_sbicg(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_sprec_type + type(psb_s_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_spk_), intent(in) :: b(:) + real(psb_spk_), intent(inout) :: x(:) + real(psb_spk_), intent(in) :: eps + class(psb_sprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_spk_), optional, intent(out) :: err + end subroutine psb_sbicg + + subroutine psb_scgstab(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_sprec_type + type(psb_s_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_spk_), intent(in) :: b(:) + real(psb_spk_), intent(inout) :: x(:) + real(psb_spk_), intent(in) :: eps + class(psb_sprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_spk_), optional, intent(out) :: err + end subroutine psb_scgstab + + Subroutine psb_scgstabl(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err, itrace,irst,istop) + use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_sprec_type + Type(psb_s_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_sprec_type), intent(in) :: prec + Real(psb_spk_), Intent(in) :: b(:) + Real(psb_spk_), Intent(inout) :: x(:) + Real(psb_spk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err + end subroutine psb_scgstabl + Subroutine psb_srgmres(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_sprec_type + Type(psb_s_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_sprec_type), intent(in) :: prec + Real(psb_spk_), Intent(in) :: b(:) + Real(psb_spk_), Intent(inout) :: x(:) + Real(psb_spk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_spk_), Optional, Intent(out) :: err + end subroutine psb_srgmres + subroutine psb_scgs(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_s_sparse_mat, psb_spk_ + use psb_prec_mod, only : psb_sprec_type + type(psb_s_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + class(psb_sprec_type), intent(in) :: prec + real(psb_spk_), intent(in) :: b(:) + real(psb_spk_), intent(inout) :: x(:) + real(psb_spk_), intent(in) :: eps + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_spk_), optional, intent(out) :: err + end subroutine psb_scgs + end interface + + integer :: ictxt,me,np,err_act + character(len=20) :: name + + + info = psb_success_ + name = 'psb_krylov' + call psb_erractionsave(err_act) + + + ictxt=psb_cd_get_context(desc_a) + + call psb_info(ictxt, me, np) + + select case(psb_toupper(method)) + case('CG') + call psb_scg(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop,cond) + case('CGS') + call psb_scgs(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('BICG') + call psb_sbicg(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('BICGSTAB') + call psb_scgstab(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('RGMRES') + call psb_srgmres(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + case('BICGSTABL') + call psb_scgstabl(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + case default + if (me == 0) write(0,*) trim(name),': Warning: Unknown method ',method,& + & ', defaulting to BiCGSTAB' + call psb_scgstab(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + end select + + if(info /= psb_success_) then + 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 == psb_act_abort_) then + call psb_error(ictxt) + return + end if + +end subroutine psb_skrylov + + + diff --git a/krylov/psb_srgmres.f90 b/krylov/psb_srgmres.f90 index 40aa7cd6..da5a1169 100644 --- a/krylov/psb_srgmres.f90 +++ b/krylov/psb_srgmres.f90 @@ -109,7 +109,8 @@ subroutine psb_srgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_srgmres + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ Parameters diff --git a/krylov/psb_z_inner_krylov_mod.f90 b/krylov/psb_z_inner_krylov_mod.f90 new file mode 100644 index 00000000..e6dc7029 --- /dev/null +++ b/krylov/psb_z_inner_krylov_mod.f90 @@ -0,0 +1,196 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.2 +!!$ (C) Copyright 2006/2007/2008 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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_krylov_mod.f90 +! Interfaces for Krylov subspace iterative methods. +! + +Module psb_z_inner_krylov_mod + + use psb_base_inner_krylov_mod + + interface psb_init_conv + module procedure psb_z_init_conv + end interface + + interface psb_check_conv + module procedure psb_z_check_conv + end interface + + + +contains + + + subroutine psb_z_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) + use psb_sparse_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: stopc, trace, itmax + type(psb_z_sparse_mat), intent(in) :: a + complex(psb_dpk_), intent(in) :: b(:) + real(psb_dpk_), intent(in) :: eps + type(psb_desc_type), intent(in) :: desc_a + type(psb_itconv_type) :: stopdat + integer, intent(out) :: info + + integer :: ictxt, me, np, err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_init_conv' + call psb_erractionsave(err_act) + + + ictxt=psb_cd_get_context(desc_a) + + call psb_info(ictxt, me, np) + + stopdat%controls(:) = 0 + stopdat%values(:) = 0.0d0 + + stopdat%controls(psb_ik_stopc_) = stopc + stopdat%controls(psb_ik_trace_) = trace + stopdat%controls(psb_ik_itmax_) = itmax + + select case(stopdat%controls(psb_ik_stopc_)) + case (1) + stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info) + if (info == psb_success_) stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info) + + case (2) + stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info) + + case default + info=psb_err_invalid_istop_ + call psb_errpush(info,name,i_err=(/stopc,0,0,0,0/)) + goto 9999 + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err="Init conv check data") + goto 9999 + end if + + stopdat%values(psb_ik_eps_) = eps + stopdat%values(psb_ik_errnum_) = dzero + stopdat%values(psb_ik_errden_) = done + + if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))& + & call log_header(methdname) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + + end subroutine psb_z_init_conv + + + function psb_z_check_conv(methdname,it,x,r,desc_a,stopdat,info) + use psb_sparse_mod + implicit none + character(len=*), intent(in) :: methdname + integer, intent(in) :: it + complex(psb_dpk_), intent(in) :: x(:), r(:) + type(psb_desc_type), intent(in) :: desc_a + type(psb_itconv_type) :: stopdat + logical :: psb_z_check_conv + integer, intent(out) :: info + + integer :: ictxt, me, np, err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_check_conv' + call psb_erractionsave(err_act) + + ictxt = psb_cd_get_context(desc_a) + call psb_info(ictxt,me,np) + psb_z_check_conv = .false. + + select case(stopdat%controls(psb_ik_stopc_)) + case(1) + stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info) + if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_) + stopdat%values(psb_ik_errden_) = & + & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)+stopdat%values(psb_ik_bni_)) + case(2) + stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info) + stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_) + stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_) + + case default + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err="Control data in stopdat messed up!") + goto 9999 + end select + if (info /= psb_success_) then + info=psb_err_from_subroutine_non_ + call psb_errpush(info,name) + goto 9999 + end if + + if (stopdat%values(psb_ik_errden_) == dzero) then + psb_z_check_conv = (stopdat%values(psb_ik_errnum_) <= stopdat%values(psb_ik_eps_)) + else + psb_z_check_conv = & + & (stopdat%values(psb_ik_errnum_) <= stopdat%values(psb_ik_eps_)*stopdat%values(psb_ik_errden_)) + end if + + psb_z_check_conv = (psb_z_check_conv.or.(stopdat%controls(psb_ik_itmax_) <= it)) + + if ( (stopdat%controls(psb_ik_trace_) > 0).and.& + & ((mod(it,stopdat%controls(psb_ik_trace_)) == 0).or.psb_z_check_conv)) then + call log_conv(methdname,me,it,1,stopdat%values(psb_ik_errnum_),& + & stopdat%values(psb_ik_errden_),stopdat%values(psb_ik_eps_)) + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + + end function psb_z_check_conv + + +end module psb_z_inner_krylov_mod diff --git a/krylov/psb_zbicg.f90 b/krylov/psb_zbicg.f90 index 8b318aaa..f62cdb98 100644 --- a/krylov/psb_zbicg.f90 +++ b/krylov/psb_zbicg.f90 @@ -96,7 +96,8 @@ subroutine psb_zbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_zbicg + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ parameters diff --git a/krylov/psb_zcg.F90 b/krylov/psb_zcg.F90 index 8a2593bc..260bfd82 100644 --- a/krylov/psb_zcg.F90 +++ b/krylov/psb_zcg.F90 @@ -98,7 +98,8 @@ subroutine psb_zcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_zcg + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ Parameters diff --git a/krylov/psb_zcgs.f90 b/krylov/psb_zcgs.f90 index 7cd6b7a2..2f8385ae 100644 --- a/krylov/psb_zcgs.f90 +++ b/krylov/psb_zcgs.f90 @@ -95,7 +95,8 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_zcgs + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ parameters diff --git a/krylov/psb_zcgstab.f90 b/krylov/psb_zcgstab.f90 index 0ca61e6b..ddf8ee55 100644 --- a/krylov/psb_zcgstab.f90 +++ b/krylov/psb_zcgstab.f90 @@ -96,7 +96,8 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_zcgstab + use psb_inner_krylov_mod + use psb_krylov_mod Implicit None !!$ parameters Type(psb_z_sparse_mat), Intent(in) :: a diff --git a/krylov/psb_zcgstabl.f90 b/krylov/psb_zcgstabl.f90 index ef1dfff7..287e3cd9 100644 --- a/krylov/psb_zcgstabl.f90 +++ b/krylov/psb_zcgstabl.f90 @@ -106,7 +106,8 @@ Subroutine psb_zcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_zcgstabl + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ parameters diff --git a/krylov/psb_zkrylov.f90 b/krylov/psb_zkrylov.f90 new file mode 100644 index 00000000..1428e98b --- /dev/null +++ b/krylov/psb_zkrylov.f90 @@ -0,0 +1,245 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.2 +!!$ (C) Copyright 2006/2007/2008 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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_krylov_mod.f90 +! Interfaces for Krylov subspace iterative methods. + ! + ! Subroutine: psb_zkrylov + ! + ! Front-end for the Krylov subspace iterations, complexversion + ! + ! Arguments: + ! + ! methd - character The specific method; can take the values: + ! CG + ! CGS + ! BICG + ! BICGSTAB + ! BICGSTABL + ! RGMRES + ! + ! a - type(psb_z_sparse_mat) 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 + ! irst - integer(optional) Input: restart parameter for RGMRES and + ! BICGSTAB(L) methods + ! istop - integer(optional) Input: stopping criterion, or how + ! to estimate the error. + ! 1: err = |r|/|b| + ! 2: err = |r|/(|a||x|+|b|) + ! where r is the (preconditioned, recursive + ! estimate of) residual + ! +Subroutine psb_zkrylov(method,a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) + use psb_sparse_mod + use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type + use psb_krylov_mod, psb_protect_name => psb_zkrylov + + character(len=*) :: method + Type(psb_z_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_zprec_type), intent(in) :: prec + complex(psb_dpk_), Intent(in) :: b(:) + complex(psb_dpk_), Intent(inout) :: x(:) + Real(psb_dpk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err + + interface + subroutine psb_zcg(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_zprec_type + type(psb_z_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + complex(psb_dpk_), intent(in) :: b(:) + complex(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(in) :: eps + class(psb_zprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_dpk_), optional, intent(out) :: err + end subroutine psb_zcg + subroutine psb_zbicg(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_zprec_type + type(psb_z_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + complex(psb_dpk_), intent(in) :: b(:) + complex(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(in) :: eps + class(psb_zprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_dpk_), optional, intent(out) :: err + end subroutine psb_zbicg + subroutine psb_zcgstab(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_zprec_type + type(psb_z_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + complex(psb_dpk_), intent(in) :: b(:) + complex(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(in) :: eps + class(psb_zprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_dpk_), optional, intent(out) :: err + end subroutine psb_zcgstab + Subroutine psb_zcgstabl(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_zprec_type + Type(psb_z_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_zprec_type), intent(in) :: prec + complex(psb_dpk_), Intent(in) :: b(:) + complex(psb_dpk_), Intent(inout) :: x(:) + Real(psb_dpk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err + end subroutine psb_zcgstabl + Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_zprec_type + Type(psb_z_sparse_mat), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_zprec_type), intent(in) :: prec + complex(psb_dpk_), Intent(in) :: b(:) + complex(psb_dpk_), Intent(inout) :: x(:) + Real(psb_dpk_), Intent(in) :: eps + integer, intent(out) :: info + Integer, Optional, Intent(in) :: itmax, itrace, irst,istop + Integer, Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err + end subroutine psb_zrgmres + subroutine psb_zcgs(a,prec,b,x,eps,& + & desc_a,info,itmax,iter,err,itrace,istop) + use psb_sparse_mod, only : psb_desc_type, psb_z_sparse_mat, psb_dpk_ + use psb_prec_mod, only : psb_zprec_type + type(psb_z_sparse_mat), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + complex(psb_dpk_), intent(in) :: b(:) + complex(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(in) :: eps + class(psb_zprec_type), intent(in) :: prec + integer, intent(out) :: info + integer, optional, intent(in) :: itmax, itrace,istop + integer, optional, intent(out) :: iter + real(psb_dpk_), optional, intent(out) :: err + end subroutine psb_zcgs + end interface + + integer :: ictxt,me,np,err_act + character(len=20) :: name + + info = psb_success_ + name = 'psb_krylov' + call psb_erractionsave(err_act) + + + ictxt=psb_cd_get_context(desc_a) + + call psb_info(ictxt, me, np) + + + select case(psb_toupper(method)) + case('CG') + call psb_zcg(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('CGS') + call psb_zcgs(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('BICG') + call psb_zbicg(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + case('BICGSTAB') + call psb_zcgstab(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,istop) + case('RGMRES') + call psb_zrgmres(a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,irst,istop) + case('BICGSTABL') + call psb_zcgstabl(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,irst,istop) + case default + if (me == 0) write(0,*) trim(name),': Warning: Unknown method ',method,& + & ', defaulting to BiCGSTAB' + call psb_zcgstab(a,prec,b,x,eps,desc_a,info,& + &itmax,iter,err,itrace,istop) + end select + + if(info /= psb_success_) then + 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 == psb_act_abort_) then + call psb_error(ictxt) + return + end if + +end subroutine psb_zkrylov + diff --git a/krylov/psb_zrgmres.f90 b/krylov/psb_zrgmres.f90 index c9a1e270..98d94a6d 100644 --- a/krylov/psb_zrgmres.f90 +++ b/krylov/psb_zrgmres.f90 @@ -108,7 +108,8 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop) use psb_sparse_mod use psb_prec_mod - use psb_krylov_mod, psb_protect_name => psb_zrgmres + use psb_inner_krylov_mod + use psb_krylov_mod implicit none !!$ Parameters