Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_Krylov_def.hpp
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 //-----------------------------------------------------
00031 // Ifpack2::KRYLOV is based on the Krylov iterative
00032 // solvers in Belos.
00033 // written by Paul Tsuji.
00034 //-----------------------------------------------------
00035 
00036 #ifndef IFPACK2_KRYLOV_DEF_HPP
00037 #define IFPACK2_KRYLOV_DEF_HPP
00038 
00039 namespace Ifpack2 {
00040 
00041 //Definitions for the Krylov methods:
00042 //==============================================================================
00043 template <class MatrixType, class PrecType>
00044 Krylov<MatrixType,PrecType>::Krylov(const Teuchos::RCP<const Tpetra::RowMatrix<scalar_type,local_ordinal_type,global_ordinal_type,node_type> >& A) :
00045   A_(A),
00046   Comm_(A->getRowMap()->getComm()),
00047   // Default values
00048   IterationType_(1),
00049   Iterations_(5),
00050   ResidualTolerance_(0.001),
00051   BlockSize_(1),
00052   ZeroStartingSolution_(true),
00053   PreconditionerType_(1),
00054   // General
00055   Condest_(-1.0),
00056   IsInitialized_(false),
00057   IsComputed_(false),
00058   NumInitialize_(0),
00059   NumCompute_(0),
00060   NumApply_(0),
00061   InitializeTime_(0.0),
00062   ComputeTime_(0.0),
00063   ApplyTime_(0.0),
00064   Time_("Ifpack2::Krylov"),
00065   NumMyRows_(-1),
00066   NumGlobalNonzeros_(0)
00067 {
00068   TEUCHOS_TEST_FOR_EXCEPTION(A_ == Teuchos::null, std::runtime_error,
00069       Teuchos::typeName(*this) << "::Krylov(): input matrix reference was null.");
00070 }
00071 
00072 //==========================================================================
00073 template <class MatrixType, class PrecType>
00074 Krylov<MatrixType,PrecType>::~Krylov() {
00075 }
00076 
00077 //==========================================================================
00078 template <class MatrixType, class PrecType>
00079 void Krylov<MatrixType,PrecType>::setParameters(const Teuchos::ParameterList& params) {
00080   using Teuchos::as;
00081   using Teuchos::Exceptions::InvalidParameterName;
00082   using Teuchos::Exceptions::InvalidParameterType;
00083   typedef Teuchos::ScalarTraits<magnitude_type> STM;
00084   // Read in parameters
00085   Ifpack2::getParameter(params, "krylov: number of iterations",Iterations_);
00086   Ifpack2::getParameter(params, "krylov: residual tolerance",ResidualTolerance_);
00087   Ifpack2::getParameter(params, "krylov: block size",BlockSize_);
00088   Ifpack2::getParameter(params, "krylov: zero starting solution",ZeroStartingSolution_);
00089   Ifpack2::getParameter(params, "krylov: preconditioner type",PreconditionerType_);
00090   params_=params;
00091 }
00092 
00093 //==========================================================================
00094 template <class MatrixType, class PrecType>
00095 const Teuchos::RCP<const Teuchos::Comm<int> > &
00096 Krylov<MatrixType,PrecType>::getComm() const{
00097   return(Comm_);
00098 }
00099 
00100 //==========================================================================
00101 template <class MatrixType, class PrecType>
00102 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
00103 Krylov<MatrixType,PrecType>::getMatrix() const {
00104   return(A_);
00105 }
00106 
00107 //==========================================================================
00108 template <class MatrixType, class PrecType>
00109 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >&
00110 Krylov<MatrixType,PrecType>::getDomainMap() const
00111 {
00112   return A_->getDomainMap();
00113 }
00114 
00115 //==========================================================================
00116 template <class MatrixType, class PrecType>
00117 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >&
00118 Krylov<MatrixType,PrecType>::getRangeMap() const
00119 {
00120   return A_->getRangeMap();
00121 }
00122 
00123 //==============================================================================
00124 template <class MatrixType, class PrecType>
00125 bool Krylov<MatrixType,PrecType>::hasTransposeApply() const {
00126   return true;
00127 }
00128 
00129 //==========================================================================
00130 template <class MatrixType, class PrecType>
00131 int Krylov<MatrixType,PrecType>::getNumInitialize() const {
00132   return(NumInitialize_);
00133 }
00134 
00135 //==========================================================================
00136 template <class MatrixType, class PrecType>
00137 int Krylov<MatrixType,PrecType>::getNumCompute() const {
00138   return(NumCompute_);
00139 }
00140 
00141 //==========================================================================
00142 template <class MatrixType, class PrecType>
00143 int Krylov<MatrixType,PrecType>::getNumApply() const {
00144   return(NumApply_);
00145 }
00146 
00147 //==========================================================================
00148 template <class MatrixType, class PrecType>
00149 double Krylov<MatrixType,PrecType>::getInitializeTime() const {
00150   return(InitializeTime_);
00151 }
00152 
00153 //==========================================================================
00154 template <class MatrixType, class PrecType>
00155 double Krylov<MatrixType,PrecType>::getComputeTime() const {
00156   return(ComputeTime_);
00157 }
00158 
00159 //==========================================================================
00160 template <class MatrixType, class PrecType>
00161 double Krylov<MatrixType,PrecType>::getApplyTime() const {
00162   return(ApplyTime_);
00163 }
00164 
00165 //=============================================================================
00166 template <class MatrixType, class PrecType>
00167 typename Krylov<MatrixType,PrecType>::magnitude_type
00168 Krylov<MatrixType,PrecType>::
00169 computeCondEst (CondestType CT,
00170                 local_ordinal_type MaxIters,
00171                 magnitude_type Tol,
00172                 const Teuchos::Ptr<const Tpetra::RowMatrix<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > &matrix) {
00173   if (!isComputed()) { // cannot compute right now
00174     return(-1.0);
00175   }
00176   // NOTE: this is computing the *local* condest
00177   if (Condest_ == -1.0) {
00178     Condest_ = Ifpack2::Condest(*this, CT, MaxIters, Tol, matrix);
00179   }
00180   return(Condest_);
00181 }
00182 
00183 //==========================================================================
00184 template <class MatrixType, class PrecType>
00185 void Krylov<MatrixType,PrecType>::initialize() {
00186   // clear any previous allocation
00187   IsInitialized_ = false;
00188   IsComputed_ = false;
00189   Time_.start(true);
00190   // check only in serial
00191   TEUCHOS_TEST_FOR_EXCEPTION(Comm_->getSize() == 1 && A_->getNodeNumRows() != A_->getNodeNumCols(), std::runtime_error, "Ifpack2::Krylov::Initialize ERROR, matrix must be square");
00192   NumMyRows_ = A_->getNodeNumRows();
00193   // Belos parameter list
00194   belosList_ = Teuchos::rcp( new Teuchos::ParameterList("GMRES") );
00195   belosList_->set("Maximum Iterations",Iterations_);
00196   belosList_->set("Convergence Tolerance",ResidualTolerance_);
00197   //belosList_->set("Verbosity", Belos::Errors + Belos::Warnings + Belos::TimingDetails + Belos::StatusTestDetails);
00198   // nothing else to do here
00199   ifpack2_prec_=Teuchos::rcp( new PrecType(A_) );
00200   ifpack2_prec_->initialize();
00201   ifpack2_prec_->setParameters(params_);
00202   IsInitialized_ = true;
00203   ++NumInitialize_;
00204   Time_.stop();
00205   InitializeTime_ += Time_.totalElapsedTime();
00206 }
00207 
00208 //==========================================================================
00209 template <class MatrixType, class PrecType>
00210 void Krylov<MatrixType,PrecType>::compute() {
00211   using Teuchos::as;
00212   //--------------------------------------------------------------------------
00213   // Ifpack2::Krylov
00214   // Computes necessary preconditioner info
00215   //--------------------------------------------------------------------------
00216   if (!isInitialized()) {
00217     initialize();
00218   }
00219   Time_.start(true);
00220   ifpack2_prec_->compute();
00221   IsComputed_ = true;
00222   ++NumCompute_;
00223   Time_.stop();
00224   ComputeTime_ += Time_.totalElapsedTime();
00225 }
00226 
00227 //==========================================================================
00228 template <class MatrixType, class PrecType>
00229 void Krylov<MatrixType,PrecType>::apply(const Tpetra::MultiVector<typename MatrixType::scalar_type,
00230             typename MatrixType::local_ordinal_type,
00231             typename MatrixType::global_ordinal_type,
00232             typename MatrixType::node_type>& X,
00233             Tpetra::MultiVector<typename MatrixType::scalar_type,
00234             typename MatrixType::local_ordinal_type,
00235             typename MatrixType::global_ordinal_type,
00236             typename MatrixType::node_type>& Y,
00237             Teuchos::ETransp mode,
00238             typename MatrixType::scalar_type alpha,
00239             typename MatrixType::scalar_type beta) const
00240 {
00241   TEUCHOS_TEST_FOR_EXCEPTION(!isComputed(), std::runtime_error,
00242     "Ifpack2::Krylov::apply() ERROR, compute() hasn't been called yet.");
00243 
00244   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00245     "Ifpack2::Krylov::apply() ERROR, X.getNumVectors() != Y.getNumVectors().");
00246 
00247   Time_.start(true);
00248 
00249   // If X and Y are pointing to the same memory location,
00250   // we need to create an auxiliary vector, Xcopy
00251   Teuchos::RCP< const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > Xcopy;
00252   if (X.getLocalMV().getValues() == Y.getLocalMV().getValues()) {
00253     Xcopy = Teuchos::rcp( new Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>(X) );
00254   }
00255   else {
00256     Xcopy = Teuchos::rcp( &X, false );
00257   }
00258 
00259   const Teuchos::RCP< Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > Zeros = 
00260     Teuchos::rcp( new Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> (A_->getRowMap(),X.getNumVectors()) );
00261   Zeros->putScalar((scalar_type) 0.0);
00262 
00263   // Belos Linear Problem
00264   Teuchos::RCP< Belos::LinearProblem<belos_scalar_type,
00265     Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>,
00266     Tpetra::Operator<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > > belosProblem = 
00267     Teuchos::rcp( new Belos::LinearProblem<belos_scalar_type,
00268       Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>,
00269       Tpetra::Operator<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > (A_, Zeros, Xcopy) );
00270   belosProblem->setLeftPrec(ifpack2_prec_);
00271   belosProblem->setProblem(Zeros,Xcopy);
00272   // Belos Solver Manager
00273   Teuchos::RCP< Belos::SolverManager<belos_scalar_type,
00274     Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>,
00275     Tpetra::Operator<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > > belosSolver = 
00276     Teuchos::rcp( new Belos::BlockGmresSolMgr<belos_scalar_type,
00277       Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>,
00278       Tpetra::Operator<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > (belosProblem, belosList_) );
00279   // iterative solve
00280   belosSolver->solve();
00281   Y=*Zeros;
00282   ++NumApply_;
00283   Time_.stop();
00284   ApplyTime_ += Time_.totalElapsedTime();
00285 }
00286 
00287 //=============================================================================
00288 template <class MatrixType, class PrecType>
00289 std::string Krylov<MatrixType,PrecType>::description() const {
00290   std::ostringstream oss;
00291   oss << Teuchos::Describable::description();
00292   if (isInitialized()) {
00293     if (isComputed()) {
00294       oss << "{status = initialized, computed";
00295     }
00296     else {
00297       oss << "{status = initialized, not computed";
00298     }
00299   }
00300   else {
00301     oss << "{status = not initialized, not computed";
00302   }
00303   oss << ", global rows = " << A_->getGlobalNumRows()
00304       << ", global cols = " << A_->getGlobalNumCols()
00305       << "}";
00306   return oss.str();
00307 }
00308 
00309 //=============================================================================
00310 template <class MatrixType, class PrecType>
00311 void Krylov<MatrixType,PrecType>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00312   using std::endl;
00313   using std::setw;
00314   using Teuchos::VERB_DEFAULT;
00315   using Teuchos::VERB_NONE;
00316   using Teuchos::VERB_LOW;
00317   using Teuchos::VERB_MEDIUM;
00318   using Teuchos::VERB_HIGH;
00319   using Teuchos::VERB_EXTREME;
00320   Teuchos::EVerbosityLevel vl = verbLevel;
00321   if (vl == VERB_DEFAULT) vl = VERB_LOW;
00322   const int myImageID = Comm_->getRank();
00323   Teuchos::OSTab tab(out);
00324   //    none: print nothing
00325   //     low: print O(1) info from node 0
00326   //  medium:
00327   //    high:
00328   // extreme:
00329   if (vl != VERB_NONE && myImageID == 0) {
00330   }
00331 }
00332 
00333 }//namespace Ifpack2
00334 
00335 #endif /* IFPACK2_KRYLOV_DEF_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends