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