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