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