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 "BelosGmresPolyOp.hpp"
00043 #include "BelosGmresIteration.hpp"
00044 #include "BelosBlockGmresIter.hpp"
00045 #include "BelosDGKSOrthoManager.hpp"
00046 #include "BelosICGSOrthoManager.hpp"
00047 #include "BelosIMGSOrthoManager.hpp"
00048 #include "BelosStatusTestMaxIters.hpp"
00049 #include "BelosStatusTestGenResNorm.hpp"
00050 #include "BelosStatusTestImpResNorm.hpp"
00051 #include "BelosStatusTestCombo.hpp"
00052 #include "BelosStatusTestOutput.hpp"
00053 #include "BelosOutputManager.hpp"
00054 #include "Teuchos_BLAS.hpp"
00055 #include "Teuchos_LAPACK.hpp"
00056 #include "Teuchos_TimeMonitor.hpp"
00057 
00073 namespace Belos {
00074   
00076 
00077   
00084 class GmresPolySolMgrLinearProblemFailure : public BelosError {public:
00085   GmresPolySolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00086     {}};
00087   
00094 class GmresPolySolMgrPolynomialFailure : public BelosError {public:
00095   GmresPolySolMgrPolynomialFailure(const std::string& what_arg) : BelosError(what_arg)
00096     {}};
00097   
00104 class GmresPolySolMgrOrthoFailure : public BelosError {public:
00105   GmresPolySolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00106     {}};
00107   
00108 template<class ScalarType, class MV, class OP>
00109 class GmresPolySolMgr : public SolverManager<ScalarType,MV,OP> {
00110     
00111 private:
00112   typedef MultiVecTraits<ScalarType,MV> MVT;
00113   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00114   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00115   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00116   typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00117     
00118 public:
00119     
00121 
00122    
00128   GmresPolySolMgr();
00129  
00142   GmresPolySolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00143     const Teuchos::RCP<Teuchos::ParameterList> &pl );
00144     
00146   virtual ~GmresPolySolMgr() {};
00148     
00150 
00151     
00154   const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00155     return *problem_;
00156   }
00157 
00160   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00161 
00164   Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00165  
00171   Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00172     return tuple(timerSolve_, timerPoly_);
00173   }
00174   
00176   int getNumIters() const {
00177     return numIters_;
00178   }
00179  
00183   bool isLOADetected() const { return loaDetected_; }
00184  
00186     
00188 
00189   
00191   void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; isSTSet_ = false; }
00192 
00194   void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00195     
00197   
00199 
00200 
00204   void reset( const ResetType type ) { 
00205     if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) {
00206       problem_->setProblem(); 
00207       isPolyBuilt_ = false;  // Rebuild the GMRES polynomial
00208     }
00209   }
00211   
00213 
00214     
00232   ReturnType solve();
00233     
00235     
00238     
00240   std::string description() const;
00241     
00243     
00244 private:
00245 
00246   // Method to convert std::string to enumerated type for residual.
00247   Belos::ScaleType convertStringToScaleType( std::string& scaleType ) {
00248     if (scaleType == "Norm of Initial Residual") {
00249       return Belos::NormOfInitRes;
00250     } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00251       return Belos::NormOfPrecInitRes;
00252     } else if (scaleType == "Norm of RHS") {
00253       return Belos::NormOfRHS;
00254     } else if (scaleType == "None") {
00255       return Belos::None;
00256     } else 
00257       TEST_FOR_EXCEPTION( true ,std::logic_error,
00258         "Belos::GmresPolySolMgr(): Invalid residual scaling type.");
00259   }
00260   
00261   // Method for checking current status test against defined linear problem.
00262   bool checkStatusTest();
00263 
00264   // Method for generating GMRES polynomial.
00265   bool generatePoly();
00266 
00267   // Linear problem.
00268   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00269     
00270   // Output manager.
00271   Teuchos::RCP<OutputManager<ScalarType> > printer_;
00272   Teuchos::RCP<std::ostream> outputStream_;
00273 
00274   // Status test.
00275   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00276   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00277   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00278   Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
00279   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00280 
00281   // Orthogonalization manager.
00282   Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00283     
00284   // Current parameter list.
00285   Teuchos::RCP<ParameterList> params_;
00286 
00287   // Default solver values.
00288   static const MagnitudeType polytol_default_;
00289   static const MagnitudeType convtol_default_;
00290   static const MagnitudeType orthoKappa_default_;
00291   static const int maxDegree_default_;
00292   static const int maxRestarts_default_;
00293   static const int maxIters_default_;
00294   static const bool strictConvTol_default_;
00295   static const bool showMaxResNormOnly_default_;
00296   static const int blockSize_default_;
00297   static const int numBlocks_default_;
00298   static const int verbosity_default_;
00299   static const int outputFreq_default_;
00300   static const std::string impResScale_default_; 
00301   static const std::string expResScale_default_; 
00302   static const std::string label_default_;
00303   static const std::string orthoType_default_;
00304   static const Teuchos::RCP<std::ostream> outputStream_default_;
00305 
00306   // Current solver values.
00307   MagnitudeType polytol_, convtol_, orthoKappa_;
00308   int maxDegree_, maxRestarts_, maxIters_, numIters_;
00309   int blockSize_, numBlocks_, verbosity_, outputFreq_;
00310   bool strictConvTol_, showMaxResNormOnly_;
00311   std::string orthoType_; 
00312   std::string impResScale_, expResScale_;
00313 
00314   // Polynomial storage
00315   int poly_dim_;
00316   Teuchos::RCP<Teuchos::SerialDenseMatrix<int, ScalarType> > poly_H_, poly_y_;
00317   Teuchos::RCP<Teuchos::SerialDenseVector<int, ScalarType> > poly_r0_;
00318   Teuchos::RCP<Belos::GmresPolyOp<ScalarType, MV, OP> > poly_Op_;  
00319   
00320   // Timers.
00321   std::string label_;
00322   Teuchos::RCP<Teuchos::Time> timerSolve_, timerPoly_;
00323 
00324   // Internal state variables.
00325   bool isPolyBuilt_;  
00326   bool isSet_, isSTSet_, expResTest_;
00327   bool loaDetected_;
00328 };
00329 
00330 
00331 // Default solver values.
00332 template<class ScalarType, class MV, class OP>
00333 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType GmresPolySolMgr<ScalarType,MV,OP>::polytol_default_ = 1e-12;
00334 
00335 template<class ScalarType, class MV, class OP>
00336 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType GmresPolySolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00337 
00338 template<class ScalarType, class MV, class OP>
00339 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType GmresPolySolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00340 
00341 template<class ScalarType, class MV, class OP>
00342 const int GmresPolySolMgr<ScalarType,MV,OP>::maxDegree_default_ = 25;
00343 
00344 template<class ScalarType, class MV, class OP>
00345 const int GmresPolySolMgr<ScalarType,MV,OP>::maxRestarts_default_ = 20;
00346 
00347 template<class ScalarType, class MV, class OP>
00348 const int GmresPolySolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00349 
00350 template<class ScalarType, class MV, class OP>
00351 const bool GmresPolySolMgr<ScalarType,MV,OP>::strictConvTol_default_ = false;
00352 
00353 template<class ScalarType, class MV, class OP>
00354 const bool GmresPolySolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00355 
00356 template<class ScalarType, class MV, class OP>
00357 const int GmresPolySolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00358 
00359 template<class ScalarType, class MV, class OP>
00360 const int GmresPolySolMgr<ScalarType,MV,OP>::numBlocks_default_ = 300;
00361 
00362 template<class ScalarType, class MV, class OP>
00363 const int GmresPolySolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00364 
00365 template<class ScalarType, class MV, class OP>
00366 const int GmresPolySolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00367 
00368 template<class ScalarType, class MV, class OP>
00369 const std::string GmresPolySolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of RHS";
00370 
00371 template<class ScalarType, class MV, class OP>
00372 const std::string GmresPolySolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of RHS";
00373 
00374 template<class ScalarType, class MV, class OP>
00375 const std::string GmresPolySolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00376 
00377 template<class ScalarType, class MV, class OP>
00378 const std::string GmresPolySolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00379 
00380 template<class ScalarType, class MV, class OP>
00381 const Teuchos::RCP<std::ostream> GmresPolySolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00382 
00383 
00384 // Empty Constructor
00385 template<class ScalarType, class MV, class OP>
00386 GmresPolySolMgr<ScalarType,MV,OP>::GmresPolySolMgr() :
00387   outputStream_(outputStream_default_),
00388   polytol_(polytol_default_),
00389   convtol_(convtol_default_),
00390   orthoKappa_(orthoKappa_default_),
00391   maxDegree_(maxDegree_default_),
00392   maxRestarts_(maxRestarts_default_),
00393   maxIters_(maxIters_default_),
00394   blockSize_(blockSize_default_),
00395   numBlocks_(numBlocks_default_),
00396   verbosity_(verbosity_default_),
00397   outputFreq_(outputFreq_default_),
00398   strictConvTol_(strictConvTol_default_),
00399   showMaxResNormOnly_(showMaxResNormOnly_default_),
00400   orthoType_(orthoType_default_),
00401   impResScale_(impResScale_default_),
00402   expResScale_(expResScale_default_),
00403   label_(label_default_),
00404   isPolyBuilt_(false),
00405   isSet_(false),
00406   isSTSet_(false),
00407   expResTest_(false),
00408   loaDetected_(false)
00409 {}
00410 
00411 
00412 // Basic Constructor
00413 template<class ScalarType, class MV, class OP>
00414 GmresPolySolMgr<ScalarType,MV,OP>::GmresPolySolMgr( 
00415   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00416   const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00417   problem_(problem),
00418   outputStream_(outputStream_default_),
00419   polytol_(polytol_default_),
00420   convtol_(convtol_default_),
00421   orthoKappa_(orthoKappa_default_),
00422   maxDegree_(maxDegree_default_),
00423   maxRestarts_(maxRestarts_default_),
00424   maxIters_(maxIters_default_),
00425   blockSize_(blockSize_default_),
00426   numBlocks_(numBlocks_default_),
00427   verbosity_(verbosity_default_),
00428   outputFreq_(outputFreq_default_),
00429   strictConvTol_(strictConvTol_default_), 
00430   showMaxResNormOnly_(showMaxResNormOnly_default_),
00431   orthoType_(orthoType_default_),
00432   impResScale_(impResScale_default_),
00433   expResScale_(expResScale_default_),
00434   label_(label_default_),
00435   isPolyBuilt_(false),
00436   isSet_(false),
00437   isSTSet_(false),
00438   expResTest_(false),
00439   loaDetected_(false)
00440 {
00441 
00442   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00443 
00444   // If the parameter list pointer is null, then set the current parameters to the default parameter list.
00445   if ( !is_null(pl) ) {
00446     setParameters( pl );  
00447   }
00448 
00449 }
00450 
00451 
00452 template<class ScalarType, class MV, class OP>
00453 Teuchos::RCP<const Teuchos::ParameterList>
00454 GmresPolySolMgr<ScalarType,MV,OP>::getValidParameters() const
00455 {
00456   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00457   if (is_null(validPL)) {
00458     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00459     pl->set("Polynomial Tolerance", polytol_default_,
00460       "The relative residual tolerance that used to construct the GMRES polynomial.");
00461     pl->set("Maximum Degree", maxDegree_default_,
00462       "The maximum degree allowed for any GMRES polynomial.");
00463     pl->set("Convergence Tolerance", convtol_default_,
00464       "The relative residual tolerance that needs to be achieved by the\n"
00465       "iterative solver in order for the linear system to be declared converged." );
00466     pl->set("Maximum Restarts", maxRestarts_default_,
00467       "The maximum number of restarts allowed for each\n"
00468       "set of RHS solved.");
00469     pl->set("Maximum Iterations", maxIters_default_,
00470       "The maximum number of block iterations allowed for each\n"
00471       "set of RHS solved.");
00472     pl->set("Num Blocks", numBlocks_default_,
00473       "The maximum number of blocks allowed in the Krylov subspace\n"
00474       "for each set of RHS solved.");
00475     pl->set("Block Size", blockSize_default_,
00476       "The number of vectors in each block.  This number times the\n"
00477       "number of blocks is the total Krylov subspace dimension.");
00478     pl->set("Verbosity", verbosity_default_,
00479       "What type(s) of solver information should be outputted\n"
00480       "to the output stream.");
00481     pl->set("Output Frequency", outputFreq_default_,
00482       "How often convergence information should be outputted\n"
00483       "to the output stream.");  
00484     pl->set("Output Stream", outputStream_default_,
00485       "A reference-counted pointer to the output stream where all\n"
00486       "solver output is sent.");
00487     pl->set("Strict Convergence", strictConvTol_default_,
00488       "After polynomial is applied, whether solver should try to achieve\n"
00489       "the relative residual tolerance.");
00490     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00491       "When convergence information is printed, only show the maximum\n"
00492       "relative residual norm when the block size is greater than one.");
00493     pl->set("Implicit Residual Scaling", impResScale_default_,
00494       "The type of scaling used in the implicit residual convergence test.");
00495     pl->set("Explicit Residual Scaling", expResScale_default_,
00496       "The type of scaling used in the explicit residual convergence test.");
00497     pl->set("Timer Label", label_default_,
00498       "The string to use as a prefix for the timer labels.");
00499     //  pl->set("Restart Timers", restartTimers_);
00500     pl->set("Orthogonalization", orthoType_default_,
00501       "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
00502     pl->set("Orthogonalization Constant",orthoKappa_default_,
00503       "The constant used by DGKS orthogonalization to determine\n"
00504       "whether another step of classical Gram-Schmidt is necessary.");
00505     validPL = pl;
00506   }
00507   return validPL;
00508 }
00509 
00510 
00511 template<class ScalarType, class MV, class OP>
00512 void GmresPolySolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00513 {
00514 
00515   // Create the internal parameter list if ones doesn't already exist.
00516   if (params_ == Teuchos::null) {
00517     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00518   }
00519   else {
00520     params->validateParameters(*getValidParameters());
00521   }
00522 
00523   // Check for maximum polynomial degree
00524   if (params->isParameter("Maximum Degree")) {
00525     maxDegree_ = params->get("Maximum Degree",maxDegree_default_);
00526 
00527     // Update parameter in our list.
00528     params_->set("Maximum Degree", maxDegree_);
00529   }
00530 
00531   // Check for maximum number of restarts
00532   if (params->isParameter("Maximum Restarts")) {
00533     maxRestarts_ = params->get("Maximum Restarts",maxRestarts_default_);
00534 
00535     // Update parameter in our list.
00536     params_->set("Maximum Restarts", maxRestarts_);
00537   }
00538 
00539   // Check for maximum number of iterations
00540   if (params->isParameter("Maximum Iterations")) {
00541     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00542 
00543     // Update parameter in our list and in status test.
00544     params_->set("Maximum Iterations", maxIters_);
00545     if (maxIterTest_!=Teuchos::null)
00546       maxIterTest_->setMaxIters( maxIters_ );
00547   }
00548 
00549   // Check for blocksize
00550   if (params->isParameter("Block Size")) {
00551     blockSize_ = params->get("Block Size",blockSize_default_);    
00552     TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00553       "Belos::GmresPolySolMgr: \"Block Size\" must be strictly positive.");
00554 
00555     // Update parameter in our list.
00556     params_->set("Block Size", blockSize_);
00557   }
00558 
00559   // Check for the maximum number of blocks.
00560   if (params->isParameter("Num Blocks")) {
00561     numBlocks_ = params->get("Num Blocks",numBlocks_default_);
00562     TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00563       "Belos::GmresPolySolMgr: \"Num Blocks\" must be strictly positive.");
00564 
00565     // Update parameter in our list.
00566     params_->set("Num Blocks", numBlocks_);
00567   }
00568 
00569   // Check to see if the timer label changed.
00570   if (params->isParameter("Timer Label")) {
00571     std::string tempLabel = params->get("Timer Label", label_default_);
00572 
00573     // Update parameter in our list and solver timer
00574     if (tempLabel != label_) {
00575       label_ = tempLabel;
00576       params_->set("Timer Label", label_);
00577       std::string solveLabel = label_ + ": GmresPolySolMgr total solve time";
00578       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00579       std::string polyLabel = label_ + ": GmresPolySolMgr polynomial creation time";
00580       timerPoly_ = Teuchos::TimeMonitor::getNewTimer(polyLabel);
00581     }
00582   }
00583 
00584   // Check if the orthogonalization changed.
00585   if (params->isParameter("Orthogonalization")) {
00586     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00587     TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00588                         std::invalid_argument,
00589       "Belos::GmresPolySolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00590     if (tempOrthoType != orthoType_) {
00591       orthoType_ = tempOrthoType;
00592       // Create orthogonalization manager
00593       if (orthoType_=="DGKS") {
00594         if (orthoKappa_ <= 0) {
00595           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00596         }
00597         else {
00598           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00599           Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00600         }
00601       }
00602       else if (orthoType_=="ICGS") {
00603         ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00604       } 
00605       else if (orthoType_=="IMGS") {
00606         ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00607       } 
00608     }  
00609   }
00610 
00611   // Check which orthogonalization constant to use.
00612   if (params->isParameter("Orthogonalization Constant")) {
00613     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00614 
00615     // Update parameter in our list.
00616     params_->set("Orthogonalization Constant",orthoKappa_);
00617     if (orthoType_=="DGKS") {
00618       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00619         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00620       }
00621     } 
00622   }
00623 
00624   // Check for a change in verbosity level
00625   if (params->isParameter("Verbosity")) {
00626     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00627       verbosity_ = params->get("Verbosity", verbosity_default_);
00628     } else {
00629       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00630     }
00631 
00632     // Update parameter in our list.
00633     params_->set("Verbosity", verbosity_);
00634     if (printer_ != Teuchos::null)
00635       printer_->setVerbosity(verbosity_);
00636   }
00637 
00638   // output stream
00639   if (params->isParameter("Output Stream")) {
00640     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00641 
00642     // Update parameter in our list.
00643     params_->set("Output Stream", outputStream_);
00644     if (printer_ != Teuchos::null)
00645       printer_->setOStream( outputStream_ );
00646   }
00647 
00648   // frequency level
00649   if (verbosity_ & Belos::StatusTestDetails) {
00650     if (params->isParameter("Output Frequency")) {
00651       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00652     }
00653 
00654     // Update parameter in out list and output status test.
00655     params_->set("Output Frequency", outputFreq_);
00656     if (outputTest_ != Teuchos::null)
00657       outputTest_->setOutputFrequency( outputFreq_ );
00658   }
00659 
00660   // Create output manager if we need to.
00661   if (printer_ == Teuchos::null) {
00662     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00663   }  
00664   
00665   // Convergence
00666   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00667   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00668 
00669   // Check for polynomial convergence tolerance
00670   if (params->isParameter("Polynomial Tolerance")) {
00671     polytol_ = params->get("Polynomial Tolerance",polytol_default_);
00672 
00673     // Update parameter in our list and residual tests.
00674     params_->set("Polynomial Tolerance", polytol_);
00675   }
00676 
00677   // Check for convergence tolerance
00678   if (params->isParameter("Convergence Tolerance")) {
00679     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00680 
00681     // Update parameter in our list and residual tests.
00682     params_->set("Convergence Tolerance", convtol_);
00683     if (impConvTest_ != Teuchos::null)
00684       impConvTest_->setTolerance( convtol_ );
00685     if (expConvTest_ != Teuchos::null)
00686       expConvTest_->setTolerance( convtol_ );
00687   }
00688  
00689   // Check if user requires solver to reach convergence tolerance
00690   if (params->isParameter("Strict Convergence")) {
00691     strictConvTol_ = params->get("Strict Convergence",strictConvTol_default_);
00692 
00693     // Update parameter in our list and residual tests
00694     params_->set("Strict Convergence", strictConvTol_);
00695   }
00696  
00697   // Check for a change in scaling, if so we need to build new residual tests.
00698   if (params->isParameter("Implicit Residual Scaling")) {
00699     std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
00700 
00701     // Only update the scaling if it's different.
00702     if (impResScale_ != tempImpResScale) {
00703       Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
00704       impResScale_ = tempImpResScale;
00705 
00706       // Update parameter in our list and residual tests
00707       params_->set("Implicit Residual Scaling", impResScale_);
00708       if (impConvTest_ != Teuchos::null) {
00709         try { 
00710           impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
00711         }
00712         catch (std::exception& e) { 
00713           // Make sure the convergence test gets constructed again.
00714           isSTSet_ = false;
00715         }
00716       }
00717     }      
00718   }
00719   
00720   if (params->isParameter("Explicit Residual Scaling")) {
00721     std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
00722 
00723     // Only update the scaling if it's different.
00724     if (expResScale_ != tempExpResScale) {
00725       Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
00726       expResScale_ = tempExpResScale;
00727 
00728       // Update parameter in our list and residual tests
00729       params_->set("Explicit Residual Scaling", expResScale_);
00730       if (expConvTest_ != Teuchos::null) {
00731         try { 
00732           expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
00733         }
00734         catch (std::exception& e) {
00735           // Make sure the convergence test gets constructed again.
00736           isSTSet_ = false;
00737         }
00738       }
00739     }      
00740   }
00741 
00742 
00743   if (params->isParameter("Show Maximum Residual Norm Only")) {
00744     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00745 
00746     // Update parameter in our list and residual tests
00747     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00748     if (impConvTest_ != Teuchos::null)
00749       impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00750     if (expConvTest_ != Teuchos::null)
00751       expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00752   }
00753 
00754   // Create orthogonalization manager if we need to.
00755   if (ortho_ == Teuchos::null) {
00756     if (orthoType_=="DGKS") {
00757       if (orthoKappa_ <= 0) {
00758         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00759       }
00760       else {
00761         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00762         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00763       }
00764     }
00765     else if (orthoType_=="ICGS") {
00766       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00767     } 
00768     else if (orthoType_=="IMGS") {
00769       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00770     } 
00771     else {
00772       TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00773         "Belos::GmresPolySolMgr(): Invalid orthogonalization type.");
00774     }  
00775   }
00776 
00777   // Create the timers if we need to.
00778   if (timerSolve_ == Teuchos::null) {
00779     std::string solveLabel = label_ + ": GmresPolySolMgr total solve time";
00780     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00781   }
00782   
00783   if (timerPoly_ == Teuchos::null) {
00784     std::string polyLabel = label_ + ": GmresPolySolMgr polynomial creation time";
00785     timerPoly_ = Teuchos::TimeMonitor::getNewTimer(polyLabel);
00786   }
00787 
00788   // Inform the solver manager that the current parameters were set.
00789   isSet_ = true;
00790 }
00791 
00792 // Check the status test versus the defined linear problem
00793 template<class ScalarType, class MV, class OP>
00794 bool GmresPolySolMgr<ScalarType,MV,OP>::checkStatusTest() {
00795 
00796   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00797   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestGenResNorm_t;
00798   typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP>  StatusTestImpResNorm_t;
00799 
00800   // Basic test checks maximum iterations and native residual.
00801   maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00802 
00803   // If there is a left preconditioner, we create a combined status test that checks the implicit
00804   // and then explicit residual norm to see if we have convergence.
00805   if (!Teuchos::is_null(problem_->getLeftPrec())) {
00806     expResTest_ = true;
00807   }
00808 
00809   if (expResTest_) {
00810    
00811     // Implicit residual test, using the native residual to determine if convergence was achieved.
00812     Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
00813       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
00814     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00815     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00816     impConvTest_ = tmpImpConvTest;
00817 
00818     // Explicit residual test once the native residual is below the tolerance
00819     Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
00820       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
00821     tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
00822     tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
00823     tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00824     expConvTest_ = tmpExpConvTest;
00825 
00826     // The convergence test is a combination of the "cheap" implicit test and explicit test.
00827     convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
00828   }
00829   else {
00830 
00831     // Implicit residual test, using the native residual to determine if convergence was achieved.
00832     // Use test that checks for loss of accuracy.
00833     Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
00834       Teuchos::rcp( new StatusTestImpResNorm_t( convtol_ ) );
00835     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00836     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00837     impConvTest_ = tmpImpConvTest;
00838 
00839     // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
00840     expConvTest_ = impConvTest_;
00841     convTest_ = impConvTest_;
00842   }
00843 
00844   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00845 
00846   if (outputFreq_ > 0) {
00847     outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,
00848         sTest_,
00849         outputFreq_,
00850         Passed+Failed+Undefined ) );
00851   }
00852   else {
00853     outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,
00854         sTest_, 1 ) );
00855   }
00856 
00857   // The status test is now set.
00858   isSTSet_ = true;
00859 
00860   return false;
00861 }
00862 
00863 template<class ScalarType, class MV, class OP>
00864 bool GmresPolySolMgr<ScalarType,MV,OP>::generatePoly()
00865 {
00866   // Create a copy of the linear problem that has a zero initial guess and random RHS.
00867   Teuchos::RCP<MV> newX  = MVT::Clone( *(problem_->getLHS()), 1 );
00868   Teuchos::RCP<MV> newB  = MVT::Clone( *(problem_->getRHS()), 1 );
00869   MVT::MvInit( *newX, SCT::zero() );
00870   MVT::MvRandom( *newB );
00871   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > newProblem =
00872     Teuchos::rcp( new LinearProblem<ScalarType,MV,OP>( problem_->getOperator(), newX, newB ) );
00873   newProblem->setLeftPrec( problem_->getLeftPrec() );
00874   newProblem->setRightPrec( problem_->getRightPrec() );
00875   newProblem->setLabel("Belos GMRES Poly Generation");
00876   newProblem->setProblem();
00877   std::vector<int> idx(1,0);       // Must set the index to be the first vector (0)!
00878   newProblem->setLSIndex( idx );
00879 
00880   // Create a parameter list for the GMRES iteration.
00881   Teuchos::ParameterList polyList;
00882 
00883   // Tell the block solver that the block size is one.
00884   polyList.set("Num Blocks",maxDegree_);
00885   polyList.set("Block Size",1);
00886   polyList.set("Keep Hessenberg", true);
00887 
00888   // Create a simple status test that either reaches the relative residual tolerance or maximum polynomial size.
00889   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxItrTst =
00890     Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxDegree_ ) );
00891 
00892   // Implicit residual test, using the native residual to determine if convergence was achieved.
00893   Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTst = 
00894     Teuchos::rcp( new StatusTestGenResNorm<ScalarType,MV,OP>( polytol_ ) );
00895   convTst->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00896 
00897   // Convergence test that stops the iteration when either are satisfied.
00898   Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > polyTest = 
00899     Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, maxItrTst, convTst ) );
00900 
00901   // Create Gmres iteration object to perform one cycle of Gmres.
00902   Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
00903   gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(newProblem,printer_,polyTest,ortho_,polyList) );
00904 
00905   // Create the first block in the current Krylov basis (residual).
00906   Teuchos::RCP<MV> V_0 = MVT::Clone( *(newProblem->getRHS()), 1 );
00907   newProblem->computeCurrPrecResVec( &*V_0 );
00908 
00909   // Get a matrix to hold the orthonormalization coefficients.
00910   poly_r0_ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(1) );
00911 
00912   // Orthonormalize the new V_0
00913   int rank = ortho_->normalize( *V_0, poly_r0_ );
00914   TEST_FOR_EXCEPTION(rank != 1,GmresPolySolMgrOrthoFailure,
00915     "Belos::GmresPolySolMgr::generatePoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
00916 
00917   // Set the new state and initialize the solver.
00918   GmresIterationState<ScalarType,MV> newstate;
00919   newstate.V = V_0;
00920   newstate.z = poly_r0_;
00921   newstate.curDim = 0;
00922   gmres_iter->initializeGmres(newstate);
00923 
00924   // Perform Gmres iteration
00925   bool polyConverged = false;
00926   try {
00927     gmres_iter->iterate();
00928 
00929     // Check convergence first
00930     if ( convTst->getStatus() == Passed ) {
00931       // we have convergence
00932       polyConverged = true;
00933     }
00934   }
00935   catch (GmresIterationOrthoFailure e) {
00936     // Try to recover the most recent least-squares solution
00937     gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
00938 
00939     // Check to see if the most recent least-squares solution yielded convergence.
00940     polyTest->checkStatus( &*gmres_iter );
00941     if (convTst->getStatus() == Passed)
00942       polyConverged = true;
00943     }
00944     catch (std::exception e) {
00945     printer_->stream(Errors) << "Error! Caught exception in BlockGmresIter::iterate() at iteration "
00946                              << gmres_iter->getNumIters() << endl
00947                              << e.what() << endl;
00948     throw;
00949   }
00950 
00951   // Get the solution for this polynomial, use in comparison below
00952   Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
00953 
00954   // Record polynomial info, get current GMRES state
00955   GmresIterationState<ScalarType,MV> gmresState = gmres_iter->getState();
00956 
00957   // If the polynomial has no dimension, the tolerance is too low, return false
00958   poly_dim_ = gmresState.curDim;
00959   if (poly_dim_ == 0) {
00960     return false;
00961   }
00962   //
00963   //  Make a view and then copy the RHS of the least squares problem.
00964   //
00965   poly_y_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *gmresState.z, poly_dim_, 1 ) );
00966   poly_H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( *gmresState.H ) );
00967   //
00968   // Solve the least squares problem.
00969   //
00970   const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00971   Teuchos::BLAS<int,ScalarType> blas;
00972   blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00973              Teuchos::NON_UNIT_DIAG, poly_dim_, 1, one,
00974              gmresState.R->values(), gmresState.R->stride(), 
00975              poly_y_->values(), poly_y_->stride() );
00976   //
00977   // Generate the polynomial operator
00978   //
00979   poly_Op_ = Teuchos::rcp( 
00980                new Belos::GmresPolyOp<ScalarType,MV,OP>( problem_, poly_H_, poly_y_, poly_r0_ ) );
00981  
00982   return true;
00983 }
00984   
00985 // solve()
00986 template<class ScalarType, class MV, class OP>
00987 ReturnType GmresPolySolMgr<ScalarType,MV,OP>::solve() {
00988 
00989   // Set the current parameters if they were not set before.
00990   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00991   // then didn't set any parameters using setParameters().
00992   if (!isSet_) {
00993     setParameters(Teuchos::parameterList(*getValidParameters()));
00994   }
00995 
00996   Teuchos::BLAS<int,ScalarType> blas;
00997   Teuchos::LAPACK<int,ScalarType> lapack;
00998   
00999   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,GmresPolySolMgrLinearProblemFailure,
01000     "Belos::GmresPolySolMgr::solve(): Linear problem is not a valid object.");
01001 
01002   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),GmresPolySolMgrLinearProblemFailure,
01003     "Belos::GmresPolySolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
01004 
01005   if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
01006     TEST_FOR_EXCEPTION( checkStatusTest(),GmresPolySolMgrLinearProblemFailure,
01007       "Belos::GmresPolySolMgr::solve(): Linear problem and requested status tests are incompatible.");
01008   }
01009 
01010   // If the GMRES polynomial has not been constructed for this matrix, preconditioner pair, generate it
01011   if (!isPolyBuilt_) {
01012     Teuchos::TimeMonitor slvtimer(*timerPoly_);
01013     isPolyBuilt_ = generatePoly();
01014     TEST_FOR_EXCEPTION( !isPolyBuilt_ && poly_dim_==0, GmresPolySolMgrPolynomialFailure,
01015       "Belos::GmresPolySolMgr::generatePoly(): Failed to generate a non-trivial polynomial, reduce polynomial tolerance.");
01016     TEST_FOR_EXCEPTION( !isPolyBuilt_, GmresPolySolMgrPolynomialFailure,
01017       "Belos::GmresPolySolMgr::generatePoly(): Failed to generate polynomial that satisfied requirements.");
01018   } 
01019 
01020   // Assume convergence is achieved if user does not require strict convergence.
01021   bool isConverged = true;    
01022 
01023   // Solve the linear system using the polynomial
01024   {
01025     Teuchos::TimeMonitor slvtimer(*timerSolve_);
01026     
01027     // Apply the polynomial to the current linear system
01028     poly_Op_->Apply( *problem_->getRHS(), *problem_->getLHS() );
01029     
01030     // Reset the problem to acknowledge the updated solution
01031     problem_->setProblem(); 
01032     
01033     // If we have to strictly adhere to the requested convergence tolerance, set up a standard GMRES solver.
01034     if (strictConvTol_) {
01035       
01036       // Create indices for the linear systems to be solved.
01037       int startPtr = 0;
01038       int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
01039       int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01040       
01041       std::vector<int> currIdx;
01042       //  If an adaptive block size is allowed then only the linear systems that need to be solved are solved.
01043       //  Otherwise, the index set is generated that informs the linear problem that some linear systems are augmented.
01044       currIdx.resize( blockSize_ );
01045       for (int i=0; i<numCurrRHS; ++i) 
01046   { currIdx[i] = startPtr+i; }
01047       for (int i=numCurrRHS; i<blockSize_; ++i) 
01048   { currIdx[i] = -1; }
01049       
01050       // Inform the linear problem of the current linear system to solve.
01051       problem_->setLSIndex( currIdx );
01052       
01054       // Parameter list
01055       Teuchos::ParameterList plist;
01056       plist.set("Block Size",blockSize_);
01057       
01058       int dim = MVT::GetVecLength( *(problem_->getRHS()) );  
01059       if (blockSize_*numBlocks_ > dim) {
01060   int tmpNumBlocks = 0;
01061   if (blockSize_ == 1)
01062     tmpNumBlocks = dim / blockSize_;  // Allow for a good breakdown.
01063   else
01064     tmpNumBlocks = ( dim - blockSize_) / blockSize_;  // Allow for restarting.
01065   printer_->stream(Warnings) << 
01066     "Warning! Requested Krylov subspace dimension is larger that operator dimension!" << std::endl <<
01067     " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
01068   plist.set("Num Blocks",tmpNumBlocks);
01069       } 
01070       else 
01071   plist.set("Num Blocks",numBlocks_);
01072       
01073       // Reset the status test.  
01074       outputTest_->reset();
01075       loaDetected_ = false;
01076       
01077       // Assume convergence is achieved, then let any failed convergence set this to false.
01078       isConverged = true; 
01079       
01081       // BlockGmres solver
01082       
01083       Teuchos::RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter;
01084       
01085       block_gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
01086       
01087       // Enter solve() iterations
01088       while ( numRHS2Solve > 0 ) {
01089   
01090   // Set the current number of blocks and blocksize with the Gmres iteration.
01091   if (blockSize_*numBlocks_ > dim) {
01092     int tmpNumBlocks = 0;
01093     if (blockSize_ == 1)
01094       tmpNumBlocks = dim / blockSize_;  // Allow for a good breakdown.
01095     else
01096       tmpNumBlocks = ( dim - blockSize_) / blockSize_;  // Allow for restarting.
01097     block_gmres_iter->setSize( blockSize_, tmpNumBlocks );
01098   }
01099   else
01100     block_gmres_iter->setSize( blockSize_, numBlocks_ );
01101   
01102   // Reset the number of iterations.
01103   block_gmres_iter->resetNumIters();
01104   
01105   // Reset the number of calls that the status test output knows about.
01106   outputTest_->resetNumCalls();
01107   
01108   // Create the first block in the current Krylov basis.
01109   Teuchos::RCP<MV> V_0;
01110   V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
01111   
01112   
01113   // Get a matrix to hold the orthonormalization coefficients.
01114   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_0 = 
01115     rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
01116   
01117   // Orthonormalize the new V_0
01118   int rank = ortho_->normalize( *V_0, z_0 );
01119   TEST_FOR_EXCEPTION(rank != blockSize_,GmresPolySolMgrOrthoFailure,
01120          "Belos::GmresPolySolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
01121   
01122   // Set the new state and initialize the solver.
01123   GmresIterationState<ScalarType,MV> newstate;
01124   newstate.V = V_0;
01125   newstate.z = z_0;
01126   newstate.curDim = 0;
01127   block_gmres_iter->initializeGmres(newstate);
01128   int numRestarts = 0;
01129   
01130   while(1) {
01131     // tell block_gmres_iter to iterate
01132     try {
01133       block_gmres_iter->iterate();
01134       
01136       //
01137       // check convergence first
01138       //
01140       if ( convTest_->getStatus() == Passed ) {
01141         if ( expConvTest_->getLOADetected() ) {
01142               // we don't have convergence
01143     loaDetected_ = true;
01144     isConverged = false;
01145         }
01146         break;  // break from while(1){block_gmres_iter->iterate()}
01147       }
01149       //
01150       // check for maximum iterations
01151       //
01153       else if ( maxIterTest_->getStatus() == Passed ) {
01154         // we don't have convergence
01155         isConverged = false;
01156         break;  // break from while(1){block_gmres_iter->iterate()}
01157       }
01159       //
01160       // check for restarting, i.e. the subspace is full
01161       //
01163       else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
01164         
01165         if ( numRestarts >= maxRestarts_ ) {
01166     isConverged = false;
01167     break; // break from while(1){block_gmres_iter->iterate()}
01168         }
01169         numRestarts++;
01170         
01171         printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
01172         
01173         // Update the linear problem.
01174         Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01175         problem_->updateSolution( update, true );
01176         
01177         // Get the state.
01178         GmresIterationState<ScalarType,MV> oldState = block_gmres_iter->getState();
01179         
01180         // Compute the restart std::vector.
01181         // Get a view of the current Krylov basis.
01182         Teuchos::RCP<MV> V_0  = MVT::Clone( *(oldState.V), blockSize_ );
01183         problem_->computeCurrPrecResVec( &*V_0 );
01184         
01185         // Get a view of the first block of the Krylov basis.
01186         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_0 = 
01187     rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
01188         
01189         // Orthonormalize the new V_0
01190         int rank = ortho_->normalize( *V_0, z_0 );
01191         TEST_FOR_EXCEPTION(rank != blockSize_,GmresPolySolMgrOrthoFailure,
01192          "Belos::GmresPolySolMgr::solve(): Failed to compute initial block of orthonormal vectors after restart.");
01193         
01194         // Set the new state and initialize the solver.
01195         GmresIterationState<ScalarType,MV> newstate;
01196         newstate.V = V_0;
01197         newstate.z = z_0;
01198         newstate.curDim = 0;
01199         block_gmres_iter->initializeGmres(newstate);
01200         
01201       } // end of restarting
01202       
01204       //
01205       // we returned from iterate(), but none of our status tests Passed.
01206       // something is wrong, and it is probably our fault.
01207       //
01209       
01210       else {
01211         TEST_FOR_EXCEPTION(true,std::logic_error,
01212          "Belos::GmresPolySolMgr::solve(): Invalid return from BlockGmresIter::iterate().");
01213       }
01214     }
01215     catch (const GmresIterationOrthoFailure &e) {
01216       // If the block size is not one, it's not considered a lucky breakdown.
01217       if (blockSize_ != 1) {
01218         printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 
01219                << block_gmres_iter->getNumIters() << std::endl 
01220                << e.what() << std::endl;
01221         if (convTest_->getStatus() != Passed)
01222     isConverged = false;
01223         break;
01224       } 
01225       else {
01226         // If the block size is one, try to recover the most recent least-squares solution
01227         block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
01228         
01229         // Check to see if the most recent least-squares solution yielded convergence.
01230         sTest_->checkStatus( &*block_gmres_iter );
01231         if (convTest_->getStatus() != Passed)
01232     isConverged = false;
01233         break;
01234       }
01235     }
01236     catch (const std::exception &e) {
01237       printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 
01238              << block_gmres_iter->getNumIters() << std::endl 
01239              << e.what() << std::endl;
01240       throw;
01241     }
01242   }
01243   
01244   // Compute the current solution.
01245   // Update the linear problem.
01246   // Attempt to get the current solution from the residual status test, if it has one.
01247   if ( !Teuchos::is_null(expConvTest_->getSolution()) ) {
01248     Teuchos::RCP<MV> newX = expConvTest_->getSolution();
01249     Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
01250     MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
01251   }
01252   else {
01253     Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01254     problem_->updateSolution( update, true );
01255   }  
01256   
01257   // Inform the linear problem that we are finished with this block linear system.
01258   problem_->setCurrLS();
01259   
01260   // Update indices for the linear systems to be solved.
01261   startPtr += numCurrRHS;
01262   numRHS2Solve -= numCurrRHS;
01263   if ( numRHS2Solve > 0 ) {
01264     numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01265     
01266     currIdx.resize( blockSize_ );
01267     for (int i=0; i<numCurrRHS; ++i) 
01268       { currIdx[i] = startPtr+i; }
01269     for (int i=numCurrRHS; i<blockSize_; ++i) 
01270       { currIdx[i] = -1; }
01271     
01272     // Set the next indices.
01273     problem_->setLSIndex( currIdx );
01274   }
01275   else {
01276     currIdx.resize( numRHS2Solve );
01277   }
01278   
01279       }// while ( numRHS2Solve > 0 )
01280       
01281       // print final summary
01282       sTest_->print( printer_->stream(FinalSummary) );
01283       
01284     } // if (strictConvTol_)
01285   } // timing block
01286 
01287   // print timing information
01288   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01289   
01290   if (!isConverged || loaDetected_) {
01291     return Unconverged; // return from GmresPolySolMgr::solve() 
01292   }
01293   return Converged; // return from GmresPolySolMgr::solve() 
01294 }
01295   
01296 //  This method requires the solver manager to return a std::string that describes itself.
01297 template<class ScalarType, class MV, class OP>
01298 std::string GmresPolySolMgr<ScalarType,MV,OP>::description() const
01299   {
01300   std::ostringstream oss;
01301   oss << "Belos::GmresPolySolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01302   oss << "{";
01303   oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" <<blockSize_;
01304   oss << ", Num Blocks=" <<numBlocks_<< ", Max Restarts=" << maxRestarts_;
01305   oss << "}";
01306   return oss.str();
01307 }
01308   
01309 } // end Belos namespace
01310 
01311 #endif /* BELOS_GMRES_POLY_SOLMGR_HPP */

Generated on Wed May 12 21:30:08 2010 for Belos by  doxygen 1.4.7