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