BelosGmresPolySolMgr.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                 Copyright (2004) 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 #ifndef BELOS_GMRES_POLY_SOLMGR_HPP
00030 #define BELOS_GMRES_POLY_SOLMGR_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 
00039 #include "BelosLinearProblem.hpp"
00040 #include "BelosSolverManager.hpp"
00041 
00042 #include "BelosGmresIteration.hpp"
00043 #include "BelosBlockGmresIter.hpp"
00044 #include "BelosDGKSOrthoManager.hpp"
00045 #include "BelosICGSOrthoManager.hpp"
00046 #include "BelosIMGSOrthoManager.hpp"
00047 #include "BelosStatusTestMaxIters.hpp"
00048 #include "BelosStatusTestGenResNorm.hpp"
00049 #include "BelosStatusTestImpResNorm.hpp"
00050 #include "BelosStatusTestCombo.hpp"
00051 #include "BelosStatusTestOutput.hpp"
00052 #include "BelosOutputManager.hpp"
00053 #include "Teuchos_BLAS.hpp"
00054 #include "Teuchos_LAPACK.hpp"
00055 #include "Teuchos_TimeMonitor.hpp"
00056 
00072 namespace Belos {
00073   
00075 
00076   
00083 class GmresPolySolMgrLinearProblemFailure : public BelosError {public:
00084   GmresPolySolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00085     {}};
00086   
00093 class GmresPolySolMgrOrthoFailure : public BelosError {public:
00094   GmresPolySolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00095     {}};
00096   
00097 template<class ScalarType, class MV, class OP>
00098 class GmresPolySolMgr : public SolverManager<ScalarType,MV,OP> {
00099     
00100 private:
00101   typedef MultiVecTraits<ScalarType,MV> MVT;
00102   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00103   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00104   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00105   typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00106     
00107 public:
00108     
00110 
00111    
00117   GmresPolySolMgr();
00118  
00131   GmresPolySolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00132     const Teuchos::RCP<Teuchos::ParameterList> &pl );
00133     
00135   virtual ~GmresPolySolMgr() {};
00137     
00139 
00140     
00143   const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00144     return *problem_;
00145   }
00146 
00149   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00150 
00153   Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00154  
00160   Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00161     return tuple(timerSolve_);
00162   }
00163   
00165   int getNumIters() const {
00166     return numIters_;
00167   }
00168  
00172   bool isLOADetected() const { return loaDetected_; }
00173  
00175     
00177 
00178     
00179   void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; isSTSet_ = false; }
00180     
00181   void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00182     
00184     
00186 
00187     
00205   ReturnType solve();
00206     
00208     
00211     
00213   std::string description() const;
00214     
00216     
00217 private:
00218 
00219   // Method to convert std::string to enumerated type for residual.
00220   Belos::ScaleType convertStringToScaleType( std::string& scaleType ) {
00221     if (scaleType == "Norm of Initial Residual") {
00222       return Belos::NormOfInitRes;
00223     } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00224       return Belos::NormOfPrecInitRes;
00225     } else if (scaleType == "Norm of RHS") {
00226       return Belos::NormOfRHS;
00227     } else if (scaleType == "None") {
00228       return Belos::None;
00229     } else 
00230       TEST_FOR_EXCEPTION( true ,std::logic_error,
00231         "Belos::GmresPolySolMgr(): Invalid residual scaling type.");
00232   }
00233   
00234   // Method for checking current status test against defined linear problem.
00235   bool checkStatusTest();
00236 
00237   // Linear problem.
00238   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00239     
00240   // Output manager.
00241   Teuchos::RCP<OutputManager<ScalarType> > printer_;
00242   Teuchos::RCP<std::ostream> outputStream_;
00243 
00244   // Status test.
00245   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00246   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00247   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00248   Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
00249   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00250 
00251   // Orthogonalization manager.
00252   Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00253     
00254   // Current parameter list.
00255   Teuchos::RCP<ParameterList> params_;
00256 
00257   // Default solver values.
00258   static const MagnitudeType convtol_default_;
00259   static const MagnitudeType orthoKappa_default_;
00260   static const int maxRestarts_default_;
00261   static const int maxIters_default_;
00262   static const bool showMaxResNormOnly_default_;
00263   static const int blockSize_default_;
00264   static const int numBlocks_default_;
00265   static const int verbosity_default_;
00266   static const int outputFreq_default_;
00267   static const std::string impResScale_default_; 
00268   static const std::string expResScale_default_; 
00269   static const std::string label_default_;
00270   static const std::string orthoType_default_;
00271   static const Teuchos::RCP<std::ostream> outputStream_default_;
00272 
00273   // Current solver values.
00274   MagnitudeType convtol_, orthoKappa_;
00275   int maxRestarts_, maxIters_, numIters_;
00276   int blockSize_, numBlocks_, verbosity_, outputFreq_;
00277   bool showMaxResNormOnly_;
00278   std::string orthoType_; 
00279   std::string impResScale_, expResScale_;
00280     
00281   // Timers.
00282   std::string label_;
00283   Teuchos::RCP<Teuchos::Time> timerSolve_;
00284 
00285   // Internal state variables.
00286   bool isPolyBuilt_;  
00287   bool isSet_, isSTSet_, expResTest_;
00288   bool loaDetected_;
00289 };
00290 
00291 
00292 // Default solver values.
00293 template<class ScalarType, class MV, class OP>
00294 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType GmresPolySolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00295 
00296 template<class ScalarType, class MV, class OP>
00297 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType GmresPolySolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00298 
00299 template<class ScalarType, class MV, class OP>
00300 const int GmresPolySolMgr<ScalarType,MV,OP>::maxRestarts_default_ = 20;
00301 
00302 template<class ScalarType, class MV, class OP>
00303 const int GmresPolySolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00304 
00305 template<class ScalarType, class MV, class OP>
00306 const bool GmresPolySolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00307 
00308 template<class ScalarType, class MV, class OP>
00309 const int GmresPolySolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00310 
00311 template<class ScalarType, class MV, class OP>
00312 const int GmresPolySolMgr<ScalarType,MV,OP>::numBlocks_default_ = 300;
00313 
00314 template<class ScalarType, class MV, class OP>
00315 const int GmresPolySolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00316 
00317 template<class ScalarType, class MV, class OP>
00318 const int GmresPolySolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00319 
00320 template<class ScalarType, class MV, class OP>
00321 const std::string GmresPolySolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of Preconditioned Initial Residual";
00322 
00323 template<class ScalarType, class MV, class OP>
00324 const std::string GmresPolySolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of Initial Residual";
00325 
00326 template<class ScalarType, class MV, class OP>
00327 const std::string GmresPolySolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00328 
00329 template<class ScalarType, class MV, class OP>
00330 const std::string GmresPolySolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00331 
00332 template<class ScalarType, class MV, class OP>
00333 const Teuchos::RCP<std::ostream> GmresPolySolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00334 
00335 
00336 // Empty Constructor
00337 template<class ScalarType, class MV, class OP>
00338 GmresPolySolMgr<ScalarType,MV,OP>::GmresPolySolMgr() :
00339   outputStream_(outputStream_default_),
00340   convtol_(convtol_default_),
00341   orthoKappa_(orthoKappa_default_),
00342   maxRestarts_(maxRestarts_default_),
00343   maxIters_(maxIters_default_),
00344   blockSize_(blockSize_default_),
00345   numBlocks_(numBlocks_default_),
00346   verbosity_(verbosity_default_),
00347   outputFreq_(outputFreq_default_),
00348   showMaxResNormOnly_(showMaxResNormOnly_default_),
00349   orthoType_(orthoType_default_),
00350   impResScale_(impResScale_default_),
00351   expResScale_(expResScale_default_),
00352   label_(label_default_),
00353   isPolyBuilt_(false),
00354   isSet_(false),
00355   isSTSet_(false),
00356   expResTest_(false),
00357   loaDetected_(false)
00358 {}
00359 
00360 
00361 // Basic Constructor
00362 template<class ScalarType, class MV, class OP>
00363 GmresPolySolMgr<ScalarType,MV,OP>::GmresPolySolMgr( 
00364   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00365   const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00366   problem_(problem),
00367   outputStream_(outputStream_default_),
00368   convtol_(convtol_default_),
00369   orthoKappa_(orthoKappa_default_),
00370   maxRestarts_(maxRestarts_default_),
00371   maxIters_(maxIters_default_),
00372   blockSize_(blockSize_default_),
00373   numBlocks_(numBlocks_default_),
00374   verbosity_(verbosity_default_),
00375   outputFreq_(outputFreq_default_),
00376   showMaxResNormOnly_(showMaxResNormOnly_default_),
00377   orthoType_(orthoType_default_),
00378   impResScale_(impResScale_default_),
00379   expResScale_(expResScale_default_),
00380   label_(label_default_),
00381   isPolyBuilt_(false),
00382   isSet_(false),
00383   isSTSet_(false),
00384   expResTest_(false),
00385   loaDetected_(false)
00386 {
00387 
00388   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00389 
00390   // If the parameter list pointer is null, then set the current parameters to the default parameter list.
00391   if ( !is_null(pl) ) {
00392     setParameters( pl );  
00393   }
00394 
00395 }
00396 
00397 
00398 template<class ScalarType, class MV, class OP>
00399 Teuchos::RCP<const Teuchos::ParameterList>
00400 GmresPolySolMgr<ScalarType,MV,OP>::getValidParameters() const
00401 {
00402   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00403   if (is_null(validPL)) {
00404     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00405     pl->set("Convergence Tolerance", convtol_default_,
00406       "The relative residual tolerance that needs to be achieved by the\n"
00407       "iterative solver in order for the linear system to be declared converged." );
00408     pl->set("Maximum Restarts", maxRestarts_default_,
00409       "The maximum number of restarts allowed for each\n"
00410       "set of RHS solved.");
00411     pl->set(
00412       "Maximum Iterations", maxIters_default_,
00413       "The maximum number of block iterations allowed for each\n"
00414       "set of RHS solved.");
00415     pl->set("Num Blocks", numBlocks_default_,
00416       "The maximum number of blocks allowed in the Krylov subspace\n"
00417       "for each set of RHS solved.");
00418     pl->set("Block Size", blockSize_default_,
00419       "The number of vectors in each block.  This number times the\n"
00420       "number of blocks is the total Krylov subspace dimension.");
00421     pl->set("Verbosity", verbosity_default_,
00422       "What type(s) of solver information should be outputted\n"
00423       "to the output stream.");
00424     pl->set("Output Frequency", outputFreq_default_,
00425       "How often convergence information should be outputted\n"
00426       "to the output stream.");  
00427     pl->set("Output Stream", outputStream_default_,
00428       "A reference-counted pointer to the output stream where all\n"
00429       "solver output is sent.");
00430     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00431       "When convergence information is printed, only show the maximum\n"
00432       "relative residual norm when the block size is greater than one.");
00433     pl->set("Implicit Residual Scaling", impResScale_default_,
00434       "The type of scaling used in the implicit residual convergence test.");
00435     pl->set("Explicit Residual Scaling", expResScale_default_,
00436       "The type of scaling used in the explicit residual convergence test.");
00437     pl->set("Timer Label", label_default_,
00438       "The string to use as a prefix for the timer labels.");
00439     //  pl->set("Restart Timers", restartTimers_);
00440     pl->set("Orthogonalization", orthoType_default_,
00441       "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
00442     pl->set("Orthogonalization Constant",orthoKappa_default_,
00443       "The constant used by DGKS orthogonalization to determine\n"
00444       "whether another step of classical Gram-Schmidt is necessary.");
00445     validPL = pl;
00446   }
00447   return validPL;
00448 }
00449 
00450 
00451 template<class ScalarType, class MV, class OP>
00452 void GmresPolySolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00453 {
00454 
00455   // Create the internal parameter list if ones doesn't already exist.
00456   if (params_ == Teuchos::null) {
00457     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00458   }
00459   else {
00460     params->validateParameters(*getValidParameters());
00461   }
00462 
00463   // Check for maximum number of restarts
00464   if (params->isParameter("Maximum Restarts")) {
00465     maxRestarts_ = params->get("Maximum Restarts",maxRestarts_default_);
00466 
00467     // Update parameter in our list.
00468     params_->set("Maximum Restarts", maxRestarts_);
00469   }
00470 
00471   // Check for maximum number of iterations
00472   if (params->isParameter("Maximum Iterations")) {
00473     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00474 
00475     // Update parameter in our list and in status test.
00476     params_->set("Maximum Iterations", maxIters_);
00477     if (maxIterTest_!=Teuchos::null)
00478       maxIterTest_->setMaxIters( maxIters_ );
00479   }
00480 
00481   // Check for blocksize
00482   if (params->isParameter("Block Size")) {
00483     blockSize_ = params->get("Block Size",blockSize_default_);    
00484     TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00485       "Belos::GmresPolySolMgr: \"Block Size\" must be strictly positive.");
00486 
00487     // Update parameter in our list.
00488     params_->set("Block Size", blockSize_);
00489   }
00490 
00491   // Check for the maximum number of blocks.
00492   if (params->isParameter("Num Blocks")) {
00493     numBlocks_ = params->get("Num Blocks",numBlocks_default_);
00494     TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00495       "Belos::GmresPolySolMgr: \"Num Blocks\" must be strictly positive.");
00496 
00497     // Update parameter in our list.
00498     params_->set("Num Blocks", numBlocks_);
00499   }
00500 
00501   // Check to see if the timer label changed.
00502   if (params->isParameter("Timer Label")) {
00503     std::string tempLabel = params->get("Timer Label", label_default_);
00504 
00505     // Update parameter in our list and solver timer
00506     if (tempLabel != label_) {
00507       label_ = tempLabel;
00508       params_->set("Timer Label", label_);
00509       std::string solveLabel = label_ + ": GmresPolySolMgr total solve time";
00510       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00511     }
00512   }
00513 
00514   // Check if the orthogonalization changed.
00515   if (params->isParameter("Orthogonalization")) {
00516     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00517     TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00518                         std::invalid_argument,
00519       "Belos::GmresPolySolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00520     if (tempOrthoType != orthoType_) {
00521       orthoType_ = tempOrthoType;
00522       // Create orthogonalization manager
00523       if (orthoType_=="DGKS") {
00524         if (orthoKappa_ <= 0) {
00525           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00526         }
00527         else {
00528           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00529           Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00530         }
00531       }
00532       else if (orthoType_=="ICGS") {
00533         ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00534       } 
00535       else if (orthoType_=="IMGS") {
00536         ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00537       } 
00538     }  
00539   }
00540 
00541   // Check which orthogonalization constant to use.
00542   if (params->isParameter("Orthogonalization Constant")) {
00543     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00544 
00545     // Update parameter in our list.
00546     params_->set("Orthogonalization Constant",orthoKappa_);
00547     if (orthoType_=="DGKS") {
00548       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00549         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00550       }
00551     } 
00552   }
00553 
00554   // Check for a change in verbosity level
00555   if (params->isParameter("Verbosity")) {
00556     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00557       verbosity_ = params->get("Verbosity", verbosity_default_);
00558     } else {
00559       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00560     }
00561 
00562     // Update parameter in our list.
00563     params_->set("Verbosity", verbosity_);
00564     if (printer_ != Teuchos::null)
00565       printer_->setVerbosity(verbosity_);
00566   }
00567 
00568   // output stream
00569   if (params->isParameter("Output Stream")) {
00570     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00571 
00572     // Update parameter in our list.
00573     params_->set("Output Stream", outputStream_);
00574     if (printer_ != Teuchos::null)
00575       printer_->setOStream( outputStream_ );
00576   }
00577 
00578   // frequency level
00579   if (verbosity_ & Belos::StatusTestDetails) {
00580     if (params->isParameter("Output Frequency")) {
00581       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00582     }
00583 
00584     // Update parameter in out list and output status test.
00585     params_->set("Output Frequency", outputFreq_);
00586     if (outputTest_ != Teuchos::null)
00587       outputTest_->setOutputFrequency( outputFreq_ );
00588   }
00589 
00590   // Create output manager if we need to.
00591   if (printer_ == Teuchos::null) {
00592     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00593   }  
00594   
00595   // Convergence
00596   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00597   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00598 
00599   // Check for convergence tolerance
00600   if (params->isParameter("Convergence Tolerance")) {
00601     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00602 
00603     // Update parameter in our list and residual tests.
00604     params_->set("Convergence Tolerance", convtol_);
00605     if (impConvTest_ != Teuchos::null)
00606       impConvTest_->setTolerance( convtol_ );
00607     if (expConvTest_ != Teuchos::null)
00608       expConvTest_->setTolerance( convtol_ );
00609   }
00610  
00611   // Check for a change in scaling, if so we need to build new residual tests.
00612   if (params->isParameter("Implicit Residual Scaling")) {
00613     std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
00614 
00615     // Only update the scaling if it's different.
00616     if (impResScale_ != tempImpResScale) {
00617       Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
00618       impResScale_ = tempImpResScale;
00619 
00620       // Update parameter in our list and residual tests
00621       params_->set("Implicit Residual Scaling", impResScale_);
00622       if (impConvTest_ != Teuchos::null) {
00623         try { 
00624           impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
00625         }
00626         catch (std::exception& e) { 
00627           // Make sure the convergence test gets constructed again.
00628           isSTSet_ = false;
00629         }
00630       }
00631     }      
00632   }
00633   
00634   if (params->isParameter("Explicit Residual Scaling")) {
00635     std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
00636 
00637     // Only update the scaling if it's different.
00638     if (expResScale_ != tempExpResScale) {
00639       Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
00640       expResScale_ = tempExpResScale;
00641 
00642       // Update parameter in our list and residual tests
00643       params_->set("Explicit Residual Scaling", expResScale_);
00644       if (expConvTest_ != Teuchos::null) {
00645         try { 
00646           expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
00647         }
00648         catch (std::exception& e) {
00649           // Make sure the convergence test gets constructed again.
00650           isSTSet_ = false;
00651         }
00652       }
00653     }      
00654   }
00655 
00656 
00657   if (params->isParameter("Show Maximum Residual Norm Only")) {
00658     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00659 
00660     // Update parameter in our list and residual tests
00661     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00662     if (impConvTest_ != Teuchos::null)
00663       impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00664     if (expConvTest_ != Teuchos::null)
00665       expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00666   }
00667 
00668   // Create orthogonalization manager if we need to.
00669   if (ortho_ == Teuchos::null) {
00670     if (orthoType_=="DGKS") {
00671       if (orthoKappa_ <= 0) {
00672         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00673       }
00674       else {
00675         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00676         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00677       }
00678     }
00679     else if (orthoType_=="ICGS") {
00680       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00681     } 
00682     else if (orthoType_=="IMGS") {
00683       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00684     } 
00685     else {
00686       TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00687         "Belos::GmresPolySolMgr(): Invalid orthogonalization type.");
00688     }  
00689   }
00690 
00691   // Create the timer if we need to.
00692   if (timerSolve_ == Teuchos::null) {
00693     std::string solveLabel = label_ + ": GmresPolySolMgr total solve time";
00694     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00695   }
00696 
00697   // Inform the solver manager that the current parameters were set.
00698   isSet_ = true;
00699 }
00700 
00701 // Check the status test versus the defined linear problem
00702 template<class ScalarType, class MV, class OP>
00703 bool GmresPolySolMgr<ScalarType,MV,OP>::checkStatusTest() {
00704 
00705   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00706   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestGenResNorm_t;
00707   typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP>  StatusTestImpResNorm_t;
00708 
00709   // Basic test checks maximum iterations and native residual.
00710   maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00711 
00712   // If there is a left preconditioner, we create a combined status test that checks the implicit
00713   // and then explicit residual norm to see if we have convergence.
00714   if (!Teuchos::is_null(problem_->getLeftPrec())) {
00715     expResTest_ = true;
00716   }
00717 
00718   if (expResTest_) {
00719    
00720     // Implicit residual test, using the native residual to determine if convergence was achieved.
00721     Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
00722       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
00723     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00724     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00725     impConvTest_ = tmpImpConvTest;
00726 
00727     // Explicit residual test once the native residual is below the tolerance
00728     Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
00729       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
00730     tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
00731     tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
00732     tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00733     expConvTest_ = tmpExpConvTest;
00734 
00735     // The convergence test is a combination of the "cheap" implicit test and explicit test.
00736     convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
00737   }
00738   else {
00739 
00740     // Implicit residual test, using the native residual to determine if convergence was achieved.
00741     // Use test that checks for loss of accuracy.
00742     Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
00743       Teuchos::rcp( new StatusTestImpResNorm_t( convtol_ ) );
00744     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00745     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00746     impConvTest_ = tmpImpConvTest;
00747 
00748     // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
00749     expConvTest_ = impConvTest_;
00750     convTest_ = impConvTest_;
00751   }
00752 
00753   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00754 
00755   if (outputFreq_ > 0) {
00756     outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,
00757         sTest_,
00758         outputFreq_,
00759         Passed+Failed+Undefined ) );
00760   }
00761   else {
00762     outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,
00763         sTest_, 1 ) );
00764   }
00765 
00766   // The status test is now set.
00767   isSTSet_ = true;
00768 
00769   return false;
00770 }
00771   
00772 // solve()
00773 template<class ScalarType, class MV, class OP>
00774 ReturnType GmresPolySolMgr<ScalarType,MV,OP>::solve() {
00775 
00776   // Set the current parameters if they were not set before.
00777   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00778   // then didn't set any parameters using setParameters().
00779   if (!isSet_) {
00780     setParameters(Teuchos::parameterList(*getValidParameters()));
00781   }
00782 
00783   Teuchos::BLAS<int,ScalarType> blas;
00784   Teuchos::LAPACK<int,ScalarType> lapack;
00785   
00786   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,GmresPolySolMgrLinearProblemFailure,
00787     "Belos::GmresPolySolMgr::solve(): Linear problem is not a valid object.");
00788 
00789   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),GmresPolySolMgrLinearProblemFailure,
00790     "Belos::GmresPolySolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00791 
00792   if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
00793     TEST_FOR_EXCEPTION( checkStatusTest(),GmresPolySolMgrLinearProblemFailure,
00794       "Belos::GmresPolySolMgr::solve(): Linear problem and requested status tests are incompatible.");
00795   }
00796 
00797   // Create indices for the linear systems to be solved.
00798   int startPtr = 0;
00799   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00800   int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
00801 
00802   std::vector<int> currIdx;
00803   //  If an adaptive block size is allowed then only the linear systems that need to be solved are solved.
00804   //  Otherwise, the index set is generated that informs the linear problem that some linear systems are augmented.
00805   currIdx.resize( blockSize_ );
00806   for (int i=0; i<numCurrRHS; ++i) 
00807     { currIdx[i] = startPtr+i; }
00808   for (int i=numCurrRHS; i<blockSize_; ++i) 
00809     { currIdx[i] = -1; }
00810 
00811   // Inform the linear problem of the current linear system to solve.
00812   problem_->setLSIndex( currIdx );
00813 
00815   // Parameter list
00816   Teuchos::ParameterList plist;
00817   plist.set("Block Size",blockSize_);
00818 
00819   int dim = MVT::GetVecLength( *(problem_->getRHS()) );  
00820   if (blockSize_*numBlocks_ > dim) {
00821     int tmpNumBlocks = 0;
00822     if (blockSize_ == 1)
00823       tmpNumBlocks = dim / blockSize_;  // Allow for a good breakdown.
00824     else
00825       tmpNumBlocks = ( dim - blockSize_) / blockSize_;  // Allow for restarting.
00826     printer_->stream(Warnings) << 
00827       "Warning! Requested Krylov subspace dimension is larger that operator dimension!" << std::endl <<
00828       " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
00829     plist.set("Num Blocks",tmpNumBlocks);
00830   } 
00831   else 
00832     plist.set("Num Blocks",numBlocks_);
00833   
00834   // Reset the status test.  
00835   outputTest_->reset();
00836   loaDetected_ = false;
00837 
00838   // Assume convergence is achieved, then let any failed convergence set this to false.
00839   bool isConverged = true;  
00840 
00842   // BlockGmres solver
00843 
00844   Teuchos::RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter;
00845 
00846   block_gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
00847   
00848   // Enter solve() iterations
00849   {
00850     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00851     
00852     while ( numRHS2Solve > 0 ) {
00853       
00854       // Set the current number of blocks and blocksize with the Gmres iteration.
00855       if (blockSize_*numBlocks_ > dim) {
00856         int tmpNumBlocks = 0;
00857         if (blockSize_ == 1)
00858           tmpNumBlocks = dim / blockSize_;  // Allow for a good breakdown.
00859         else
00860           tmpNumBlocks = ( dim - blockSize_) / blockSize_;  // Allow for restarting.
00861         block_gmres_iter->setSize( blockSize_, tmpNumBlocks );
00862       }
00863       else
00864         block_gmres_iter->setSize( blockSize_, numBlocks_ );
00865       
00866       // Reset the number of iterations.
00867       block_gmres_iter->resetNumIters();
00868       
00869       // Reset the number of calls that the status test output knows about.
00870       outputTest_->resetNumCalls();
00871       
00872       // Create the first block in the current Krylov basis.
00873       Teuchos::RCP<MV> V_0;
00874       V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
00875       
00876       
00877       // Get a matrix to hold the orthonormalization coefficients.
00878       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_0 = 
00879         rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
00880       
00881       // Orthonormalize the new V_0
00882       int rank = ortho_->normalize( *V_0, z_0 );
00883       TEST_FOR_EXCEPTION(rank != blockSize_,GmresPolySolMgrOrthoFailure,
00884        "Belos::GmresPolySolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
00885       
00886       // Set the new state and initialize the solver.
00887       GmresIterationState<ScalarType,MV> newstate;
00888       newstate.V = V_0;
00889       newstate.z = z_0;
00890       newstate.curDim = 0;
00891       block_gmres_iter->initializeGmres(newstate);
00892       int numRestarts = 0;
00893       
00894       while(1) {
00895         // tell block_gmres_iter to iterate
00896         try {
00897           block_gmres_iter->iterate();
00898     
00900           //
00901           // check convergence first
00902           //
00904           if ( convTest_->getStatus() == Passed ) {
00905             if ( expConvTest_->getLOADetected() ) {
00906               // we don't have convergence
00907               loaDetected_ = true;
00908               isConverged = false;
00909             }
00910             break;  // break from while(1){block_gmres_iter->iterate()}
00911           }
00913           //
00914           // check for maximum iterations
00915           //
00917           else if ( maxIterTest_->getStatus() == Passed ) {
00918             // we don't have convergence
00919             isConverged = false;
00920             break;  // break from while(1){block_gmres_iter->iterate()}
00921           }
00923           //
00924           // check for restarting, i.e. the subspace is full
00925           //
00927           else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
00928       
00929             if ( numRestarts >= maxRestarts_ ) {
00930               isConverged = false;
00931               break; // break from while(1){block_gmres_iter->iterate()}
00932             }
00933             numRestarts++;
00934       
00935             printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
00936       
00937             // Update the linear problem.
00938             Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
00939       problem_->updateSolution( update, true );
00940 
00941             // Get the state.
00942             GmresIterationState<ScalarType,MV> oldState = block_gmres_iter->getState();
00943       
00944             // Compute the restart std::vector.
00945             // Get a view of the current Krylov basis.
00946             Teuchos::RCP<MV> V_0  = MVT::Clone( *(oldState.V), blockSize_ );
00947       problem_->computeCurrPrecResVec( &*V_0 );
00948       
00949             // Get a view of the first block of the Krylov basis.
00950             Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_0 = 
00951               rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
00952       
00953             // Orthonormalize the new V_0
00954             int rank = ortho_->normalize( *V_0, z_0 );
00955             TEST_FOR_EXCEPTION(rank != blockSize_,GmresPolySolMgrOrthoFailure,
00956              "Belos::GmresPolySolMgr::solve(): Failed to compute initial block of orthonormal vectors after restart.");
00957       
00958             // Set the new state and initialize the solver.
00959             GmresIterationState<ScalarType,MV> newstate;
00960             newstate.V = V_0;
00961             newstate.z = z_0;
00962             newstate.curDim = 0;
00963             block_gmres_iter->initializeGmres(newstate);
00964       
00965           } // end of restarting
00966     
00968           //
00969           // we returned from iterate(), but none of our status tests Passed.
00970           // something is wrong, and it is probably our fault.
00971           //
00973     
00974           else {
00975             TEST_FOR_EXCEPTION(true,std::logic_error,
00976              "Belos::GmresPolySolMgr::solve(): Invalid return from BlockGmresIter::iterate().");
00977           }
00978         }
00979         catch (GmresIterationOrthoFailure e) {
00980           // If the block size is not one, it's not considered a lucky breakdown.
00981           if (blockSize_ != 1) {
00982             printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 
00983                                      << block_gmres_iter->getNumIters() << std::endl 
00984                                      << e.what() << std::endl;
00985             if (convTest_->getStatus() != Passed)
00986               isConverged = false;
00987             break;
00988           } 
00989           else {
00990             // If the block size is one, try to recover the most recent least-squares solution
00991             block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
00992       
00993             // Check to see if the most recent least-squares solution yielded convergence.
00994             sTest_->checkStatus( &*block_gmres_iter );
00995             if (convTest_->getStatus() != Passed)
00996               isConverged = false;
00997             break;
00998           }
00999         }
01000         catch (std::exception e) {
01001           printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 
01002                                    << block_gmres_iter->getNumIters() << std::endl 
01003                                    << e.what() << std::endl;
01004           throw;
01005         }
01006       }
01007       
01008       // Compute the current solution.
01009       // Update the linear problem.
01010       // Attempt to get the current solution from the residual status test, if it has one.
01011       if ( !Teuchos::is_null(expConvTest_->getSolution()) ) {
01012   Teuchos::RCP<MV> newX = expConvTest_->getSolution();
01013   Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
01014   MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
01015       }
01016       else {
01017   Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01018   problem_->updateSolution( update, true );
01019       }  
01020       
01021       // Inform the linear problem that we are finished with this block linear system.
01022       problem_->setCurrLS();
01023       
01024       // Update indices for the linear systems to be solved.
01025       startPtr += numCurrRHS;
01026       numRHS2Solve -= numCurrRHS;
01027       if ( numRHS2Solve > 0 ) {
01028         numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01029   
01030   currIdx.resize( blockSize_ );
01031   for (int i=0; i<numCurrRHS; ++i) 
01032           { currIdx[i] = startPtr+i; }
01033   for (int i=numCurrRHS; i<blockSize_; ++i) 
01034           { currIdx[i] = -1; }
01035 
01036         // Set the next indices.
01037         problem_->setLSIndex( currIdx );
01038       }
01039       else {
01040         currIdx.resize( numRHS2Solve );
01041       }
01042       
01043     }// while ( numRHS2Solve > 0 )
01044     
01045   }
01046   
01047   // print final summary
01048   sTest_->print( printer_->stream(FinalSummary) );
01049   
01050   // print timing information
01051   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01052   
01053   if (!isConverged || loaDetected_) {
01054     return Unconverged; // return from GmresPolySolMgr::solve() 
01055   }
01056   return Converged; // return from GmresPolySolMgr::solve() 
01057 }
01058   
01059 //  This method requires the solver manager to return a std::string that describes itself.
01060 template<class ScalarType, class MV, class OP>
01061 std::string GmresPolySolMgr<ScalarType,MV,OP>::description() const
01062   {
01063   std::ostringstream oss;
01064   oss << "Belos::GmresPolySolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01065   oss << "{";
01066   oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" <<blockSize_;
01067   oss << ", Num Blocks=" <<numBlocks_<< ", Max Restarts=" << maxRestarts_;
01068   oss << "}";
01069   return oss.str();
01070 }
01071   
01072 } // end Belos namespace
01073 
01074 #endif /* BELOS_GMRES_POLY_SOLMGR_HPP */

Generated on Wed May 12 21:45:51 2010 for Belos by  doxygen 1.4.7