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