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   blockSize_ (blockSize_default_),
00443   numBlocks_ (numBlocks_default_),
00444   verbosity_ (verbosity_default_),
00445   outputStyle_ (outputStyle_default_),
00446   outputFreq_ (outputFreq_default_),
00447   strictConvTol_ (strictConvTol_default_),
00448   showMaxResNormOnly_ (showMaxResNormOnly_default_),
00449   orthoType_ (orthoType_default_),
00450   impResScale_ (impResScale_default_),
00451   expResScale_ (expResScale_default_),
00452   label_ (label_default_),
00453   isPolyBuilt_ (false),
00454   isSet_ (false),
00455   isSTSet_ (false),
00456   expResTest_ (false),
00457   loaDetected_ (false)
00458 {}
00459 
00460 
00461 template<class ScalarType, class MV, class OP>
00462 GmresPolySolMgr<ScalarType,MV,OP>::
00463 GmresPolySolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00464                  const Teuchos::RCP<Teuchos::ParameterList> &pl) :
00465   problem_ (problem),
00466   outputStream_ (outputStream_default_),
00467   polytol_ (polytol_default_),
00468   convtol_ (convtol_default_),
00469   orthoKappa_ (orthoKappa_default_),
00470   maxDegree_ (maxDegree_default_),
00471   maxRestarts_ (maxRestarts_default_),
00472   maxIters_ (maxIters_default_),
00473   blockSize_ (blockSize_default_),
00474   numBlocks_ (numBlocks_default_),
00475   verbosity_ (verbosity_default_),
00476   outputStyle_ (outputStyle_default_),
00477   outputFreq_ (outputFreq_default_),
00478   strictConvTol_ (strictConvTol_default_),
00479   showMaxResNormOnly_ (showMaxResNormOnly_default_),
00480   orthoType_ (orthoType_default_),
00481   impResScale_ (impResScale_default_),
00482   expResScale_ (expResScale_default_),
00483   label_ (label_default_),
00484   isPolyBuilt_ (false),
00485   isSet_ (false),
00486   isSTSet_ (false),
00487   expResTest_ (false),
00488   loaDetected_ (false)
00489 {
00490   TEUCHOS_TEST_FOR_EXCEPTION(
00491     problem_.is_null (), std::invalid_argument,
00492     "Belos::GmresPolySolMgr: The given linear problem is null.  "
00493     "Please call this constructor with a nonnull LinearProblem argument, "
00494     "or call the constructor that does not take a LinearProblem.");
00495 
00496   // If the input parameter list is null, then the parameters take
00497   // default values.
00498   if (! pl.is_null ()) {
00499     setParameters (pl);
00500   }
00501 }
00502 
00503 
00504 template<class ScalarType, class MV, class OP>
00505 Teuchos::RCP<const Teuchos::ParameterList>
00506 GmresPolySolMgr<ScalarType,MV,OP>::getValidParameters() const
00507 {
00508   if (validPL_.is_null ()) {
00509     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList ();
00510     pl->set("Polynomial Tolerance", polytol_default_,
00511       "The relative residual tolerance that used to construct the GMRES polynomial.");
00512     pl->set("Maximum Degree", maxDegree_default_,
00513       "The maximum degree allowed for any GMRES polynomial.");
00514     pl->set("Convergence Tolerance", convtol_default_,
00515       "The relative residual tolerance that needs to be achieved by the\n"
00516       "iterative solver in order for the linear system to be declared converged." );
00517     pl->set("Maximum Restarts", maxRestarts_default_,
00518       "The maximum number of restarts allowed for each\n"
00519       "set of RHS solved.");
00520     pl->set("Maximum Iterations", maxIters_default_,
00521       "The maximum number of block iterations allowed for each\n"
00522       "set of RHS solved.");
00523     pl->set("Num Blocks", numBlocks_default_,
00524       "The maximum number of blocks allowed in the Krylov subspace\n"
00525       "for each set of RHS solved.");
00526     pl->set("Block Size", blockSize_default_,
00527       "The number of vectors in each block.  This number times the\n"
00528       "number of blocks is the total Krylov subspace dimension.");
00529     pl->set("Verbosity", verbosity_default_,
00530       "What type(s) of solver information should be outputted\n"
00531       "to the output stream.");
00532     pl->set("Output Style", outputStyle_default_,
00533       "What style is used for the solver information outputted\n"
00534       "to the output stream.");
00535     pl->set("Output Frequency", outputFreq_default_,
00536       "How often convergence information should be outputted\n"
00537       "to the output stream.");
00538     pl->set("Output Stream", outputStream_default_,
00539       "A reference-counted pointer to the output stream where all\n"
00540       "solver output is sent.");
00541     pl->set("Strict Convergence", strictConvTol_default_,
00542       "After polynomial is applied, whether solver should try to achieve\n"
00543       "the relative residual tolerance.");
00544     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00545       "When convergence information is printed, only show the maximum\n"
00546       "relative residual norm when the block size is greater than one.");
00547     pl->set("Implicit Residual Scaling", impResScale_default_,
00548       "The type of scaling used in the implicit residual convergence test.");
00549     pl->set("Explicit Residual Scaling", expResScale_default_,
00550       "The type of scaling used in the explicit residual convergence test.");
00551     pl->set("Timer Label", label_default_,
00552       "The string to use as a prefix for the timer labels.");
00553     //  pl->set("Restart Timers", restartTimers_);
00554     pl->set("Orthogonalization", orthoType_default_,
00555       "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
00556     pl->set("Orthogonalization Constant",orthoKappa_default_,
00557       "The constant used by DGKS orthogonalization to determine\n"
00558       "whether another step of classical Gram-Schmidt is necessary.");
00559     validPL_ = pl;
00560   }
00561   return validPL_;
00562 }
00563 
00564 
00565 template<class ScalarType, class MV, class OP>
00566 void GmresPolySolMgr<ScalarType,MV,OP>::
00567 setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
00568 {
00569   // Create the internal parameter list if ones doesn't already exist.
00570   if (params_.is_null ()) {
00571     params_ = Teuchos::parameterList (*getValidParameters ());
00572   }
00573   else {
00574     params->validateParameters (*getValidParameters ());
00575   }
00576 
00577   // Check for maximum polynomial degree
00578   if (params->isParameter("Maximum Degree")) {
00579     maxDegree_ = params->get("Maximum Degree",maxDegree_default_);
00580 
00581     // Update parameter in our list.
00582     params_->set("Maximum Degree", maxDegree_);
00583   }
00584 
00585   // Check for maximum number of restarts
00586   if (params->isParameter("Maximum Restarts")) {
00587     maxRestarts_ = params->get("Maximum Restarts",maxRestarts_default_);
00588 
00589     // Update parameter in our list.
00590     params_->set("Maximum Restarts", maxRestarts_);
00591   }
00592 
00593   // Check for maximum number of iterations
00594   if (params->isParameter("Maximum Iterations")) {
00595     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00596 
00597     // Update parameter in our list and in status test.
00598     params_->set("Maximum Iterations", maxIters_);
00599     if (maxIterTest_!=Teuchos::null)
00600       maxIterTest_->setMaxIters( maxIters_ );
00601   }
00602 
00603   // Check for blocksize
00604   if (params->isParameter("Block Size")) {
00605     blockSize_ = params->get("Block Size",blockSize_default_);
00606     TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00607       "Belos::GmresPolySolMgr: \"Block Size\" must be strictly positive.");
00608 
00609     // Update parameter in our list.
00610     params_->set("Block Size", blockSize_);
00611   }
00612 
00613   // Check for the maximum number of blocks.
00614   if (params->isParameter("Num Blocks")) {
00615     numBlocks_ = params->get("Num Blocks",numBlocks_default_);
00616     TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00617       "Belos::GmresPolySolMgr: \"Num Blocks\" must be strictly positive.");
00618 
00619     // Update parameter in our list.
00620     params_->set("Num Blocks", numBlocks_);
00621   }
00622 
00623   // Check to see if the timer label changed.
00624   if (params->isParameter("Timer Label")) {
00625     std::string tempLabel = params->get("Timer Label", label_default_);
00626 
00627     // Update parameter in our list and solver timer
00628     if (tempLabel != label_) {
00629       label_ = tempLabel;
00630       params_->set("Timer Label", label_);
00631       std::string solveLabel = label_ + ": GmresPolySolMgr total solve time";
00632 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00633       timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
00634 #endif
00635       std::string polyLabel = label_ + ": GmresPolySolMgr polynomial creation time";
00636 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00637       timerPoly_ = Teuchos::TimeMonitor::getNewCounter(polyLabel);
00638 #endif
00639       if (ortho_ != Teuchos::null) {
00640         ortho_->setLabel( label_ );
00641       }
00642     }
00643   }
00644 
00645   // Check if the orthogonalization changed.
00646   if (params->isParameter("Orthogonalization")) {
00647     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00648     TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS",
00649                         std::invalid_argument,
00650                         "Belos::GmresPolySolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00651     if (tempOrthoType != orthoType_) {
00652       orthoType_ = tempOrthoType;
00653       // Create orthogonalization manager
00654       if (orthoType_=="DGKS") {
00655         if (orthoKappa_ <= 0) {
00656           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00657         }
00658         else {
00659           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00660           Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00661         }
00662       }
00663       else if (orthoType_=="ICGS") {
00664         ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00665       }
00666       else if (orthoType_=="IMGS") {
00667         ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00668       }
00669     }
00670   }
00671 
00672   // Check which orthogonalization constant to use.
00673   if (params->isParameter("Orthogonalization Constant")) {
00674     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00675 
00676     // Update parameter in our list.
00677     params_->set("Orthogonalization Constant",orthoKappa_);
00678     if (orthoType_=="DGKS") {
00679       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00680         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00681       }
00682     }
00683   }
00684 
00685   // Check for a change in verbosity level
00686   if (params->isParameter("Verbosity")) {
00687     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00688       verbosity_ = params->get("Verbosity", verbosity_default_);
00689     } else {
00690       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00691     }
00692 
00693     // Update parameter in our list.
00694     params_->set("Verbosity", verbosity_);
00695     if (printer_ != Teuchos::null)
00696       printer_->setVerbosity(verbosity_);
00697   }
00698 
00699   // Check for a change in output style
00700   if (params->isParameter("Output Style")) {
00701     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00702       outputStyle_ = params->get("Output Style", outputStyle_default_);
00703     } else {
00704       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00705     }
00706 
00707     // Reconstruct the convergence test if the explicit residual test is not being used.
00708     params_->set("Output Style", outputStyle_);
00709     if (outputTest_ != Teuchos::null) {
00710       isSTSet_ = false;
00711     }
00712   }
00713 
00714   // output stream
00715   if (params->isParameter("Output Stream")) {
00716     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00717 
00718     // Update parameter in our list.
00719     params_->set("Output Stream", outputStream_);
00720     if (printer_ != Teuchos::null)
00721       printer_->setOStream( outputStream_ );
00722   }
00723 
00724   // frequency level
00725   if (verbosity_ & Belos::StatusTestDetails) {
00726     if (params->isParameter("Output Frequency")) {
00727       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00728     }
00729 
00730     // Update parameter in out list and output status test.
00731     params_->set("Output Frequency", outputFreq_);
00732     if (outputTest_ != Teuchos::null)
00733       outputTest_->setOutputFrequency( outputFreq_ );
00734   }
00735 
00736   // Create output manager if we need to.
00737   if (printer_ == Teuchos::null) {
00738     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00739   }
00740 
00741   // Convergence
00742   //typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t; // unused
00743   //typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t; // unused
00744 
00745   // Check for polynomial convergence tolerance
00746   if (params->isParameter("Polynomial Tolerance")) {
00747     polytol_ = params->get("Polynomial Tolerance",polytol_default_);
00748 
00749     // Update parameter in our list and residual tests.
00750     params_->set("Polynomial Tolerance", polytol_);
00751   }
00752 
00753   // Check for convergence tolerance
00754   if (params->isParameter("Convergence Tolerance")) {
00755     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00756 
00757     // Update parameter in our list and residual tests.
00758     params_->set("Convergence Tolerance", convtol_);
00759     if (impConvTest_ != Teuchos::null)
00760       impConvTest_->setTolerance( convtol_ );
00761     if (expConvTest_ != Teuchos::null)
00762       expConvTest_->setTolerance( convtol_ );
00763   }
00764 
00765   // Check if user requires solver to reach convergence tolerance
00766   if (params->isParameter("Strict Convergence")) {
00767     strictConvTol_ = params->get("Strict Convergence",strictConvTol_default_);
00768 
00769     // Update parameter in our list and residual tests
00770     params_->set("Strict Convergence", strictConvTol_);
00771   }
00772 
00773   // Check for a change in scaling, if so we need to build new residual tests.
00774   if (params->isParameter("Implicit Residual Scaling")) {
00775     std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
00776 
00777     // Only update the scaling if it's different.
00778     if (impResScale_ != tempImpResScale) {
00779       Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
00780       impResScale_ = tempImpResScale;
00781 
00782       // Update parameter in our list and residual tests
00783       params_->set("Implicit Residual Scaling", impResScale_);
00784       if (impConvTest_ != Teuchos::null) {
00785         try {
00786           impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
00787         }
00788         catch (std::exception& e) {
00789           // Make sure the convergence test gets constructed again.
00790           isSTSet_ = false;
00791         }
00792       }
00793     }
00794   }
00795 
00796   if (params->isParameter("Explicit Residual Scaling")) {
00797     std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
00798 
00799     // Only update the scaling if it's different.
00800     if (expResScale_ != tempExpResScale) {
00801       Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
00802       expResScale_ = tempExpResScale;
00803 
00804       // Update parameter in our list and residual tests
00805       params_->set("Explicit Residual Scaling", expResScale_);
00806       if (expConvTest_ != Teuchos::null) {
00807         try {
00808           expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
00809         }
00810         catch (std::exception& e) {
00811           // Make sure the convergence test gets constructed again.
00812           isSTSet_ = false;
00813         }
00814       }
00815     }
00816   }
00817 
00818 
00819   if (params->isParameter("Show Maximum Residual Norm Only")) {
00820     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00821 
00822     // Update parameter in our list and residual tests
00823     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00824     if (impConvTest_ != Teuchos::null)
00825       impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00826     if (expConvTest_ != Teuchos::null)
00827       expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00828   }
00829 
00830   // Create orthogonalization manager if we need to.
00831   if (ortho_ == Teuchos::null) {
00832     if (orthoType_=="DGKS") {
00833       if (orthoKappa_ <= 0) {
00834         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00835       }
00836       else {
00837         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00838         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00839       }
00840     }
00841     else if (orthoType_=="ICGS") {
00842       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00843     }
00844     else if (orthoType_=="IMGS") {
00845       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00846     }
00847     else {
00848       TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00849         "Belos::GmresPolySolMgr(): Invalid orthogonalization type.");
00850     }
00851   }
00852 
00853   // Create the timers if we need to.
00854   if (timerSolve_ == Teuchos::null) {
00855     std::string solveLabel = label_ + ": GmresPolySolMgr total solve time";
00856 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00857     timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
00858 #endif
00859   }
00860 
00861   if (timerPoly_ == Teuchos::null) {
00862     std::string polyLabel = label_ + ": GmresPolySolMgr polynomial creation time";
00863 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00864     timerPoly_ = Teuchos::TimeMonitor::getNewCounter(polyLabel);
00865 #endif
00866   }
00867 
00868   // Inform the solver manager that the current parameters were set.
00869   isSet_ = true;
00870 }
00871 
00872 // Check the status test versus the defined linear problem
00873 template<class ScalarType, class MV, class OP>
00874 bool GmresPolySolMgr<ScalarType,MV,OP>::checkStatusTest() {
00875 
00876   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00877   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestGenResNorm_t;
00878   typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP>  StatusTestImpResNorm_t;
00879 
00880   // Basic test checks maximum iterations and native residual.
00881   maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00882 
00883   // If there is a left preconditioner, we create a combined status test that checks the implicit
00884   // and then explicit residual norm to see if we have convergence.
00885   if (!Teuchos::is_null(problem_->getLeftPrec())) {
00886     expResTest_ = true;
00887   }
00888 
00889   if (expResTest_) {
00890 
00891     // Implicit residual test, using the native residual to determine if convergence was achieved.
00892     Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
00893       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
00894     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00895     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00896     impConvTest_ = tmpImpConvTest;
00897 
00898     // Explicit residual test once the native residual is below the tolerance
00899     Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
00900       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
00901     tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
00902     tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
00903     tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00904     expConvTest_ = tmpExpConvTest;
00905 
00906     // The convergence test is a combination of the "cheap" implicit test and explicit test.
00907     convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
00908   }
00909   else {
00910 
00911     // Implicit residual test, using the native residual to determine if convergence was achieved.
00912     // Use test that checks for loss of accuracy.
00913     Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
00914       Teuchos::rcp( new StatusTestImpResNorm_t( convtol_ ) );
00915     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00916     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00917     impConvTest_ = tmpImpConvTest;
00918 
00919     // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
00920     expConvTest_ = impConvTest_;
00921     convTest_ = impConvTest_;
00922   }
00923 
00924   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00925 
00926   // Create the status test output class.
00927   // This class manages and formats the output from the status test.
00928   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00929   outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00930 
00931   // Set the solver string for the output test
00932   std::string solverDesc = " Gmres Polynomial ";
00933   outputTest_->setSolverDesc( solverDesc );
00934 
00935 
00936   // The status test is now set.
00937   isSTSet_ = true;
00938 
00939   return false;
00940 }
00941 
00942 template<class ScalarType, class MV, class OP>
00943 bool GmresPolySolMgr<ScalarType,MV,OP>::generatePoly()
00944 {
00945   // Create a copy of the linear problem that has a zero initial guess and random RHS.
00946   Teuchos::RCP<MV> newX  = MVT::Clone( *(problem_->getLHS()), 1 );
00947   Teuchos::RCP<MV> newB  = MVT::Clone( *(problem_->getRHS()), 1 );
00948   MVT::MvInit( *newX, STS::zero() );
00949   MVT::MvRandom( *newB );
00950   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > newProblem =
00951     Teuchos::rcp( new LinearProblem<ScalarType,MV,OP>( problem_->getOperator(), newX, newB ) );
00952   newProblem->setLeftPrec( problem_->getLeftPrec() );
00953   newProblem->setRightPrec( problem_->getRightPrec() );
00954   newProblem->setLabel("Belos GMRES Poly Generation");
00955   newProblem->setProblem();
00956   std::vector<int> idx(1,0);       // Must set the index to be the first vector (0)!
00957   newProblem->setLSIndex( idx );
00958 
00959   // Create a parameter list for the GMRES iteration.
00960   Teuchos::ParameterList polyList;
00961 
00962   // Tell the block solver that the block size is one.
00963   polyList.set("Num Blocks",maxDegree_);
00964   polyList.set("Block Size",1);
00965   polyList.set("Keep Hessenberg", true);
00966 
00967   // Create a simple status test that either reaches the relative residual tolerance or maximum polynomial size.
00968   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxItrTst =
00969     Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxDegree_ ) );
00970 
00971   // Implicit residual test, using the native residual to determine if convergence was achieved.
00972   Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTst =
00973     Teuchos::rcp( new StatusTestGenResNorm<ScalarType,MV,OP>( polytol_ ) );
00974   convTst->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00975 
00976   // Convergence test that stops the iteration when either are satisfied.
00977   Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > polyTest =
00978     Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, maxItrTst, convTst ) );
00979 
00980   // Create Gmres iteration object to perform one cycle of Gmres.
00981   Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
00982   gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(newProblem,printer_,polyTest,ortho_,polyList) );
00983 
00984   // Create the first block in the current Krylov basis (residual).
00985   Teuchos::RCP<MV> V_0 = MVT::Clone( *(newProblem->getRHS()), 1 );
00986   newProblem->computeCurrPrecResVec( &*V_0 );
00987 
00988   // Get a matrix to hold the orthonormalization coefficients.
00989   poly_r0_ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(1) );
00990 
00991   // Orthonormalize the new V_0
00992   int rank = ortho_->normalize( *V_0, poly_r0_ );
00993   TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GmresPolySolMgrOrthoFailure,
00994     "Belos::GmresPolySolMgr::generatePoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
00995 
00996   // Set the new state and initialize the solver.
00997   GmresIterationState<ScalarType,MV> newstate;
00998   newstate.V = V_0;
00999   newstate.z = poly_r0_;
01000   newstate.curDim = 0;
01001   gmres_iter->initializeGmres(newstate);
01002 
01003   // Perform Gmres iteration
01004   bool polyConverged = false;
01005   try {
01006     gmres_iter->iterate();
01007 
01008     // Check convergence first
01009     if ( convTst->getStatus() == Passed ) {
01010       // we have convergence
01011       polyConverged = true;
01012     }
01013   }
01014   catch (GmresIterationOrthoFailure e) {
01015     // Try to recover the most recent least-squares solution
01016     gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
01017 
01018     // Check to see if the most recent least-squares solution yielded convergence.
01019     polyTest->checkStatus( &*gmres_iter );
01020     if (convTst->getStatus() == Passed)
01021       polyConverged = true;
01022     }
01023     catch (std::exception e) {
01024     printer_->stream(Errors) << "Error! Caught exception in BlockGmresIter::iterate() at iteration "
01025                              << gmres_iter->getNumIters() << std::endl
01026                              << e.what() << std::endl;
01027     throw;
01028   }
01029 
01030   // FIXME (mfh 27 Apr 2013) Why aren't we using polyConverged after
01031   // setting it?  I'm tired of the compile warning so I'll silence it
01032   // here, but I'm curious why we aren't using the variable.
01033   (void) polyConverged;
01034 
01035   // Get the solution for this polynomial, use in comparison below
01036   Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
01037 
01038   // Record polynomial info, get current GMRES state
01039   GmresIterationState<ScalarType,MV> gmresState = gmres_iter->getState();
01040 
01041   // If the polynomial has no dimension, the tolerance is too low, return false
01042   poly_dim_ = gmresState.curDim;
01043   if (poly_dim_ == 0) {
01044     return false;
01045   }
01046   //
01047   //  Make a view and then copy the RHS of the least squares problem.
01048   //
01049   poly_y_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *gmresState.z, poly_dim_, 1 ) );
01050   poly_H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( *gmresState.H ) );
01051   //
01052   // Solve the least squares problem.
01053   //
01054   const ScalarType one = STS::one ();
01055   Teuchos::BLAS<int,ScalarType> blas;
01056   blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
01057              Teuchos::NON_UNIT_DIAG, poly_dim_, 1, one,
01058              gmresState.R->values(), gmresState.R->stride(),
01059              poly_y_->values(), poly_y_->stride() );
01060   //
01061   // Generate the polynomial operator
01062   //
01063   poly_Op_ = Teuchos::rcp(
01064                new Belos::GmresPolyOp<ScalarType,MV,OP>( problem_, poly_H_, poly_y_, poly_r0_ ) );
01065 
01066   return true;
01067 }
01068 
01069 
01070 template<class ScalarType, class MV, class OP>
01071 ReturnType GmresPolySolMgr<ScalarType,MV,OP>::solve ()
01072 {
01073   using Teuchos::RCP;
01074   using Teuchos::rcp;
01075   typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
01076 
01077   // Set the current parameters if they were not set before.  NOTE:
01078   // This may occur if the user generated the solver manager with the
01079   // default constructor and then didn't set any parameters using
01080   // setParameters().
01081   if (! isSet_) {
01082     setParameters (Teuchos::parameterList (*getValidParameters ()));
01083   }
01084 
01085   Teuchos::BLAS<int,ScalarType> blas;
01086   Teuchos::LAPACK<int,ScalarType> lapack;
01087 
01088   TEUCHOS_TEST_FOR_EXCEPTION(
01089     problem_.is_null (), GmresPolySolMgrLinearProblemFailure,
01090     "Belos::GmresPolySolMgr::solve: The linear problem has not been set yet, "
01091     "or was set to null.  Please call setProblem() with a nonnull input before "
01092     "calling solve().");
01093 
01094   TEUCHOS_TEST_FOR_EXCEPTION(
01095     ! problem_->isProblemSet (), GmresPolySolMgrLinearProblemFailure,
01096     "Belos::GmresPolySolMgr::solve: The linear problem is not ready.  Please "
01097     "call setProblem() on the LinearProblem object before calling solve().");
01098 
01099   if (! isSTSet_ || (! expResTest_ && ! problem_->getLeftPrec ().is_null ())) {
01100     // checkStatusTest() shouldn't have side effects, but it's still
01101     // not such a good idea to put possibly side-effect-y function
01102     // calls in a macro invocation.  It could be expensive if the
01103     // macro evaluates the argument more than once, for example.
01104     const bool check = checkStatusTest ();
01105     TEUCHOS_TEST_FOR_EXCEPTION(
01106       check, GmresPolySolMgrLinearProblemFailure,
01107       "Belos::GmresPolySolMgr::solve: Linear problem and requested status "
01108       "tests are incompatible.");
01109   }
01110 
01111   // If the GMRES polynomial has not been constructed for this
01112   // (nmatrix, preconditioner) pair, generate it.
01113   if (! isPolyBuilt_) {
01114 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01115     Teuchos::TimeMonitor slvtimer(*timerPoly_);
01116 #endif
01117     isPolyBuilt_ = generatePoly();
01118     TEUCHOS_TEST_FOR_EXCEPTION( !isPolyBuilt_ && poly_dim_==0, GmresPolySolMgrPolynomialFailure,
01119       "Belos::GmresPolySolMgr::generatePoly(): Failed to generate a non-trivial polynomial, reduce polynomial tolerance.");
01120     TEUCHOS_TEST_FOR_EXCEPTION( !isPolyBuilt_, GmresPolySolMgrPolynomialFailure,
01121       "Belos::GmresPolySolMgr::generatePoly(): Failed to generate polynomial that satisfied requirements.");
01122   }
01123 
01124   // Assume convergence is achieved if user does not require strict convergence.
01125   bool isConverged = true;
01126 
01127   // Solve the linear system using the polynomial
01128   {
01129 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01130     Teuchos::TimeMonitor slvtimer(*timerSolve_);
01131 #endif
01132 
01133     // Apply the polynomial to the current linear system
01134     poly_Op_->Apply( *problem_->getRHS(), *problem_->getLHS() );
01135 
01136     // Reset the problem to acknowledge the updated solution
01137     problem_->setProblem ();
01138 
01139     // If we have to strictly adhere to the requested convergence tolerance, set up a standard GMRES solver.
01140     if (strictConvTol_) {
01141 
01142       // Create indices for the linear systems to be solved.
01143       int startPtr = 0;
01144       int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
01145       int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01146 
01147 
01148       // If an adaptive block size is allowed then only the linear
01149       // systems that need to be solved are solved.  Otherwise, the
01150       // index set is generated that informs the linear problem that
01151       // some linear systems are augmented.
01152 
01153       std::vector<int> currIdx (blockSize_);
01154       for (int i = 0; i < numCurrRHS; ++i) {
01155         currIdx[i] = startPtr+i;
01156       }
01157       for (int i = numCurrRHS; i < blockSize_; ++i) {
01158         currIdx[i] = -1;
01159       }
01160 
01161       // Inform the linear problem of the current linear system to solve.
01162       problem_->setLSIndex (currIdx);
01163 
01165       // Parameter list
01166       Teuchos::ParameterList plist;
01167       plist.set ("Block Size", blockSize_);
01168 
01169       ptrdiff_t dim = MVText::GetGlobalLength( *(problem_->getRHS()) );
01170       if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
01171         ptrdiff_t tmpNumBlocks = 0;
01172         if (blockSize_ == 1) {
01173           tmpNumBlocks = dim / blockSize_;  // Allow for a good breakdown.
01174         } else {
01175           tmpNumBlocks = ( dim - blockSize_) / blockSize_;  // Allow for restarting.
01176         }
01177         printer_->stream(Warnings)
01178           << "Warning! Requested Krylov subspace dimension is larger than "
01179           << "operator dimension!" << std::endl << "The maximum number of "
01180           << "blocks allowed for the Krylov subspace will be adjusted to "
01181           << tmpNumBlocks << std::endl;
01182         plist.set ("Num Blocks", Teuchos::asSafe<int> (tmpNumBlocks));
01183       }
01184       else {
01185         plist.set ("Num Blocks", numBlocks_);
01186       }
01187 
01188       // Reset the status test.
01189       outputTest_->reset ();
01190       loaDetected_ = false;
01191 
01192       // Assume convergence is achieved, then let any failed
01193       // convergence set this to false.
01194       isConverged = true;
01195 
01196       //
01197       // Solve using BlockGmres
01198       //
01199 
01200       RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter =
01201         rcp (new BlockGmresIter<ScalarType,MV,OP> (problem_, printer_,
01202                                                    outputTest_, ortho_, plist));
01203 
01204       // Enter solve() iterations
01205       while (numRHS2Solve > 0) {
01206         // Set the current number of blocks and block size with the
01207         // Gmres iteration.
01208         if (blockSize_*numBlocks_ > dim) {
01209           int tmpNumBlocks = 0;
01210           if (blockSize_ == 1) {
01211             // Leave space for happy breakdown.
01212             tmpNumBlocks = dim / blockSize_;
01213           } else {
01214             // Leave space for restarting.
01215             tmpNumBlocks = (dim - blockSize_) / blockSize_;
01216           }
01217           block_gmres_iter->setSize (blockSize_, tmpNumBlocks);
01218         }
01219         else {
01220           block_gmres_iter->setSize (blockSize_, numBlocks_);
01221         }
01222 
01223         // Reset the number of iterations.
01224         block_gmres_iter->resetNumIters ();
01225 
01226         // Reset the number of calls that the status test output knows about.
01227         outputTest_->resetNumCalls ();
01228 
01229         // Create the first block in the current Krylov basis.
01230         RCP<MV> V_0 = MVT::CloneCopy (*(problem_->getInitPrecResVec ()), currIdx);
01231 
01232         // Get a matrix to hold the orthonormalization coefficients.
01233         RCP<SDM> z_0 = rcp (new SDM (blockSize_, blockSize_));
01234 
01235         // Orthonormalize the new V_0
01236         int rank = ortho_->normalize (*V_0, z_0);
01237         TEUCHOS_TEST_FOR_EXCEPTION(
01238           rank != blockSize_, GmresPolySolMgrOrthoFailure,
01239           "Belos::GmresPolySolMgr::solve: Failed to compute initial block of "
01240           "orthonormal vectors.");
01241 
01242         // Set the new state and initialize the solver.
01243         GmresIterationState<ScalarType,MV> newstate;
01244         newstate.V = V_0;
01245         newstate.z = z_0;
01246         newstate.curDim = 0;
01247         block_gmres_iter->initializeGmres(newstate);
01248         int numRestarts = 0;
01249 
01250         while(1) {
01251           try {
01252             block_gmres_iter->iterate();
01253 
01254             // check convergence first
01255             if ( convTest_->getStatus() == Passed ) {
01256               if ( expConvTest_->getLOADetected() ) {
01257               // we don't have convergence
01258                 loaDetected_ = true;
01259                 isConverged = false;
01260               }
01261               break;  // break from while(1){block_gmres_iter->iterate()}
01262             }
01263 
01264             // check for maximum iterations
01265             else if ( maxIterTest_->getStatus() == Passed ) {
01266               // we don't have convergence
01267               isConverged = false;
01268               break;  // break from while(1){block_gmres_iter->iterate()}
01269             }
01270             // check for restarting, i.e. the subspace is full
01271             else if (block_gmres_iter->getCurSubspaceDim () ==
01272                      block_gmres_iter->getMaxSubspaceDim ()) {
01273               if (numRestarts >= maxRestarts_) {
01274                 isConverged = false;
01275                 break; // break from while(1){block_gmres_iter->iterate()}
01276               }
01277               numRestarts++;
01278 
01279               printer_->stream(Debug)
01280                 << " Performing restart number " << numRestarts << " of "
01281                 << maxRestarts_ << std::endl << std::endl;
01282 
01283               // Update the linear problem.
01284               RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01285               problem_->updateSolution (update, true);
01286 
01287               // Get the state.
01288               GmresIterationState<ScalarType,MV> oldState = block_gmres_iter->getState();
01289 
01290               // Compute the restart vector.
01291               // Get a view of the current Krylov basis.
01292               //
01293               // We call this VV_0 to avoid shadowing the previously
01294               // defined V_0 above.
01295               RCP<MV> VV_0  = MVT::Clone (*(oldState.V), blockSize_);
01296               problem_->computeCurrPrecResVec (&*VV_0);
01297 
01298               // Get a view of the first block of the Krylov basis.
01299               //
01300               // We call this zz_0 to avoid shadowing the previously
01301               // defined z_0 above.
01302               RCP<SDM> zz_0 = rcp (new SDM (blockSize_, blockSize_));
01303 
01304               // Orthonormalize the new VV_0
01305               const int theRank = ortho_->normalize( *VV_0, zz_0 );
01306               TEUCHOS_TEST_FOR_EXCEPTION(
01307                 theRank != blockSize_, GmresPolySolMgrOrthoFailure,
01308                 "Belos::GmresPolySolMgr::solve: Failed to compute initial "
01309                 "block of orthonormal vectors after restart.");
01310 
01311               // Set the new state and initialize the solver.
01312               GmresIterationState<ScalarType,MV> theNewState;
01313               theNewState.V = VV_0;
01314               theNewState.z = zz_0;
01315               theNewState.curDim = 0;
01316               block_gmres_iter->initializeGmres (theNewState);
01317             } // end of restarting
01318             //
01319             // We returned from iterate(), but none of our status
01320             // tests Passed.  Something is wrong, and it is probably
01321             // our fault.
01322             //
01323             else {
01324               TEUCHOS_TEST_FOR_EXCEPTION(
01325                 true, std::logic_error,
01326                 "Belos::GmresPolySolMgr::solve: Invalid return from "
01327                 "BlockGmresIter::iterate().  Please report this bug "
01328                 "to the Belos developers.");
01329             }
01330           }
01331           catch (const GmresIterationOrthoFailure& e) {
01332             // If the block size is not one, it's not considered a lucky breakdown.
01333             if (blockSize_ != 1) {
01334               printer_->stream(Errors)
01335                 << "Error! Caught std::exception in BlockGmresIter::iterate() "
01336                 << "at iteration " << block_gmres_iter->getNumIters()
01337                 << std::endl << e.what() << std::endl;
01338               if (convTest_->getStatus() != Passed) {
01339                 isConverged = false;
01340               }
01341               break;
01342             }
01343             else {
01344               // If the block size is one, try to recover the most
01345               // recent least-squares solution
01346               block_gmres_iter->updateLSQR (block_gmres_iter->getCurSubspaceDim ());
01347 
01348               // Check to see if the most recent least-squares
01349               // solution yielded convergence.
01350               sTest_->checkStatus (&*block_gmres_iter);
01351               if (convTest_->getStatus() != Passed) {
01352                 isConverged = false;
01353               }
01354               break;
01355             }
01356           }
01357           catch (const std::exception &e) {
01358             printer_->stream(Errors)
01359               << "Error! Caught std::exception in BlockGmresIter::iterate() "
01360               << "at iteration " << block_gmres_iter->getNumIters() << std::endl
01361               << e.what() << std::endl;
01362             throw;
01363           }
01364         }
01365 
01366         // Compute the current solution.  Update the linear problem.
01367         // Attempt to get the current solution from the residual
01368         // status test, if it has one.
01369         if (! Teuchos::is_null (expConvTest_->getSolution ()) ) {
01370           RCP<MV> newX = expConvTest_->getSolution ();
01371           RCP<MV> curX = problem_->getCurrLHSVec ();
01372           MVT::MvAddMv (STS::zero (), *newX, STS::one (), *newX, *curX);
01373         }
01374         else {
01375           RCP<MV> update = block_gmres_iter->getCurrentUpdate ();
01376           problem_->updateSolution (update, true);
01377         }
01378 
01379         // Inform the linear problem that we are finished with this block linear system.
01380         problem_->setCurrLS ();
01381 
01382         // Update indices for the linear systems to be solved.
01383         startPtr += numCurrRHS;
01384         numRHS2Solve -= numCurrRHS;
01385         if (numRHS2Solve > 0) {
01386           numCurrRHS = (numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01387 
01388           currIdx.resize (blockSize_);
01389           for (int i=0; i<numCurrRHS; ++i) {
01390             currIdx[i] = startPtr+i;
01391           }
01392           for (int i=numCurrRHS; i<blockSize_; ++i) {
01393             currIdx[i] = -1;
01394           }
01395 
01396           // Set the next indices.
01397           problem_->setLSIndex( currIdx );
01398         }
01399         else {
01400           currIdx.resize( numRHS2Solve );
01401         }
01402 
01403       }// while ( numRHS2Solve > 0 )
01404 
01405       // print final summary
01406       sTest_->print( printer_->stream(FinalSummary) );
01407 
01408     } // if (strictConvTol_)
01409   } // timing block
01410 
01411   // print timing information
01412 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01413   // Calling summarize() can be expensive, so don't call unless the
01414   // user wants to print out timing details.  summarize() will do all
01415   // the work even if it's passed a "black hole" output stream.
01416   if (verbosity_ & TimingDetails)
01417     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01418 #endif
01419 
01420   if (!isConverged || loaDetected_) {
01421     return Unconverged; // return from GmresPolySolMgr::solve()
01422   }
01423   return Converged; // return from GmresPolySolMgr::solve()
01424 }
01425 
01426 
01427 template<class ScalarType, class MV, class OP>
01428 std::string GmresPolySolMgr<ScalarType,MV,OP>::description () const
01429 {
01430   std::ostringstream out;
01431 
01432   out << "\"Belos::GmresPolySolMgr\": {"
01433       << "ScalarType: " << Teuchos::TypeNameTraits<ScalarType>::name ()
01434       << ", Ortho Type: " << orthoType_
01435       << ", Block Size: " << blockSize_
01436       << ", Num Blocks: " << numBlocks_
01437       << ", Max Restarts: " << maxRestarts_;
01438   out << "}";
01439   return out.str ();
01440 }
01441 
01442 } // namespace Belos
01443 
01444 #endif // BELOS_GMRES_POLY_SOLMGR_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines