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