diff --git a/krylov/Makefile b/krylov/Makefile index b8579593..b2bdc2d0 100644 --- a/krylov/Makefile +++ b/krylov/Makefile @@ -35,8 +35,6 @@ symlink: (/bin/ln -fs $(PSBKRYLDIR)/*.f90 . ) (/bin/ln -fs $(PSBKRYLDIR)/*.F90 . ) -mld_krylov_mod.o: mld_krylov_mod.f90 $(LIBDIR)/mld_prec_mod$(.mod) - $(F90) $(F90COPT) $(FINCLUDES) -c $< psb_krylov_mod.o: psb_krylov_mod.f90 $(LIBDIR)/mld_prec_mod$(.mod) $(F90) $(F90COPT) $(FINCLUDES) -c $< diff --git a/krylov/mld_krylov_mod.f90 b/krylov/mld_krylov_mod.f90 deleted file mode 100644 index 1d2d6582..00000000 --- a/krylov/mld_krylov_mod.f90 +++ /dev/null @@ -1,487 +0,0 @@ -!!$ -!!$ -!!$ MLD2P4 -!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS v.2.0) -!!$ -!!$ (C) Copyright 2006 Alfredo Buttari University of Rome Tor Vergata -!!$ Pasqua D'Ambra ICAR-CNR, Naples -!!$ Daniela di Serafino Second University of Naples -!!$ Salvatore Filippone 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 MLD2P4 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 MLD2P4 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: mld_krylov_mod.f90. -! -! Module: mld_krylov_mod. -! Contains: mld_dkrylov, mld_zkrylov. -! -! This module defines the interfaces to the real and complex versions of the -! routines implementing Krylov solvers: -! - mld_krylov : driver providing a general interface to the Krylov -! solvers available in PSBLAS; -! - psb_cg : Conjugate Gradient (CG) solver (real version only); -! - psb_bicg : BiConjugate Gradient (BiCG) solver (real version only); -! - psb_bicgstab : BiConjugate Gradient Stabilized (BiCGStab) solver; -! - psb_bicgstabl : BiCGStab(l) solver (real version only); -! - psb_rgmres : Generalized Minimal Residual (GMRES) solver with -! Restarting; -! - psb_cgs : Conjugate Gradient Squared (CGS) solver. -! -! This is a partial duplication of material already present in the Krylov -! part of the PSBLAS subtree. This is not bad since -! 1. there is no actual computation here, these are just wrappers around -! the actual methods, so very little code is necessary; -! 2. these wrappers enable issuing error messages with a subroutine name -! coherent with what the user sees; -! 3. at the very least we would need a layer renaming the PSBLAS krylov -! interfaces, so we are doing approximately the same amount of (re)coding. -! -Module mld_krylov_mod - - interface mld_krylov - module procedure mld_dkrylov, mld_zkrylov - end interface - - interface psb_cg - subroutine psb_dcg(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_base_mod - use psb_prec_mod - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(kind(1.d0)), intent(in) :: b(:) - real(kind(1.d0)), intent(inout) :: x(:) - real(kind(1.d0)), intent(in) :: eps - type(psb_dprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(kind(1.d0)), optional, intent(out) :: err - end subroutine psb_dcg - end interface - - interface psb_bicg - subroutine psb_dbicg(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_base_mod - use psb_prec_mod - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(kind(1.d0)), intent(in) :: b(:) - real(kind(1.d0)), intent(inout) :: x(:) - real(kind(1.d0)), intent(in) :: eps - type(psb_dprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(kind(1.d0)), optional, intent(out) :: err - end subroutine psb_dbicg - end interface - - interface psb_bicgstab - subroutine psb_dcgstab(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_base_mod - use psb_prec_mod - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(kind(1.d0)), intent(in) :: b(:) - real(kind(1.d0)), intent(inout) :: x(:) - real(kind(1.d0)), intent(in) :: eps - type(psb_dprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(kind(1.d0)), optional, intent(out) :: err - end subroutine psb_dcgstab - subroutine psb_zcgstab(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_base_mod - use psb_prec_mod - type(psb_zspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - complex(kind(1.d0)), intent(in) :: b(:) - complex(kind(1.d0)), intent(inout) :: x(:) - real(kind(1.d0)), intent(in) :: eps - type(psb_zprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(kind(1.d0)), optional, intent(out) :: err - end subroutine psb_zcgstab - end interface - - interface psb_bicgstabl - Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err, itrace,irst,istop) - use psb_base_mod - use psb_prec_mod - Type(psb_dspmat_type), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - type(psb_dprec_type), intent(in) :: prec - Real(Kind(1.d0)), Intent(in) :: b(:) - Real(Kind(1.d0)), Intent(inout) :: x(:) - Real(Kind(1.d0)), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(Kind(1.d0)), Optional, Intent(out) :: err - end subroutine psb_dcgstabl - end interface - - interface psb_rgmres - Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - use psb_base_mod - use psb_prec_mod - Type(psb_dspmat_type), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - type(psb_dprec_type), intent(in) :: prec - Real(Kind(1.d0)), Intent(in) :: b(:) - Real(Kind(1.d0)), Intent(inout) :: x(:) - Real(Kind(1.d0)), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(Kind(1.d0)), Optional, Intent(out) :: err - end subroutine psb_drgmres - Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - use psb_base_mod - use psb_prec_mod - Type(psb_zspmat_type), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - type(psb_zprec_type), intent(in) :: prec - complex(Kind(1.d0)), Intent(in) :: b(:) - complex(Kind(1.d0)), Intent(inout) :: x(:) - Real(Kind(1.d0)), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(Kind(1.d0)), Optional, Intent(out) :: err - end subroutine psb_zrgmres - end interface - - interface psb_cgs - subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - use psb_base_mod - use psb_prec_mod - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_dprec_type), intent(in) :: prec - real(kind(1.d0)), intent(in) :: b(:) - real(kind(1.d0)), intent(inout) :: x(:) - real(kind(1.d0)), intent(in) :: eps - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(kind(1.d0)), optional, intent(out) :: err - end subroutine psb_dcgs - subroutine psb_zcgs(a,prec,b,x,eps,& - & desc_a,info,itmax,iter,err,itrace,istop) - use psb_base_mod - use psb_prec_mod - type(psb_zspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - complex(kind(1.d0)), intent(in) :: b(:) - complex(kind(1.d0)), intent(inout) :: x(:) - real(kind(1.d0)), intent(in) :: eps - type(psb_zprec_type), intent(in) :: prec - integer, intent(out) :: info - integer, optional, intent(in) :: itmax, itrace,istop - integer, optional, intent(out) :: iter - real(kind(1.d0)), optional, intent(out) :: err - end subroutine psb_zcgs - end interface - -contains - - ! - ! Subroutine: mld_dkrylov/mld_zkrylov. - ! Version: real/complex. - ! Note: internal subroutine of module mld_krylov_mod. - ! - ! This routine solves a linear system - ! - ! Ax = b - ! - ! by using one of following preconditioned Krylov solvers: CG, BiCG, BiCGStab, - ! BiCGStab(l), Restarted GMRES, CGS. The preconditioning is applied as right - ! preconditioning. - ! - ! The user can choose among two stopping criteria: - ! - ! - one based on the normwise backward error in the infinity norm, i.e. - ! - ! ||r_i||_inf < eps * (||A||_inf * ||x_i||_inf + ||b||_inf); - ! - ! - the other based on the relative residual reduction in the 2-norm, i.e. - ! - ! ||r_i||_2 < eps * ||b||_2. - ! - ! In the above formulae, x_i is the approximate solution and r_i = b - A*x_i - ! is the corresponding residual at the i-th iteration. - ! - ! This routine is a driver that provides a general interface to the Krylov - ! solvers implemented in PSBLAS. These solvers can use all the preconditioners - ! available in MLD2P4. - ! - ! - ! Arguments: - ! method - character(len=*), input. - ! The Krylov method to be applied. The following values are - ! are allowed in the real case: 'cg', 'bicg', 'bicgstab', - ! 'bicgstabl', 'rgmres', 'cgs' (and all the upper/lower case - ! variants). The same values except 'cg', 'bicg' and 'bicgstabl' - ! are allowed in the complex case. - ! a - type(psb_dspmat_type)/type(psb_zspmat_type), input. - ! The sparse matrix structure containing the local part of the - ! matrix A. - ! prec - type(mld_dprec_type)/type(mld_zprec_type), input. - ! The preconditioner data structure containing the local part - ! of the preconditioner to be applied. - ! b - real(kind(1.d0))/complex(kind(1.d0)), dimension(:), input. - ! The local part of the right-hand side. - ! x - real(kind(1.d0))/complex(kind(1.d0)), dimension(:), input/output. - ! The local part of the approximate solution. In input it - ! contains an initial guess of the solution, while, in output, - ! the approximation computed by the selected Krylov solver. - ! eps - real(kind(1.d0)), input. - ! The tolerance used in the stopping criterion. - ! desc_a - type(psb_desc_type), input. - ! The communication descriptor associated to a. - ! info - integer, output. - ! Error code. - ! itmax - integer, optional, input. - ! The maximum number of iterations to be performed. If itmax - ! is not present, a maximum number of 1000 iterations is - ! assumed. - ! iter - integer, optional, output. - ! The number of iterations actually performed. - ! err - real(kind(1.d0)), optional, output. - ! An estimate of the error in the computed solution, i.e. the - ! relative backward error ||r|| / (||A||*||x||+||b||), where - ! || || is the infinity norm, or the relative residual - ! ||r||/||b||, where || || is the 2-norm, according to the - ! selected stopping criterion (see istop below). - ! itrace - integer, optional, input. - ! If itrace>0 then information about convergence is printed - ! every itrace iterations; otherwise it is ignored. If itrace - ! is not present, itrace=0 is assumed. - ! irst - integer, optional, input. - ! The restart parameter used in the GMRES and BiCGSTAB(l) - ! solvers. It is ignored if other solvers are chosen. If irst - ! is not present, irst=10 is assumed for Restarded GMRES and - ! irst=l=1 for BiCGStab(l) (note that BiCGStab(1)=BiCGStab). - ! istop - integer, optional, input. - ! The normwise backward error stopping criterion is used if - ! istop=1; the relative residual one if istop=2. If istop - ! is not present, istop=1 is assumed. - ! - - Subroutine mld_dkrylov(method,a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - - use psb_base_mod - use psb_prec_mod - - ! Arguments - character(len=*), intent(in) :: method - Type(psb_dspmat_type), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - type(psb_dprec_type), intent(in) :: prec - Real(Kind(1.d0)), Intent(in) :: b(:) - Real(Kind(1.d0)), Intent(inout) :: x(:) - Real(Kind(1.d0)), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(Kind(1.d0)), Optional, Intent(out) :: err - - ! Local variables - integer :: ictxt,me,np,err_act - character(len=20) :: name - - info = 0 - name = 'mld_krylov' - call psb_erractionsave(err_act) - - ictxt=psb_cd_get_context(desc_a) - - call psb_info(ictxt, me, np) - - select case(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,*) 'Warning: Unknown method ',method,& - & ' in ',name,' defaulting to BiCGSTAB' - call psb_bicgstab(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - end select - - if(info/=0) 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.eq.psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end subroutine mld_dkrylov - - Subroutine mld_zkrylov(method,a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,irst,istop) - - use psb_base_mod - use psb_prec_mod - - ! Arguments - character(len=*), intent(in) :: method - Type(psb_zspmat_type), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - type(psb_zprec_type), intent(in) :: prec - complex(Kind(1.d0)), Intent(in) :: b(:) - complex(Kind(1.d0)), Intent(inout) :: x(:) - Real(Kind(1.d0)), Intent(in) :: eps - integer, intent(out) :: info - Integer, Optional, Intent(in) :: itmax, itrace, irst,istop - Integer, Optional, Intent(out) :: iter - Real(Kind(1.d0)), Optional, Intent(out) :: err - - ! Local variables - integer :: ictxt,me,np,err_act - character(len=20) :: name - - info = 0 - name = 'mld_krylov' - call psb_erractionsave(err_act) - - - ictxt=psb_cd_get_context(desc_a) - - call psb_info(ictxt, me, np) - - - select case(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,*) 'Warning: Unknown method ',method,& - & ' in ',name,' defaulting to BiCGSTAB' - call psb_bicgstab(a,prec,b,x,eps,desc_a,info,& - &itmax,iter,err,itrace,istop) - end select - - if(info/=0) 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.eq.psb_act_abort_) then - call psb_error(ictxt) - return - end if - - end subroutine mld_zkrylov - - subroutine log_conv(methdname,me,itx,itrace,errnum,errden,eps) - character(len=*), intent(in) :: methdname - integer, intent(in) :: me, itx, itrace - real(kind(1.d0)), intent(in) :: errnum, errden, eps - character(len=*), parameter :: fmt='(a,i4,3(2x,es10.4))' - - if ((mod(itx,itrace)==0).and.(me == 0)) then - write(*,fmt) methdname//': ',itx,errnum,eps*errden - Endif - - end subroutine log_conv - - subroutine end_log(methdname,me,itx,errnum,errden,eps) - character(len=*), intent(in) :: methdname - integer, intent(in) :: me, itx - real(kind(1.d0)), intent(in) :: errnum, errden, eps - character(len=*), parameter :: fmt='(a,2x,es10.4,1x,a,1x,i4,1x,a)' - character(len=*), parameter :: fmt1='(a,3(2x,es10.4))' - - if (me==0) then - write(*,fmt) methdname//' failed to converge to ',eps,& - & ' in ',ITX,' iterations. ' - write(*,fmt1) 'Last iteration convergence indicators: ',& - & errnum,eps*errden,errnum/errden - end if - end subroutine end_log - -end module mld_krylov_mod - - -