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