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