Belos Version of the Day
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 the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef BELOS_GMRES_POLY_SOLMGR_HPP
00043 #define BELOS_GMRES_POLY_SOLMGR_HPP
00044 
00048 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 
00052 #include "BelosLinearProblem.hpp"
00053 #include "BelosSolverManager.hpp"
00054 
00055 #include "BelosGmresPolyOp.hpp"
00056 #include "BelosGmresIteration.hpp"
00057 #include "BelosBlockGmresIter.hpp"
00058 #include "BelosDGKSOrthoManager.hpp"
00059 #include "BelosICGSOrthoManager.hpp"
00060 #include "BelosIMGSOrthoManager.hpp"
00061 #include "BelosStatusTestMaxIters.hpp"
00062 #include "BelosStatusTestGenResNorm.hpp"
00063 #include "BelosStatusTestImpResNorm.hpp"
00064 #include "BelosStatusTestCombo.hpp"
00065 #include "BelosStatusTestOutputFactory.hpp"
00066 #include "BelosOutputManager.hpp"
00067 #include "Teuchos_BLAS.hpp"
00068 #include "Teuchos_LAPACK.hpp"
00069 #include "Teuchos_as.hpp"
00070 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00071 #include "Teuchos_TimeMonitor.hpp"
00072 #endif
00073 
00074 
00075 namespace Belos {
00076 
00078 
00079 
00086 class GmresPolySolMgrLinearProblemFailure : public BelosError {public:
00087   GmresPolySolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00088     {}};
00089 
00096 class GmresPolySolMgrPolynomialFailure : public BelosError {public:
00097   GmresPolySolMgrPolynomialFailure(const std::string& what_arg) : BelosError(what_arg)
00098     {}};
00099 
00106 class GmresPolySolMgrOrthoFailure : public BelosError {public:
00107   GmresPolySolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00108     {}};
00109 
00155 template<class ScalarType, class MV, class OP>
00156 class GmresPolySolMgr : public SolverManager<ScalarType,MV,OP> {
00157 private:
00158   typedef MultiVecTraits<ScalarType,MV> MVT;
00159   typedef MultiVecTraitsExt<ScalarType,MV> MVText;
00160   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00161   typedef Teuchos::ScalarTraits<ScalarType> STS;
00162   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00163   typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00164 
00165 public:
00166 
00168 
00169 
00175   GmresPolySolMgr();
00176 
00190   GmresPolySolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00191     const Teuchos::RCP<Teuchos::ParameterList> &pl );
00192 
00194   virtual ~GmresPolySolMgr() {};
00196 
00198 
00199 
00202   const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00203     return *problem_;
00204   }
00205 
00208   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00209 
00212   Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00213 
00219   Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00220     return Teuchos::tuple(timerSolve_, timerPoly_);
00221   }
00222 
00224   int getNumIters() const {
00225     return numIters_;
00226   }
00227 
00231   bool isLOADetected() const { return loaDetected_; }
00232 
00234 
00236 
00237 
00239   void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; isSTSet_ = false; }
00240 
00242   void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00243 
00245 
00246 
00247 
00256   void reset( const ResetType type ) {
00257     if ((type & Belos::Problem) && ! problem_.is_null ()) {
00258       problem_->setProblem ();
00259       isPolyBuilt_ = false;  // Rebuild the GMRES polynomial
00260     }
00261   }
00262 
00264 
00265 
00266 
00284   ReturnType solve();
00285 
00287 
00290 
00292   std::string description() const;
00293 
00295 
00296 private:
00297 
00298   // Method for checking current status test against defined linear problem.
00299   bool checkStatusTest();
00300 
00301   // Method for generating GMRES polynomial.
00302   bool generatePoly();
00303 
00304   // Linear problem.
00305   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00306 
00307   // Output manager.
00308   Teuchos::RCP<OutputManager<ScalarType> > printer_;
00309   Teuchos::RCP<std::ostream> outputStream_;
00310 
00311   // Status test.
00312   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00313   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00314   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00315   Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
00316   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00317 
00318   // Orthogonalization manager.
00319   Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
00320 
00321   // Current parameter list.
00322   Teuchos::RCP<Teuchos::ParameterList> params_;
00323 
00324   // Default solver values.
00325   static const MagnitudeType polytol_default_;
00326   static const MagnitudeType convtol_default_;
00327   static const MagnitudeType orthoKappa_default_;
00328   static const int maxDegree_default_;
00329   static const int maxRestarts_default_;
00330   static const int maxIters_default_;
00331   static const bool strictConvTol_default_;
00332   static const bool showMaxResNormOnly_default_;
00333   static const int blockSize_default_;
00334   static const int numBlocks_default_;
00335   static const int verbosity_default_;
00336   static const int outputStyle_default_;
00337   static const int outputFreq_default_;
00338   static const std::string impResScale_default_;
00339   static const std::string expResScale_default_;
00340   static const std::string label_default_;
00341   static const std::string orthoType_default_;
00342   static const Teuchos::RCP<std::ostream> outputStream_default_;
00343 
00344   // Current solver values.
00345   MagnitudeType polytol_, convtol_, orthoKappa_;
00346   int maxDegree_, maxRestarts_, maxIters_, numIters_;
00347   int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_;
00348   bool strictConvTol_, showMaxResNormOnly_;
00349   std::string orthoType_;
00350   std::string impResScale_, expResScale_;
00351 
00352   // Polynomial storage
00353   int poly_dim_;
00354   Teuchos::RCP<Teuchos::SerialDenseMatrix<int, ScalarType> > poly_H_, poly_y_;
00355   Teuchos::RCP<Teuchos::SerialDenseVector<int, ScalarType> > poly_r0_;
00356   Teuchos::RCP<Belos::GmresPolyOp<ScalarType, MV, OP> > poly_Op_;
00357 
00358   // Timers.
00359   std::string label_;
00360   Teuchos::RCP<Teuchos::Time> timerSolve_, timerPoly_;
00361 
00362   // Internal state variables.
00363   bool isPolyBuilt_;
00364   bool isSet_, isSTSet_, expResTest_;
00365   bool loaDetected_;
00366 
00368   mutable Teuchos::RCP<const Teuchos::ParameterList> validPL_;
00369 };
00370 
00371 
00372 // Default solver values.
00373 template<class ScalarType, class MV, class OP>
00374 const typename GmresPolySolMgr<ScalarType,MV,OP>::MagnitudeType
00375 GmresPolySolMgr<ScalarType,MV,OP>::polytol_default_ = 1e-12;
00376 
00377 template<class ScalarType, class MV, class OP>
00378 const typename GmresPolySolMgr<ScalarType,MV,OP>::MagnitudeType
00379 GmresPolySolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00380 
00381 template<class ScalarType, class MV, class OP>
00382 const typename GmresPolySolMgr<ScalarType,MV,OP>::MagnitudeType
00383 GmresPolySolMgr<ScalarType,MV,OP>::orthoKappa_default_ =
00384   -Teuchos::ScalarTraits<MagnitudeType>::one();
00385 
00386 template<class ScalarType, class MV, class OP>
00387 const int GmresPolySolMgr<ScalarType,MV,OP>::maxDegree_default_ = 25;
00388 
00389 template<class ScalarType, class MV, class OP>
00390 const int GmresPolySolMgr<ScalarType,MV,OP>::maxRestarts_default_ = 20;
00391 
00392 template<class ScalarType, class MV, class OP>
00393 const int GmresPolySolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00394 
00395 template<class ScalarType, class MV, class OP>
00396 const bool GmresPolySolMgr<ScalarType,MV,OP>::strictConvTol_default_ = false;
00397 
00398 template<class ScalarType, class MV, class OP>
00399 const bool GmresPolySolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00400 
00401 template<class ScalarType, class MV, class OP>
00402 const int GmresPolySolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00403 
00404 template<class ScalarType, class MV, class OP>
00405 const int GmresPolySolMgr<ScalarType,MV,OP>::numBlocks_default_ = 300;
00406 
00407 template<class ScalarType, class MV, class OP>
00408 const int GmresPolySolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00409 
00410 template<class ScalarType, class MV, class OP>
00411 const int GmresPolySolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00412 
00413 template<class ScalarType, class MV, class OP>
00414 const int GmresPolySolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00415 
00416 template<class ScalarType, class MV, class OP>
00417 const std::string GmresPolySolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of RHS";
00418 
00419 template<class ScalarType, class MV, class OP>
00420 const std::string GmresPolySolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of RHS";
00421 
00422 template<class ScalarType, class MV, class OP>
00423 const std::string GmresPolySolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00424 
00425 template<class ScalarType, class MV, class OP>
00426 const std::string GmresPolySolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00427 
00428 template<class ScalarType, class MV, class OP>
00429 const Teuchos::RCP<std::ostream>
00430 GmresPolySolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcpFromRef (std::cout);
00431 
00432 
00433 template<class ScalarType, class MV, class OP>
00434 GmresPolySolMgr<ScalarType,MV,OP>::GmresPolySolMgr () :
00435   outputStream_ (outputStream_default_),
00436   polytol_ (polytol_default_),
00437   convtol_ (convtol_default_),
00438   orthoKappa_ (orthoKappa_default_),
00439   maxDegree_ (maxDegree_default_),
00440   maxRestarts_ (maxRestarts_default_),
00441   maxIters_ (maxIters_default_),
00442   numIters_ (0),
00443   blockSize_ (blockSize_default_),
00444   numBlocks_ (numBlocks_default_),
00445   verbosity_ (verbosity_default_),
00446   outputStyle_ (outputStyle_default_),
00447   outputFreq_ (outputFreq_default_),
00448   strictConvTol_ (strictConvTol_default_),
00449   showMaxResNormOnly_ (showMaxResNormOnly_default_),
00450   orthoType_ (orthoType_default_),
00451   impResScale_ (impResScale_default_),
00452   expResScale_ (expResScale_default_),
00453   poly_dim_ (0),
00454   label_ (label_default_),
00455   isPolyBuilt_ (false),
00456   isSet_ (false),
00457   isSTSet_ (false),
00458   expResTest_ (false),
00459   loaDetected_ (false)
00460 {}
00461 
00462 
00463 template<class ScalarType, class MV, class OP>
00464 GmresPolySolMgr<ScalarType,MV,OP>::
00465 GmresPolySolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00466                  const Teuchos::RCP<Teuchos::ParameterList> &pl) :
00467   problem_ (problem),
00468   outputStream_ (outputStream_default_),
00469   polytol_ (polytol_default_),
00470   convtol_ (convtol_default_),
00471   orthoKappa_ (orthoKappa_default_),
00472   maxDegree_ (maxDegree_default_),
00473   maxRestarts_ (maxRestarts_default_),
00474   maxIters_ (maxIters_default_),
00475   numIters_ (0),
00476   blockSize_ (blockSize_default_),
00477   numBlocks_ (numBlocks_default_),
00478   verbosity_ (verbosity_default_),
00479   outputStyle_ (outputStyle_default_),
00480   outputFreq_ (outputFreq_default_),
00481   strictConvTol_ (strictConvTol_default_),
00482   showMaxResNormOnly_ (showMaxResNormOnly_default_),
00483   orthoType_ (orthoType_default_),
00484   impResScale_ (impResScale_default_),
00485   expResScale_ (expResScale_default_),
00486   poly_dim_ (0),
00487   label_ (label_default_),
00488   isPolyBuilt_ (false),
00489   isSet_ (false),
00490   isSTSet_ (false),
00491   expResTest_ (false),
00492   loaDetected_ (false)
00493 {
00494   TEUCHOS_TEST_FOR_EXCEPTION(
00495     problem_.is_null (), std::invalid_argument,
00496     "Belos::GmresPolySolMgr: The given linear problem is null.  "
00497     "Please call this constructor with a nonnull LinearProblem argument, "
00498     "or call the constructor that does not take a LinearProblem.");
00499 
00500   // If the input parameter list is null, then the parameters take
00501   // default values.
00502   if (! pl.is_null ()) {
00503     setParameters (pl);
00504   }
00505 }
00506 
00507 
00508 template<class ScalarType, class MV, class OP>
00509 Teuchos::RCP<const Teuchos::ParameterList>
00510 GmresPolySolMgr<ScalarType,MV,OP>::getValidParameters() const
00511 {
00512   if (validPL_.is_null ()) {
00513     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList ();
00514     pl->set("Polynomial Tolerance", polytol_default_,
00515       "The relative residual tolerance that used to construct the GMRES polynomial.");
00516     pl->set("Maximum Degree", maxDegree_default_,
00517       "The maximum degree allowed for any GMRES polynomial.");
00518     pl->set("Convergence Tolerance", convtol_default_,
00519       "The relative residual tolerance that needs to be achieved by the\n"
00520       "iterative solver in order for the linear system to be declared converged." );
00521     pl->set("Maximum Restarts", maxRestarts_default_,
00522       "The maximum number of restarts allowed for each\n"
00523       "set of RHS solved.");
00524     pl->set("Maximum Iterations", maxIters_default_,
00525       "The maximum number of block iterations allowed for each\n"
00526       "set of RHS solved.");
00527     pl->set("Num Blocks", numBlocks_default_,
00528       "The maximum number of blocks allowed in the Krylov subspace\n"
00529       "for each set of RHS solved.");
00530     pl->set("Block Size", blockSize_default_,
00531       "The number of vectors in each block.  This number times the\n"
00532       "number of blocks is the total Krylov subspace dimension.");
00533     pl->set("Verbosity", verbosity_default_,
00534       "What type(s) of solver information should be outputted\n"
00535       "to the output stream.");
00536     pl->set("Output Style", outputStyle_default_,
00537       "What style is used for the solver information outputted\n"
00538       "to the output stream.");
00539     pl->set("Output Frequency", outputFreq_default_,
00540       "How often convergence information should be outputted\n"
00541       "to the output stream.");
00542     pl->set("Output Stream", outputStream_default_,
00543       "A reference-counted pointer to the output stream where all\n"
00544       "solver output is sent.");
00545     pl->set("Strict Convergence", strictConvTol_default_,
00546       "After polynomial is applied, whether solver should try to achieve\n"
00547       "the relative residual tolerance.");
00548     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00549       "When convergence information is printed, only show the maximum\n"
00550       "relative residual norm when the block size is greater than one.");
00551     pl->set("Implicit Residual Scaling", impResScale_default_,
00552       "The type of scaling used in the implicit residual convergence test.");
00553     pl->set("Explicit Residual Scaling", expResScale_default_,
00554       "The type of scaling used in the explicit residual convergence test.");
00555     pl->set("Timer Label", label_default_,
00556       "The string to use as a prefix for the timer labels.");
00557     //  pl->set("Restart Timers", restartTimers_);
00558     pl->set("Orthogonalization", orthoType_default_,
00559       "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
00560     pl->set("Orthogonalization Constant",orthoKappa_default_,
00561       "The constant used by DGKS orthogonalization to determine\n"
00562       "whether another step of classical Gram-Schmidt is necessary.");
00563     validPL_ = pl;
00564   }
00565   return validPL_;
00566 }
00567 
00568 
00569 template<class ScalarType, class MV, class OP>
00570 void GmresPolySolMgr<ScalarType,MV,OP>::
00571 setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
00572 {
00573   // Create the internal parameter list if ones doesn't already exist.
00574   if (params_.is_null ()) {
00575     params_ = Teuchos::parameterList (*getValidParameters ());
00576   }
00577   else {
00578     params->validateParameters (*getValidParameters ());
00579   }
00580 
00581   // Check for maximum polynomial degree
00582   if (params->isParameter("Maximum Degree")) {
00583     maxDegree_ = params->get("Maximum Degree",maxDegree_default_);
00584 
00585     // Update parameter in our list.
00586     params_->set("Maximum Degree", maxDegree_);
00587   }
00588 
00589   // Check for maximum number of restarts
00590   if (params->isParameter("Maximum Restarts")) {
00591     maxRestarts_ = params->get("Maximum Restarts",maxRestarts_default_);
00592 
00593     // Update parameter in our list.
00594     params_->set("Maximum Restarts", maxRestarts_);
00595   }
00596 
00597   // Check for maximum number of iterations
00598   if (params->isParameter("Maximum Iterations")) {
00599     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00600 
00601     // Update parameter in our list and in status test.
00602     params_->set("Maximum Iterations", maxIters_);
00603     if (maxIterTest_!=Teuchos::null)
00604       maxIterTest_->setMaxIters( maxIters_ );
00605   }
00606 
00607   // Check for blocksize
00608   if (params->isParameter("Block Size")) {
00609     blockSize_ = params->get("Block Size",blockSize_default_);
00610     TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00611       "Belos::GmresPolySolMgr: \"Block Size\" must be strictly positive.");
00612 
00613     // Update parameter in our list.
00614     params_->set("Block Size", blockSize_);
00615   }
00616 
00617   // Check for the maximum number of blocks.
00618   if (params->isParameter("Num Blocks")) {
00619     numBlocks_ = params->get("Num Blocks",numBlocks_default_);
00620     TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00621       "Belos::GmresPolySolMgr: \"Num Blocks\" must be strictly positive.");
00622 
00623     // Update parameter in our list.
00624     params_->set("Num Blocks", numBlocks_);
00625   }
00626 
00627   // Check to see if the timer label changed.
00628   if (params->isParameter("Timer Label")) {
00629     std::string tempLabel = params->get("Timer Label", label_default_);
00630 
00631     // Update parameter in our list and solver timer
00632     if (tempLabel != label_) {
00633       label_ = tempLabel;
00634       params_->set("Timer Label", label_);
00635       std::string solveLabel = label_ + ": GmresPolySolMgr total solve time";
00636 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00637       timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
00638 #endif
00639       std::string polyLabel = label_ + ": GmresPolySolMgr polynomial creation time";
00640 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00641       timerPoly_ = Teuchos::TimeMonitor::getNewCounter(polyLabel);
00642 #endif
00643       if (ortho_ != Teuchos::null) {
00644         ortho_->setLabel( label_ );
00645       }
00646     }
00647   }
00648 
00649   // Check if the orthogonalization changed.
00650   if (params->isParameter("Orthogonalization")) {
00651     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00652     TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS",
00653                         std::invalid_argument,
00654                         "Belos::GmresPolySolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00655     if (tempOrthoType != orthoType_) {
00656       orthoType_ = tempOrthoType;
00657       // Create orthogonalization manager
00658       if (orthoType_=="DGKS") {
00659         if (orthoKappa_ <= 0) {
00660           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00661         }
00662         else {
00663           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00664           Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00665         }
00666       }
00667       else if (orthoType_=="ICGS") {
00668         ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00669       }
00670       else if (orthoType_=="IMGS") {
00671         ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00672       }
00673     }
00674   }
00675 
00676   // Check which orthogonalization constant to use.
00677   if (params->isParameter("Orthogonalization Constant")) {
00678     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00679 
00680     // Update parameter in our list.
00681     params_->set("Orthogonalization Constant",orthoKappa_);
00682     if (orthoType_=="DGKS") {
00683       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00684         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00685       }
00686     }
00687   }
00688 
00689   // Check for a change in verbosity level
00690   if (params->isParameter("Verbosity")) {
00691     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00692       verbosity_ = params->get("Verbosity", verbosity_default_);
00693     } else {
00694       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00695     }
00696 
00697     // Update parameter in our list.
00698     params_->set("Verbosity", verbosity_);
00699     if (printer_ != Teuchos::null)
00700       printer_->setVerbosity(verbosity_);
00701   }
00702 
00703   // Check for a change in output style
00704   if (params->isParameter("Output Style")) {
00705     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00706       outputStyle_ = params->get("Output Style", outputStyle_default_);
00707     } else {
00708       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00709     }
00710 
00711     // Reconstruct the convergence test if the explicit residual test is not being used.
00712     params_->set("Output Style", outputStyle_);
00713     if (outputTest_ != Teuchos::null) {
00714       isSTSet_ = false;
00715     }
00716   }
00717 
00718   // output stream
00719   if (params->isParameter("Output Stream")) {
00720     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00721 
00722     // Update parameter in our list.
00723     params_->set("Output Stream", outputStream_);
00724     if (printer_ != Teuchos::null)
00725       printer_->setOStream( outputStream_ );
00726   }
00727 
00728   // frequency level
00729   if (verbosity_ & Belos::StatusTestDetails) {
00730     if (params->isParameter("Output Frequency")) {
00731       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00732     }
00733 
00734     // Update parameter in out list and output status test.
00735     params_->set("Output Frequency", outputFreq_);
00736     if (outputTest_ != Teuchos::null)
00737       outputTest_->setOutputFrequency( outputFreq_ );
00738   }
00739 
00740   // Create output manager if we need to.
00741   if (printer_ == Teuchos::null) {
00742     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00743   }
00744 
00745   // Convergence
00746   //typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t; // unused
00747   //typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t; // unused
00748 
00749   // Check for polynomial convergence tolerance
00750   if (params->isParameter("Polynomial Tolerance")) {
00751     polytol_ = params->get("Polynomial Tolerance",polytol_default_);
00752 
00753     // Update parameter in our list and residual tests.
00754     params_->set("Polynomial Tolerance", polytol_);
00755   }
00756 
00757   // Check for convergence tolerance
00758   if (params->isParameter("Convergence Tolerance")) {
00759     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00760 
00761     // Update parameter in our list and residual tests.
00762     params_->set("Convergence Tolerance", convtol_);
00763     if (impConvTest_ != Teuchos::null)
00764       impConvTest_->setTolerance( convtol_ );
00765     if (expConvTest_ != Teuchos::null)
00766       expConvTest_->setTolerance( convtol_ );
00767   }
00768 
00769   // Check if user requires solver to reach convergence tolerance
00770   if (params->isParameter("Strict Convergence")) {
00771     strictConvTol_ = params->get("Strict Convergence",strictConvTol_default_);
00772 
00773     // Update parameter in our list and residual tests
00774     params_->set("Strict Convergence", strictConvTol_);
00775   }
00776 
00777   // Check for a change in scaling, if so we need to build new residual tests.
00778   if (params->isParameter("Implicit Residual Scaling")) {
00779     std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
00780 
00781     // Only update the scaling if it's different.
00782     if (impResScale_ != tempImpResScale) {
00783       Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
00784       impResScale_ = tempImpResScale;
00785 
00786       // Update parameter in our list and residual tests
00787       params_->set("Implicit Residual Scaling", impResScale_);
00788       if (impConvTest_ != Teuchos::null) {
00789         try {
00790           impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
00791         }
00792         catch (std::exception& e) {
00793           // Make sure the convergence test gets constructed again.
00794           isSTSet_ = false;
00795         }
00796       }
00797     }
00798   }
00799 
00800   if (params->isParameter("Explicit Residual Scaling")) {
00801     std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
00802 
00803     // Only update the scaling if it's different.
00804     if (expResScale_ != tempExpResScale) {
00805       Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
00806       expResScale_ = tempExpResScale;
00807 
00808       // Update parameter in our list and residual tests
00809       params_->set("Explicit Residual Scaling", expResScale_);
00810       if (expConvTest_ != Teuchos::null) {
00811         try {
00812           expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
00813         }
00814         catch (std::exception& e) {
00815           // Make sure the convergence test gets constructed again.
00816           isSTSet_ = false;
00817         }
00818       }
00819     }
00820   }
00821 
00822 
00823   if (params->isParameter("Show Maximum Residual Norm Only")) {
00824     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00825 
00826     // Update parameter in our list and residual tests
00827     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00828     if (impConvTest_ != Teuchos::null)
00829       impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00830     if (expConvTest_ != Teuchos::null)
00831       expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00832   }
00833 
00834   // Create orthogonalization manager if we need to.
00835   if (ortho_ == Teuchos::null) {
00836     if (orthoType_=="DGKS") {
00837       if (orthoKappa_ <= 0) {
00838         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00839       }
00840       else {
00841         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00842         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00843       }
00844     }
00845     else if (orthoType_=="ICGS") {
00846       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00847     }
00848     else if (orthoType_=="IMGS") {
00849       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00850     }
00851     else {
00852       TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00853         "Belos::GmresPolySolMgr(): Invalid orthogonalization type.");
00854     }
00855   }
00856 
00857   // Create the timers if we need to.
00858   if (timerSolve_ == Teuchos::null) {
00859     std::string solveLabel = label_ + ": GmresPolySolMgr total solve time";
00860 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00861     timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
00862 #endif
00863   }
00864 
00865   if (timerPoly_ == Teuchos::null) {
00866     std::string polyLabel = label_ + ": GmresPolySolMgr polynomial creation time";
00867 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00868     timerPoly_ = Teuchos::TimeMonitor::getNewCounter(polyLabel);
00869 #endif
00870   }
00871 
00872   // Inform the solver manager that the current parameters were set.
00873   isSet_ = true;
00874 }
00875 
00876 // Check the status test versus the defined linear problem
00877 template<class ScalarType, class MV, class OP>
00878 bool GmresPolySolMgr<ScalarType,MV,OP>::checkStatusTest() {
00879 
00880   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00881   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestGenResNorm_t;
00882   typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP>  StatusTestImpResNorm_t;
00883 
00884   // Basic test checks maximum iterations and native residual.
00885   maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00886 
00887   // If there is a left preconditioner, we create a combined status test that checks the implicit
00888   // and then explicit residual norm to see if we have convergence.
00889   if (!Teuchos::is_null(problem_->getLeftPrec())) {
00890     expResTest_ = true;
00891   }
00892 
00893   if (expResTest_) {
00894 
00895     // Implicit residual test, using the native residual to determine if convergence was achieved.
00896     Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
00897       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
00898     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00899     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00900     impConvTest_ = tmpImpConvTest;
00901 
00902     // Explicit residual test once the native residual is below the tolerance
00903     Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
00904       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
00905     tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
00906     tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
00907     tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00908     expConvTest_ = tmpExpConvTest;
00909 
00910     // The convergence test is a combination of the "cheap" implicit test and explicit test.
00911     convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
00912   }
00913   else {
00914 
00915     // Implicit residual test, using the native residual to determine if convergence was achieved.
00916     // Use test that checks for loss of accuracy.
00917     Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
00918       Teuchos::rcp( new StatusTestImpResNorm_t( convtol_ ) );
00919     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00920     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00921     impConvTest_ = tmpImpConvTest;
00922 
00923     // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
00924     expConvTest_ = impConvTest_;
00925     convTest_ = impConvTest_;
00926   }
00927 
00928   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00929 
00930   // Create the status test output class.
00931   // This class manages and formats the output from the status test.
00932   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00933   outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00934 
00935   // Set the solver string for the output test
00936   std::string solverDesc = " Gmres Polynomial ";
00937   outputTest_->setSolverDesc( solverDesc );
00938 
00939 
00940   // The status test is now set.
00941   isSTSet_ = true;
00942 
00943   return false;
00944 }
00945 
00946 template<class ScalarType, class MV, class OP>
00947 bool GmresPolySolMgr<ScalarType,MV,OP>::generatePoly()
00948 {
00949   // Create a copy of the linear problem that has a zero initial guess and random RHS.
00950   Teuchos::RCP<MV> newX  = MVT::Clone( *(problem_->getLHS()), 1 );
00951   Teuchos::RCP<MV> newB  = MVT::Clone( *(problem_->getRHS()), 1 );
00952   MVT::MvInit( *newX, STS::zero() );
00953   MVT::MvRandom( *newB );
00954   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > newProblem =
00955     Teuchos::rcp( new LinearProblem<ScalarType,MV,OP>( problem_->getOperator(), newX, newB ) );
00956   newProblem->setLeftPrec( problem_->getLeftPrec() );
00957   newProblem->setRightPrec( problem_->getRightPrec() );
00958   newProblem->setLabel("Belos GMRES Poly Generation");
00959   newProblem->setProblem();
00960   std::vector<int> idx(1,0);       // Must set the index to be the first vector (0)!
00961   newProblem->setLSIndex( idx );
00962 
00963   // Create a parameter list for the GMRES iteration.
00964   Teuchos::ParameterList polyList;
00965 
00966   // Tell the block solver that the block size is one.
00967   polyList.set("Num Blocks",maxDegree_);
00968   polyList.set("Block Size",1);
00969   polyList.set("Keep Hessenberg", true);
00970 
00971   // Create a simple status test that either reaches the relative residual tolerance or maximum polynomial size.
00972   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxItrTst =
00973     Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxDegree_ ) );
00974 
00975   // Implicit residual test, using the native residual to determine if convergence was achieved.
00976   Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTst =
00977     Teuchos::rcp( new StatusTestGenResNorm<ScalarType,MV,OP>( polytol_ ) );
00978   convTst->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00979 
00980   // Convergence test that stops the iteration when either are satisfied.
00981   Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > polyTest =
00982     Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, maxItrTst, convTst ) );
00983 
00984   // Create Gmres iteration object to perform one cycle of Gmres.
00985   Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
00986   gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(newProblem,printer_,polyTest,ortho_,polyList) );
00987 
00988   // Create the first block in the current Krylov basis (residual).
00989   Teuchos::RCP<MV> V_0 = MVT::Clone( *(newProblem->getRHS()), 1 );
00990   newProblem->computeCurrPrecResVec( &*V_0 );
00991 
00992   // Get a matrix to hold the orthonormalization coefficients.
00993   poly_r0_ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(1) );
00994 
00995   // Orthonormalize the new V_0
00996   int rank = ortho_->normalize( *V_0, poly_r0_ );
00997   TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GmresPolySolMgrOrthoFailure,
00998     "Belos::GmresPolySolMgr::generatePoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
00999 
01000   // Set the new state and initialize the solver.
01001   GmresIterationState<ScalarType,MV> newstate;
01002   newstate.V = V_0;
01003   newstate.z = poly_r0_;
01004   newstate.curDim = 0;
01005   gmres_iter->initializeGmres(newstate);
01006 
01007   // Perform Gmres iteration
01008   bool polyConverged = false;
01009   try {
01010     gmres_iter->iterate();
01011 
01012     // Check convergence first
01013     if ( convTst->getStatus() == Passed ) {
01014       // we have convergence
01015       polyConverged = true;
01016     }
01017   }
01018   catch (GmresIterationOrthoFailure e) {
01019     // Try to recover the most recent least-squares solution
01020     gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
01021 
01022     // Check to see if the most recent least-squares solution yielded convergence.
01023     polyTest->checkStatus( &*gmres_iter );
01024     if (convTst->getStatus() == Passed)
01025       polyConverged = true;
01026     }
01027     catch (std::exception e) {
01028     printer_->stream(Errors) << "Error! Caught exception in BlockGmresIter::iterate() at iteration "
01029                              << gmres_iter->getNumIters() << std::endl
01030                              << e.what() << std::endl;
01031     throw;
01032   }
01033 
01034   // FIXME (mfh 27 Apr 2013) Why aren't we using polyConverged after
01035   // setting it?  I'm tired of the compile warning so I'll silence it
01036   // here, but I'm curious why we aren't using the variable.
01037   (void) polyConverged;
01038 
01039   // Get the solution for this polynomial, use in comparison below
01040   Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
01041 
01042   // Record polynomial info, get current GMRES state
01043   GmresIterationState<ScalarType,MV> gmresState = gmres_iter->getState();
01044 
01045   // If the polynomial has no dimension, the tolerance is too low, return false
01046   poly_dim_ = gmresState.curDim;
01047   if (poly_dim_ == 0) {
01048     return false;
01049   }
01050   //
01051   //  Make a view and then copy the RHS of the least squares problem.
01052   //
01053   poly_y_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *gmresState.z, poly_dim_, 1 ) );
01054   poly_H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( *gmresState.H ) );
01055   //
01056   // Solve the least squares problem.
01057   //
01058   const ScalarType one = STS::one ();
01059   Teuchos::BLAS<int,ScalarType> blas;
01060   blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
01061              Teuchos::NON_UNIT_DIAG, poly_dim_, 1, one,
01062              gmresState.R->values(), gmresState.R->stride(),
01063              poly_y_->values(), poly_y_->stride() );
01064   //
01065   // Generate the polynomial operator
01066   //
01067   poly_Op_ = Teuchos::rcp(
01068                new Belos::GmresPolyOp<ScalarType,MV,OP>( problem_, poly_H_, poly_y_, poly_r0_ ) );
01069 
01070   return true;
01071 }
01072 
01073 
01074 template<class ScalarType, class MV, class OP>
01075 ReturnType GmresPolySolMgr<ScalarType,MV,OP>::solve ()
01076 {
01077   using Teuchos::RCP;
01078   using Teuchos::rcp;
01079   typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
01080 
01081   // Set the current parameters if they were not set before.  NOTE:
01082   // This may occur if the user generated the solver manager with the
01083   // default constructor and then didn't set any parameters using
01084   // setParameters().
01085   if (! isSet_) {
01086     setParameters (Teuchos::parameterList (*getValidParameters ()));
01087   }
01088 
01089   Teuchos::BLAS<int,ScalarType> blas;
01090   Teuchos::LAPACK<int,ScalarType> lapack;
01091 
01092   TEUCHOS_TEST_FOR_EXCEPTION(
01093     problem_.is_null (), GmresPolySolMgrLinearProblemFailure,
01094     "Belos::GmresPolySolMgr::solve: The linear problem has not been set yet, "
01095     "or was set to null.  Please call setProblem() with a nonnull input before "
01096     "calling solve().");
01097 
01098   TEUCHOS_TEST_FOR_EXCEPTION(
01099     ! problem_->isProblemSet (), GmresPolySolMgrLinearProblemFailure,
01100     "Belos::GmresPolySolMgr::solve: The linear problem is not ready.  Please "
01101     "call setProblem() on the LinearProblem object before calling solve().");
01102 
01103   if (! isSTSet_ || (! expResTest_ && ! problem_->getLeftPrec ().is_null ())) {
01104     // checkStatusTest() shouldn't have side effects, but it's still
01105     // not such a good idea to put possibly side-effect-y function
01106     // calls in a macro invocation.  It could be expensive if the
01107     // macro evaluates the argument more than once, for example.
01108     const bool check = checkStatusTest ();
01109     TEUCHOS_TEST_FOR_EXCEPTION(
01110       check, GmresPolySolMgrLinearProblemFailure,
01111       "Belos::GmresPolySolMgr::solve: Linear problem and requested status "
01112       "tests are incompatible.");
01113   }
01114 
01115   // If the GMRES polynomial has not been constructed for this
01116   // (nmatrix, preconditioner) pair, generate it.
01117   if (! isPolyBuilt_) {
01118 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01119     Teuchos::TimeMonitor slvtimer(*timerPoly_);
01120 #endif
01121     isPolyBuilt_ = generatePoly();
01122     TEUCHOS_TEST_FOR_EXCEPTION( !isPolyBuilt_ && poly_dim_==0, GmresPolySolMgrPolynomialFailure,
01123       "Belos::GmresPolySolMgr::generatePoly(): Failed to generate a non-trivial polynomial, reduce polynomial tolerance.");
01124     TEUCHOS_TEST_FOR_EXCEPTION( !isPolyBuilt_, GmresPolySolMgrPolynomialFailure,
01125       "Belos::GmresPolySolMgr::generatePoly(): Failed to generate polynomial that satisfied requirements.");
01126   }
01127 
01128   // Assume convergence is achieved if user does not require strict convergence.
01129   bool isConverged = true;
01130 
01131   // Solve the linear system using the polynomial
01132   {
01133 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01134     Teuchos::TimeMonitor slvtimer(*timerSolve_);
01135 #endif
01136 
01137     // Apply the polynomial to the current linear system
01138     poly_Op_->Apply( *problem_->getRHS(), *problem_->getLHS() );
01139 
01140     // Reset the problem to acknowledge the updated solution
01141     problem_->setProblem ();
01142 
01143     // If we have to strictly adhere to the requested convergence tolerance, set up a standard GMRES solver.
01144     if (strictConvTol_) {
01145 
01146       // Create indices for the linear systems to be solved.
01147       int startPtr = 0;
01148       int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
01149       int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01150 
01151 
01152       // If an adaptive block size is allowed then only the linear
01153       // systems that need to be solved are solved.  Otherwise, the
01154       // index set is generated that informs the linear problem that
01155       // some linear systems are augmented.
01156 
01157       std::vector<int> currIdx (blockSize_);
01158       for (int i = 0; i < numCurrRHS; ++i) {
01159         currIdx[i] = startPtr+i;
01160       }
01161       for (int i = numCurrRHS; i < blockSize_; ++i) {
01162         currIdx[i] = -1;
01163       }
01164 
01165       // Inform the linear problem of the current linear system to solve.
01166       problem_->setLSIndex (currIdx);
01167 
01169       // Parameter list
01170       Teuchos::ParameterList plist;
01171       plist.set ("Block Size", blockSize_);
01172 
01173       ptrdiff_t dim = MVText::GetGlobalLength( *(problem_->getRHS()) );
01174       if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
01175         ptrdiff_t tmpNumBlocks = 0;
01176         if (blockSize_ == 1) {
01177           tmpNumBlocks = dim / blockSize_;  // Allow for a good breakdown.
01178         } else {
01179           tmpNumBlocks = ( dim - blockSize_) / blockSize_;  // Allow for restarting.
01180         }
01181         printer_->stream(Warnings)
01182           << "Warning! Requested Krylov subspace dimension is larger than "
01183           << "operator dimension!" << std::endl << "The maximum number of "
01184           << "blocks allowed for the Krylov subspace will be adjusted to "
01185           << tmpNumBlocks << std::endl;
01186         plist.set ("Num Blocks", Teuchos::asSafe<int> (tmpNumBlocks));
01187       }
01188       else {
01189         plist.set ("Num Blocks", numBlocks_);
01190       }
01191 
01192       // Reset the status test.
01193       outputTest_->reset ();
01194       loaDetected_ = false;
01195 
01196       // Assume convergence is achieved, then let any failed
01197       // convergence set this to false.
01198       isConverged = true;
01199 
01200       //
01201       // Solve using BlockGmres
01202       //
01203 
01204       RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter =
01205         rcp (new BlockGmresIter<ScalarType,MV,OP> (problem_, printer_,
01206                                                    outputTest_, ortho_, plist));
01207 
01208       // Enter solve() iterations
01209       while (numRHS2Solve > 0) {
01210         // Set the current number of blocks and block size with the
01211         // Gmres iteration.
01212         if (blockSize_*numBlocks_ > dim) {
01213           int tmpNumBlocks = 0;
01214           if (blockSize_ == 1) {
01215             // Leave space for happy breakdown.
01216             tmpNumBlocks = dim / blockSize_;
01217           } else {
01218             // Leave space for restarting.
01219             tmpNumBlocks = (dim - blockSize_) / blockSize_;
01220           }
01221           block_gmres_iter->setSize (blockSize_, tmpNumBlocks);
01222         }
01223         else {
01224           block_gmres_iter->setSize (blockSize_, numBlocks_);
01225         }
01226 
01227         // Reset the number of iterations.
01228         block_gmres_iter->resetNumIters ();
01229 
01230         // Reset the number of calls that the status test output knows about.
01231         outputTest_->resetNumCalls ();
01232 
01233         // Create the first block in the current Krylov basis.
01234         RCP<MV> V_0 = MVT::CloneCopy (*(problem_->getInitPrecResVec ()), currIdx);
01235 
01236         // Get a matrix to hold the orthonormalization coefficients.
01237         RCP<SDM> z_0 = rcp (new SDM (blockSize_, blockSize_));
01238 
01239         // Orthonormalize the new V_0
01240         int rank = ortho_->normalize (*V_0, z_0);
01241         TEUCHOS_TEST_FOR_EXCEPTION(
01242           rank != blockSize_, GmresPolySolMgrOrthoFailure,
01243           "Belos::GmresPolySolMgr::solve: Failed to compute initial block of "
01244           "orthonormal vectors.");
01245 
01246         // Set the new state and initialize the solver.
01247         GmresIterationState<ScalarType,MV> newstate;
01248         newstate.V = V_0;
01249         newstate.z = z_0;
01250         newstate.curDim = 0;
01251         block_gmres_iter->initializeGmres(newstate);
01252         int numRestarts = 0;
01253 
01254         while(1) {
01255           try {
01256             block_gmres_iter->iterate();
01257 
01258             // check convergence first
01259             if ( convTest_->getStatus() == Passed ) {
01260               if ( expConvTest_->getLOADetected() ) {
01261               // we don't have convergence
01262                 loaDetected_ = true;
01263                 isConverged = false;
01264               }
01265               break;  // break from while(1){block_gmres_iter->iterate()}
01266             }
01267 
01268             // check for maximum iterations
01269             else if ( maxIterTest_->getStatus() == Passed ) {
01270               // we don't have convergence
01271               isConverged = false;
01272               break;  // break from while(1){block_gmres_iter->iterate()}
01273             }
01274             // check for restarting, i.e. the subspace is full
01275             else if (block_gmres_iter->getCurSubspaceDim () ==
01276                      block_gmres_iter->getMaxSubspaceDim ()) {
01277               if (numRestarts >= maxRestarts_) {
01278                 isConverged = false;
01279                 break; // break from while(1){block_gmres_iter->iterate()}
01280               }
01281               numRestarts++;
01282 
01283               printer_->stream(Debug)
01284                 << " Performing restart number " << numRestarts << " of "
01285                 << maxRestarts_ << std::endl << std::endl;
01286 
01287               // Update the linear problem.
01288               RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01289               problem_->updateSolution (update, true);
01290 
01291               // Get the state.
01292               GmresIterationState<ScalarType,MV> oldState = block_gmres_iter->getState();
01293 
01294               // Compute the restart vector.
01295               // Get a view of the current Krylov basis.
01296               //
01297               // We call this VV_0 to avoid shadowing the previously
01298               // defined V_0 above.
01299               RCP<MV> VV_0  = MVT::Clone (*(oldState.V), blockSize_);
01300               problem_->computeCurrPrecResVec (&*VV_0);
01301 
01302               // Get a view of the first block of the Krylov basis.
01303               //
01304               // We call this zz_0 to avoid shadowing the previously
01305               // defined z_0 above.
01306               RCP<SDM> zz_0 = rcp (new SDM (blockSize_, blockSize_));
01307 
01308               // Orthonormalize the new VV_0
01309               const int theRank = ortho_->normalize( *VV_0, zz_0 );
01310               TEUCHOS_TEST_FOR_EXCEPTION(
01311                 theRank != blockSize_, GmresPolySolMgrOrthoFailure,
01312                 "Belos::GmresPolySolMgr::solve: Failed to compute initial "
01313                 "block of orthonormal vectors after restart.");
01314 
01315               // Set the new state and initialize the solver.
01316               GmresIterationState<ScalarType,MV> theNewState;
01317               theNewState.V = VV_0;
01318               theNewState.z = zz_0;
01319               theNewState.curDim = 0;
01320               block_gmres_iter->initializeGmres (theNewState);
01321             } // end of restarting
01322             //
01323             // We returned from iterate(), but none of our status
01324             // tests Passed.  Something is wrong, and it is probably
01325             // our fault.
01326             //
01327             else {
01328               TEUCHOS_TEST_FOR_EXCEPTION(
01329                 true, std::logic_error,
01330                 "Belos::GmresPolySolMgr::solve: Invalid return from "
01331                 "BlockGmresIter::iterate().  Please report this bug "
01332                 "to the Belos developers.");
01333             }
01334           }
01335           catch (const GmresIterationOrthoFailure& e) {
01336             // If the block size is not one, it's not considered a lucky breakdown.
01337             if (blockSize_ != 1) {
01338               printer_->stream(Errors)
01339                 << "Error! Caught std::exception in BlockGmresIter::iterate() "
01340                 << "at iteration " << block_gmres_iter->getNumIters()
01341                 << std::endl << e.what() << std::endl;
01342               if (convTest_->getStatus() != Passed) {
01343                 isConverged = false;
01344               }
01345               break;
01346             }
01347             else {
01348               // If the block size is one, try to recover the most
01349               // recent least-squares solution
01350               block_gmres_iter->updateLSQR (block_gmres_iter->getCurSubspaceDim ());
01351 
01352               // Check to see if the most recent least-squares
01353               // solution yielded convergence.
01354               sTest_->checkStatus (&*block_gmres_iter);
01355               if (convTest_->getStatus() != Passed) {
01356                 isConverged = false;
01357               }
01358               break;
01359             }
01360           }
01361           catch (const std::exception &e) {
01362             printer_->stream(Errors)
01363               << "Error! Caught std::exception in BlockGmresIter::iterate() "
01364               << "at iteration " << block_gmres_iter->getNumIters() << std::endl
01365               << e.what() << std::endl;
01366             throw;
01367           }
01368         }
01369 
01370         // Compute the current solution.  Update the linear problem.
01371         // Attempt to get the current solution from the residual
01372         // status test, if it has one.
01373         if (! Teuchos::is_null (expConvTest_->getSolution ()) ) {
01374           RCP<MV> newX = expConvTest_->getSolution ();
01375           RCP<MV> curX = problem_->getCurrLHSVec ();
01376           MVT::MvAddMv (STS::zero (), *newX, STS::one (), *newX, *curX);
01377         }
01378         else {
01379           RCP<MV> update = block_gmres_iter->getCurrentUpdate ();
01380           problem_->updateSolution (update, true);
01381         }
01382 
01383         // Inform the linear problem that we are finished with this block linear system.
01384         problem_->setCurrLS ();
01385 
01386         // Update indices for the linear systems to be solved.
01387         startPtr += numCurrRHS;
01388         numRHS2Solve -= numCurrRHS;
01389         if (numRHS2Solve > 0) {
01390           numCurrRHS = (numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01391 
01392           currIdx.resize (blockSize_);
01393           for (int i=0; i<numCurrRHS; ++i) {
01394             currIdx[i] = startPtr+i;
01395           }
01396           for (int i=numCurrRHS; i<blockSize_; ++i) {
01397             currIdx[i] = -1;
01398           }
01399 
01400           // Set the next indices.
01401           problem_->setLSIndex( currIdx );
01402         }
01403         else {
01404           currIdx.resize( numRHS2Solve );
01405         }
01406 
01407       }// while ( numRHS2Solve > 0 )
01408 
01409       // print final summary
01410       sTest_->print( printer_->stream(FinalSummary) );
01411 
01412     } // if (strictConvTol_)
01413   } // timing block
01414 
01415   // print timing information
01416 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01417   // Calling summarize() can be expensive, so don't call unless the
01418   // user wants to print out timing details.  summarize() will do all
01419   // the work even if it's passed a "black hole" output stream.
01420   if (verbosity_ & TimingDetails)
01421     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01422 #endif
01423 
01424   if (!isConverged || loaDetected_) {
01425     return Unconverged; // return from GmresPolySolMgr::solve()
01426   }
01427   return Converged; // return from GmresPolySolMgr::solve()
01428 }
01429 
01430 
01431 template<class ScalarType, class MV, class OP>
01432 std::string GmresPolySolMgr<ScalarType,MV,OP>::description () const
01433 {
01434   std::ostringstream out;
01435 
01436   out << "\"Belos::GmresPolySolMgr\": {"
01437       << "ScalarType: " << Teuchos::TypeNameTraits<ScalarType>::name ()
01438       << ", Ortho Type: " << orthoType_
01439       << ", Block Size: " << blockSize_
01440       << ", Num Blocks: " << numBlocks_
01441       << ", Max Restarts: " << maxRestarts_;
01442   out << "}";
01443   return out.str ();
01444 }
01445 
01446 } // namespace Belos
01447 
01448 #endif // BELOS_GMRES_POLY_SOLMGR_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines