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