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