Belos Version of the Day
BelosPseudoBlockGmresSolMgr.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_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
00043 #define BELOS_PSEUDO_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 "BelosPseudoBlockGmresIter.hpp"
00056 #include "BelosDGKSOrthoManager.hpp"
00057 #include "BelosICGSOrthoManager.hpp"
00058 #include "BelosIMGSOrthoManager.hpp"
00059 #ifdef HAVE_BELOS_TSQR
00060 #  include "BelosTsqrOrthoManager.hpp"
00061 #endif // HAVE_BELOS_TSQR
00062 #include "BelosStatusTestMaxIters.hpp"
00063 #include "BelosStatusTestGenResNorm.hpp"
00064 #include "BelosStatusTestImpResNorm.hpp"
00065 #include "BelosStatusTestCombo.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 
00080 namespace Belos {
00081 
00083 
00084 
00091   class PseudoBlockGmresSolMgrLinearProblemFailure : public BelosError {public:
00092     PseudoBlockGmresSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00093     {}};
00094 
00101   class PseudoBlockGmresSolMgrOrthoFailure : public BelosError {public:
00102     PseudoBlockGmresSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00103     {}};
00104 
00128   template<class ScalarType, class MV, class OP>
00129   class PseudoBlockGmresSolMgr : public SolverManager<ScalarType,MV,OP> {
00130 
00131   private:
00132     typedef MultiVecTraits<ScalarType,MV> MVT;
00133     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00134     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00135     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00136     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00137 
00138   public:
00139 
00141 
00142 
00150     PseudoBlockGmresSolMgr();
00151 
00255     PseudoBlockGmresSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00256                             const Teuchos::RCP<Teuchos::ParameterList> &pl );
00257 
00259     virtual ~PseudoBlockGmresSolMgr() {};
00261 
00263 
00264 
00265     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00266       return *problem_;
00267     }
00268 
00270     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00271 
00273     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00274 
00280     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00281       return Teuchos::tuple(timerSolve_);
00282     }
00283 
00294     MagnitudeType achievedTol() const {
00295       return achievedTol_;
00296     }
00297 
00299     int getNumIters() const {
00300       return numIters_;
00301     }
00302 
00358     bool isLOADetected() const { return loaDetected_; }
00359 
00361 
00363 
00364 
00366     void setProblem (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem) {
00367       problem_ = problem;
00368     }
00369 
00371     void setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params);
00372 
00379     virtual void setUserConvStatusTest(
00380       const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &userConvStatusTest
00381       );
00382 
00384 
00386 
00387 
00391     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00393 
00395 
00396 
00414     ReturnType solve();
00415 
00417 
00420 
00427     void
00428     describe (Teuchos::FancyOStream& out,
00429               const Teuchos::EVerbosityLevel verbLevel =
00430               Teuchos::Describable::verbLevel_default) const;
00431 
00433     std::string description () const;
00434 
00436 
00437   private:
00438 
00453     bool checkStatusTest();
00454 
00456     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00457 
00458     // Output manager.
00459     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00460     Teuchos::RCP<std::ostream> outputStream_;
00461 
00462     // Status tests.
00463     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > userConvStatusTest_;
00464     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00465     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00466     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00467     Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > impConvTest_, expConvTest_;
00468     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00469 
00470     // Orthogonalization manager.
00471     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
00472 
00473      // Current parameter list.
00474     Teuchos::RCP<Teuchos::ParameterList> params_;
00475 
00476     // Default solver values.
00477     static const MagnitudeType convtol_default_;
00478     static const MagnitudeType orthoKappa_default_;
00479     static const int maxRestarts_default_;
00480     static const int maxIters_default_;
00481     static const bool showMaxResNormOnly_default_;
00482     static const int blockSize_default_;
00483     static const int numBlocks_default_;
00484     static const int verbosity_default_;
00485     static const int outputStyle_default_;
00486     static const int outputFreq_default_;
00487     static const int defQuorum_default_;
00488     static const std::string impResScale_default_;
00489     static const std::string expResScale_default_;
00490     static const std::string label_default_;
00491     static const std::string orthoType_default_;
00492     static const Teuchos::RCP<std::ostream> outputStream_default_;
00493 
00494     // Current solver values.
00495     MagnitudeType convtol_, orthoKappa_, achievedTol_;
00496     int maxRestarts_, maxIters_, numIters_;
00497     int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_, defQuorum_;
00498     bool showMaxResNormOnly_;
00499     std::string orthoType_;
00500     std::string impResScale_, expResScale_;
00501 
00502     // Timers.
00503     std::string label_;
00504     Teuchos::RCP<Teuchos::Time> timerSolve_;
00505 
00506     // Internal state variables.
00507     bool isSet_, isSTSet_, expResTest_;
00508     bool loaDetected_;
00509   };
00510 
00511 
00512 // Default solver values.
00513 template<class ScalarType, class MV, class OP>
00514 const typename PseudoBlockGmresSolMgr<ScalarType,MV,OP>::MagnitudeType PseudoBlockGmresSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00515 
00516 template<class ScalarType, class MV, class OP>
00517 const typename PseudoBlockGmresSolMgr<ScalarType,MV,OP>::MagnitudeType PseudoBlockGmresSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00518 
00519 template<class ScalarType, class MV, class OP>
00520 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::maxRestarts_default_ = 20;
00521 
00522 template<class ScalarType, class MV, class OP>
00523 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00524 
00525 template<class ScalarType, class MV, class OP>
00526 const bool PseudoBlockGmresSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00527 
00528 template<class ScalarType, class MV, class OP>
00529 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00530 
00531 template<class ScalarType, class MV, class OP>
00532 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::numBlocks_default_ = 300;
00533 
00534 template<class ScalarType, class MV, class OP>
00535 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00536 
00537 template<class ScalarType, class MV, class OP>
00538 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00539 
00540 template<class ScalarType, class MV, class OP>
00541 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00542 
00543 template<class ScalarType, class MV, class OP>
00544 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::defQuorum_default_ = 1;
00545 
00546 template<class ScalarType, class MV, class OP>
00547 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of Preconditioned Initial Residual";
00548 
00549 template<class ScalarType, class MV, class OP>
00550 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of Initial Residual";
00551 
00552 template<class ScalarType, class MV, class OP>
00553 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00554 
00555 template<class ScalarType, class MV, class OP>
00556 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00557 
00558 template<class ScalarType, class MV, class OP>
00559 const Teuchos::RCP<std::ostream> PseudoBlockGmresSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00560 
00561 
00562 // Empty Constructor
00563 template<class ScalarType, class MV, class OP>
00564 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::PseudoBlockGmresSolMgr() :
00565   outputStream_(outputStream_default_),
00566   convtol_(convtol_default_),
00567   orthoKappa_(orthoKappa_default_),
00568   achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
00569   maxRestarts_(maxRestarts_default_),
00570   maxIters_(maxIters_default_),
00571   numIters_(0),
00572   blockSize_(blockSize_default_),
00573   numBlocks_(numBlocks_default_),
00574   verbosity_(verbosity_default_),
00575   outputStyle_(outputStyle_default_),
00576   outputFreq_(outputFreq_default_),
00577   defQuorum_(defQuorum_default_),
00578   showMaxResNormOnly_(showMaxResNormOnly_default_),
00579   orthoType_(orthoType_default_),
00580   impResScale_(impResScale_default_),
00581   expResScale_(expResScale_default_),
00582   label_(label_default_),
00583   isSet_(false),
00584   isSTSet_(false),
00585   expResTest_(false),
00586   loaDetected_(false)
00587 {}
00588 
00589 // Basic Constructor
00590 template<class ScalarType, class MV, class OP>
00591 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::
00592 PseudoBlockGmresSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00593                         const Teuchos::RCP<Teuchos::ParameterList> &pl) :
00594   problem_(problem),
00595   outputStream_(outputStream_default_),
00596   convtol_(convtol_default_),
00597   orthoKappa_(orthoKappa_default_),
00598   achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
00599   maxRestarts_(maxRestarts_default_),
00600   maxIters_(maxIters_default_),
00601   numIters_(0),
00602   blockSize_(blockSize_default_),
00603   numBlocks_(numBlocks_default_),
00604   verbosity_(verbosity_default_),
00605   outputStyle_(outputStyle_default_),
00606   outputFreq_(outputFreq_default_),
00607   defQuorum_(defQuorum_default_),
00608   showMaxResNormOnly_(showMaxResNormOnly_default_),
00609   orthoType_(orthoType_default_),
00610   impResScale_(impResScale_default_),
00611   expResScale_(expResScale_default_),
00612   label_(label_default_),
00613   isSet_(false),
00614   isSTSet_(false),
00615   expResTest_(false),
00616   loaDetected_(false)
00617 {
00618   TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00619 
00620   // If the parameter list pointer is null, then set the current parameters to the default parameter list.
00621   if (!is_null(pl)) {
00622     // Set the parameters using the list that was passed in.
00623     setParameters( pl );
00624   }
00625 }
00626 
00627 template<class ScalarType, class MV, class OP>
00628 void
00629 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::
00630 setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
00631 {
00632   using Teuchos::ParameterList;
00633   using Teuchos::parameterList;
00634   using Teuchos::rcp;
00635   using Teuchos::rcp_dynamic_cast;
00636 
00637   // Create the internal parameter list if one doesn't already exist.
00638   if (params_ == Teuchos::null) {
00639     params_ = parameterList (*getValidParameters ());
00640   } else {
00641     params->validateParameters (*getValidParameters ());
00642   }
00643 
00644   // Check for maximum number of restarts
00645   if (params->isParameter ("Maximum Restarts")) {
00646     maxRestarts_ = params->get ("Maximum Restarts", maxRestarts_default_);
00647 
00648     // Update parameter in our list.
00649     params_->set ("Maximum Restarts", maxRestarts_);
00650   }
00651 
00652   // Check for maximum number of iterations
00653   if (params->isParameter ("Maximum Iterations")) {
00654     maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
00655 
00656     // Update parameter in our list and in status test.
00657     params_->set ("Maximum Iterations", maxIters_);
00658     if (! maxIterTest_.is_null ()) {
00659       maxIterTest_->setMaxIters (maxIters_);
00660     }
00661   }
00662 
00663   // Check for blocksize
00664   if (params->isParameter ("Block Size")) {
00665     blockSize_ = params->get ("Block Size", blockSize_default_);
00666     TEUCHOS_TEST_FOR_EXCEPTION(
00667       blockSize_ <= 0, std::invalid_argument,
00668       "Belos::PseudoBlockGmresSolMgr::setParameters: "
00669       "The \"Block Size\" parameter must be strictly positive, "
00670       "but you specified a value of " << blockSize_ << ".");
00671 
00672     // Update parameter in our list.
00673     params_->set ("Block Size", blockSize_);
00674   }
00675 
00676   // Check for the maximum number of blocks.
00677   if (params->isParameter ("Num Blocks")) {
00678     numBlocks_ = params->get ("Num Blocks", numBlocks_default_);
00679     TEUCHOS_TEST_FOR_EXCEPTION(
00680       numBlocks_ <= 0, std::invalid_argument,
00681       "Belos::PseudoBlockGmresSolMgr::setParameters: "
00682       "The \"Num Blocks\" parameter must be strictly positive, "
00683       "but you specified a value of " << numBlocks_ << ".");
00684 
00685     // Update parameter in our list.
00686     params_->set ("Num Blocks", numBlocks_);
00687   }
00688 
00689   // Check to see if the timer label changed.
00690   if (params->isParameter ("Timer Label")) {
00691     const std::string tempLabel = params->get ("Timer Label", label_default_);
00692 
00693     // Update parameter in our list and solver timer
00694     if (tempLabel != label_) {
00695       label_ = tempLabel;
00696       params_->set ("Timer Label", label_);
00697       const std::string solveLabel =
00698         label_ + ": PseudoBlockGmresSolMgr total solve time";
00699 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00700       timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
00701 #endif // BELOS_TEUCHOS_TIME_MONITOR
00702       if (ortho_ != Teuchos::null) {
00703         ortho_->setLabel( label_ );
00704       }
00705     }
00706   }
00707 
00708   // Check if the orthogonalization changed.
00709   if (params->isParameter ("Orthogonalization")) {
00710     std::string tempOrthoType = params->get ("Orthogonalization", orthoType_default_);
00711 #ifdef HAVE_BELOS_TSQR
00712     TEUCHOS_TEST_FOR_EXCEPTION(
00713       tempOrthoType != "DGKS" && tempOrthoType != "ICGS" &&
00714       tempOrthoType != "IMGS" && tempOrthoType != "TSQR",
00715       std::invalid_argument,
00716       "Belos::PseudoBlockGmresSolMgr::setParameters: "
00717       "The \"Orthogonalization\" parameter must be one of \"DGKS\", \"ICGS\", "
00718       "\"IMGS\", or \"TSQR\".");
00719 #else
00720     TEUCHOS_TEST_FOR_EXCEPTION(
00721       tempOrthoType != "DGKS" && tempOrthoType != "ICGS" &&
00722       tempOrthoType != "IMGS",
00723       std::invalid_argument,
00724       "Belos::PseudoBlockGmresSolMgr::setParameters: "
00725       "The \"Orthogonalization\" parameter must be one of \"DGKS\", \"ICGS\", "
00726       "or \"IMGS\".");
00727 #endif // HAVE_BELOS_TSQR
00728 
00729     if (tempOrthoType != orthoType_) {
00730       orthoType_ = tempOrthoType;
00731       // Create orthogonalization manager
00732       if (orthoType_ == "DGKS") {
00733         typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_type;
00734         if (orthoKappa_ <= 0) {
00735           ortho_ = rcp (new ortho_type (label_));
00736         }
00737         else {
00738           ortho_ = rcp (new ortho_type (label_));
00739           rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
00740         }
00741       }
00742       else if (orthoType_ == "ICGS") {
00743         typedef ICGSOrthoManager<ScalarType, MV, OP> ortho_type;
00744         ortho_ = rcp (new ortho_type (label_));
00745       }
00746       else if (orthoType_ == "IMGS") {
00747         typedef IMGSOrthoManager<ScalarType, MV, OP> ortho_type;
00748         ortho_ = rcp (new ortho_type (label_));
00749       }
00750 #ifdef HAVE_BELOS_TSQR
00751       else if (orthoType_ == "TSQR") {
00752         typedef TsqrMatOrthoManager<ScalarType, MV, OP> ortho_type;
00753         ortho_ = rcp (new ortho_type (label_));
00754       }
00755 #endif // HAVE_BELOS_TSQR
00756     }
00757   }
00758 
00759   // Check which orthogonalization constant to use.
00760   if (params->isParameter ("Orthogonalization Constant")) {
00761     orthoKappa_ = params->get ("Orthogonalization Constant", orthoKappa_default_);
00762 
00763     // Update parameter in our list.
00764     params_->set ("Orthogonalization Constant", orthoKappa_);
00765     if (orthoType_ == "DGKS") {
00766       if (orthoKappa_ > 0 && ! ortho_.is_null ()) {
00767         typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_type;
00768         rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
00769       }
00770     }
00771   }
00772 
00773   // Check for a change in verbosity level
00774   if (params->isParameter ("Verbosity")) {
00775     if (Teuchos::isParameterType<int> (*params, "Verbosity")) {
00776       verbosity_ = params->get ("Verbosity", verbosity_default_);
00777     } else {
00778       verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params, "Verbosity");
00779     }
00780 
00781     // Update parameter in our list.
00782     params_->set ("Verbosity", verbosity_);
00783     if (! printer_.is_null ()) {
00784       printer_->setVerbosity (verbosity_);
00785     }
00786   }
00787 
00788   // Check for a change in output style.
00789   if (params->isParameter ("Output Style")) {
00790     if (Teuchos::isParameterType<int> (*params, "Output Style")) {
00791       outputStyle_ = params->get ("Output Style", outputStyle_default_);
00792     } else {
00793       outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params, "Output Style");
00794     }
00795 
00796     // Update parameter in our list.
00797     params_->set ("Output Style", verbosity_);
00798     if (! outputTest_.is_null ()) {
00799       isSTSet_ = false;
00800     }
00801   }
00802 
00803   // output stream
00804   if (params->isParameter ("Output Stream")) {
00805     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> > (*params, "Output Stream");
00806 
00807     // Update parameter in our list.
00808     params_->set("Output Stream", outputStream_);
00809     if (! printer_.is_null ()) {
00810       printer_->setOStream (outputStream_);
00811     }
00812   }
00813 
00814   // frequency level
00815   if (verbosity_ & Belos::StatusTestDetails) {
00816     if (params->isParameter ("Output Frequency")) {
00817       outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
00818     }
00819 
00820     // Update parameter in out list and output status test.
00821     params_->set ("Output Frequency", outputFreq_);
00822     if (! outputTest_.is_null ()) {
00823       outputTest_->setOutputFrequency (outputFreq_);
00824     }
00825   }
00826 
00827   // Create output manager if we need to.
00828   if (printer_.is_null ()) {
00829     printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00830   }
00831 
00832   // Convergence
00833 
00834   // Check for convergence tolerance
00835   if (params->isParameter ("Convergence Tolerance")) {
00836     convtol_ = params->get ("Convergence Tolerance", convtol_default_);
00837 
00838     // Update parameter in our list and residual tests.
00839     params_->set ("Convergence Tolerance", convtol_);
00840     if (! impConvTest_.is_null ()) {
00841       impConvTest_->setTolerance (convtol_);
00842     }
00843     if (! expConvTest_.is_null ()) {
00844       expConvTest_->setTolerance (convtol_);
00845     }
00846   }
00847 
00848   // Check for a change in scaling, if so we need to build new residual tests.
00849   if (params->isParameter ("Implicit Residual Scaling")) {
00850     const std::string tempImpResScale =
00851       Teuchos::getParameter<std::string> (*params, "Implicit Residual Scaling");
00852 
00853     // Only update the scaling if it's different.
00854     if (impResScale_ != tempImpResScale) {
00855       Belos::ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
00856       impResScale_ = tempImpResScale;
00857 
00858       // Update parameter in our list and residual tests
00859       params_->set ("Implicit Residual Scaling", impResScale_);
00860       if (! impConvTest_.is_null ()) {
00861         try {
00862           impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
00863         }
00864         catch (std::exception& e) {
00865           // Make sure the convergence test gets constructed again.
00866           isSTSet_ = false;
00867         }
00868       }
00869     }
00870   }
00871 
00872   if (params->isParameter ("Explicit Residual Scaling")) {
00873     const std::string tempExpResScale =
00874       Teuchos::getParameter<std::string> (*params, "Explicit Residual Scaling");
00875 
00876     // Only update the scaling if it's different.
00877     if (expResScale_ != tempExpResScale) {
00878       Belos::ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
00879       expResScale_ = tempExpResScale;
00880 
00881       // Update parameter in our list and residual tests
00882       params_->set ("Explicit Residual Scaling", expResScale_);
00883       if (! expConvTest_.is_null ()) {
00884         try {
00885           expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
00886         }
00887         catch (std::exception& e) {
00888           // Make sure the convergence test gets constructed again.
00889           isSTSet_ = false;
00890         }
00891       }
00892     }
00893   }
00894 
00895 
00896   if (params->isParameter ("Show Maximum Residual Norm Only")) {
00897     showMaxResNormOnly_ =
00898       Teuchos::getParameter<bool> (*params, "Show Maximum Residual Norm Only");
00899 
00900     // Update parameter in our list and residual tests
00901     params_->set ("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00902     if (! impConvTest_.is_null ()) {
00903       impConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
00904     }
00905     if (! expConvTest_.is_null ()) {
00906       expConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
00907     }
00908   }
00909 
00910   // Create status tests if we need to.
00911 
00912   // Get the deflation quorum, or number of converged systems before deflation is allowed
00913   if (params->isParameter("Deflation Quorum")) {
00914     defQuorum_ = params->get("Deflation Quorum", defQuorum_);
00915     TEUCHOS_TEST_FOR_EXCEPTION(
00916       defQuorum_ > blockSize_, std::invalid_argument,
00917       "Belos::PseudoBlockGmresSolMgr::setParameters: "
00918       "The \"Deflation Quorum\" parameter (= " << defQuorum_ << ") must not be "
00919       "larger than \"Block Size\" (= " << blockSize_ << ").");
00920     params_->set ("Deflation Quorum", defQuorum_);
00921     if (! impConvTest_.is_null ()) {
00922       impConvTest_->setQuorum (defQuorum_);
00923     }
00924     if (! expConvTest_.is_null ()) {
00925       expConvTest_->setQuorum (defQuorum_);
00926     }
00927   }
00928 
00929   // Create orthogonalization manager if we need to.
00930   if (ortho_.is_null ()) {
00931     if (orthoType_ == "DGKS") {
00932       typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_type;
00933       if (orthoKappa_ <= 0) {
00934         ortho_ = rcp (new ortho_type (label_));
00935       }
00936       else {
00937         ortho_ = rcp (new ortho_type (label_));
00938         rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
00939       }
00940     }
00941     else if (orthoType_ == "ICGS") {
00942       typedef ICGSOrthoManager<ScalarType, MV, OP> ortho_type;
00943       ortho_ = rcp (new ortho_type (label_));
00944     }
00945     else if (orthoType_ == "IMGS") {
00946       typedef IMGSOrthoManager<ScalarType, MV, OP> ortho_type;
00947       ortho_ = rcp (new ortho_type (label_));
00948     }
00949 #ifdef HAVE_BELOS_TSQR
00950     else if (orthoType_ == "TSQR") {
00951       typedef TsqrMatOrthoManager<ScalarType, MV, OP> ortho_type;
00952       ortho_ = rcp (new ortho_type (label_));
00953     }
00954 #endif // HAVE_BELOS_TSQR
00955     else {
00956 #ifdef HAVE_BELOS_TSQR
00957       TEUCHOS_TEST_FOR_EXCEPTION(
00958         orthoType_ != "ICGS" && orthoType_ != "DGKS" &&
00959         orthoType_ != "IMGS" && orthoType_ != "TSQR",
00960         std::logic_error,
00961         "Belos::PseudoBlockGmresSolMgr::setParameters(): "
00962         "Invalid orthogonalization type \"" << orthoType_ << "\".");
00963 #else
00964       TEUCHOS_TEST_FOR_EXCEPTION(
00965         orthoType_ != "ICGS" && orthoType_ != "DGKS" &&
00966         orthoType_ != "IMGS",
00967         std::logic_error,
00968         "Belos::PseudoBlockGmresSolMgr::setParameters(): "
00969         "Invalid orthogonalization type \"" << orthoType_ << "\".");
00970 #endif // HAVE_BELOS_TSQR
00971     }
00972   }
00973 
00974   // Create the timer if we need to.
00975   if (timerSolve_ == Teuchos::null) {
00976     std::string solveLabel = label_ + ": PseudoBlockGmresSolMgr total solve time";
00977 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00978     timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
00979 #endif
00980   }
00981 
00982   // Inform the solver manager that the current parameters were set.
00983   isSet_ = true;
00984 }
00985 
00986 
00987 template<class ScalarType, class MV, class OP>
00988 void
00989 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::setUserConvStatusTest(
00990   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &userConvStatusTest
00991   )
00992 {
00993   userConvStatusTest_ = userConvStatusTest;
00994 }
00995 
00996 
00997 template<class ScalarType, class MV, class OP>
00998 Teuchos::RCP<const Teuchos::ParameterList>
00999 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::getValidParameters() const
01000 {
01001   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
01002   if (is_null(validPL)) {
01003     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
01004   // Set all the valid parameters and their default values.
01005     pl= Teuchos::rcp( new Teuchos::ParameterList() );
01006     pl->set("Convergence Tolerance", convtol_default_,
01007       "The relative residual tolerance that needs to be achieved by the\n"
01008       "iterative solver in order for the linear system to be declared converged.");
01009     pl->set("Maximum Restarts", maxRestarts_default_,
01010       "The maximum number of restarts allowed for each\n"
01011       "set of RHS solved.");
01012     pl->set("Maximum Iterations", maxIters_default_,
01013       "The maximum number of block iterations allowed for each\n"
01014       "set of RHS solved.");
01015     pl->set("Num Blocks", numBlocks_default_,
01016       "The maximum number of vectors allowed in the Krylov subspace\n"
01017       "for each set of RHS solved.");
01018     pl->set("Block Size", blockSize_default_,
01019       "The number of RHS solved simultaneously.");
01020     pl->set("Verbosity", verbosity_default_,
01021       "What type(s) of solver information should be outputted\n"
01022       "to the output stream.");
01023     pl->set("Output Style", outputStyle_default_,
01024       "What style is used for the solver information outputted\n"
01025       "to the output stream.");
01026     pl->set("Output Frequency", outputFreq_default_,
01027       "How often convergence information should be outputted\n"
01028       "to the output stream.");
01029     pl->set("Deflation Quorum", defQuorum_default_,
01030       "The number of linear systems that need to converge before\n"
01031       "they are deflated.  This number should be <= block size.");
01032     pl->set("Output Stream", outputStream_default_,
01033       "A reference-counted pointer to the output stream where all\n"
01034       "solver output is sent.");
01035     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
01036       "When convergence information is printed, only show the maximum\n"
01037       "relative residual norm when the block size is greater than one.");
01038     pl->set("Implicit Residual Scaling", impResScale_default_,
01039       "The type of scaling used in the implicit residual convergence test.");
01040     pl->set("Explicit Residual Scaling", expResScale_default_,
01041       "The type of scaling used in the explicit residual convergence test.");
01042     pl->set("Timer Label", label_default_,
01043       "The string to use as a prefix for the timer labels.");
01044     //  defaultParams_->set("Restart Timers", restartTimers_);
01045     pl->set("Orthogonalization", orthoType_default_,
01046       "The type of orthogonalization to use: DGKS, ICGS, IMGS.");
01047     pl->set("Orthogonalization Constant",orthoKappa_default_,
01048       "The constant used by DGKS orthogonalization to determine\n"
01049       "whether another step of classical Gram-Schmidt is necessary.");
01050     validPL = pl;
01051   }
01052   return validPL;
01053 }
01054 
01055 // Check the status test versus the defined linear problem
01056 template<class ScalarType, class MV, class OP>
01057 bool PseudoBlockGmresSolMgr<ScalarType,MV,OP>::checkStatusTest() {
01058 
01059   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
01060   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestGenResNorm_t;
01061   typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP>  StatusTestImpResNorm_t;
01062 
01063   // Basic test checks maximum iterations and native residual.
01064   maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
01065 
01066   // If there is a left preconditioner, we create a combined status test that checks the implicit
01067   // and then explicit residual norm to see if we have convergence.
01068   if ( !Teuchos::is_null(problem_->getLeftPrec()) ) {
01069     expResTest_ = true;
01070   }
01071 
01072   if (expResTest_) {
01073 
01074     // Implicit residual test, using the native residual to determine if convergence was achieved.
01075     Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
01076       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
01077     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
01078     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
01079     impConvTest_ = tmpImpConvTest;
01080 
01081     // Explicit residual test once the native residual is below the tolerance
01082     Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
01083       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
01084     tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
01085     tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
01086     tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
01087     expConvTest_ = tmpExpConvTest;
01088 
01089     // The convergence test is a combination of the "cheap" implicit test and explicit test.
01090     convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
01091   }
01092   else {
01093 
01094     // Implicit residual test, using the native residual to determine if convergence was achieved.
01095     // Use test that checks for loss of accuracy.
01096     Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
01097       Teuchos::rcp( new StatusTestImpResNorm_t( convtol_, defQuorum_ ) );
01098     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
01099     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
01100     impConvTest_ = tmpImpConvTest;
01101 
01102     // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
01103     expConvTest_ = impConvTest_;
01104     convTest_ = impConvTest_;
01105   }
01106 
01107   if (nonnull(userConvStatusTest_) ) {
01108     // Override the overall convergence test with the users convergence test
01109     convTest_ = Teuchos::rcp(
01110       new StatusTestCombo_t( StatusTestCombo_t::SEQ, convTest_, userConvStatusTest_ ) );
01111     // NOTE: Above, you have to run the other convergence tests also because
01112     // the logic in this class depends on it.  This is very unfortunate.
01113   }
01114 
01115   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
01116 
01117   // Create the status test output class.
01118   // This class manages and formats the output from the status test.
01119   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
01120   outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
01121 
01122   // Set the solver string for the output test
01123   std::string solverDesc = " Pseudo Block Gmres ";
01124   outputTest_->setSolverDesc( solverDesc );
01125 
01126 
01127   // The status test is now set.
01128   isSTSet_ = true;
01129 
01130   return false;
01131 }
01132 
01133 
01134 // solve()
01135 template<class ScalarType, class MV, class OP>
01136 ReturnType PseudoBlockGmresSolMgr<ScalarType,MV,OP>::solve() {
01137 
01138   // Set the current parameters if they were not set before.
01139   // NOTE:  This may occur if the user generated the solver manager with the default constructor and
01140   // then didn't set any parameters using setParameters().
01141   if (!isSet_) { setParameters( params_ ); }
01142 
01143   Teuchos::BLAS<int,ScalarType> blas;
01144 
01145   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockGmresSolMgrLinearProblemFailure,
01146                      "Belos::PseudoBlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
01147 
01148   // Check if we have to create the status tests.
01149   if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
01150     TEUCHOS_TEST_FOR_EXCEPTION( checkStatusTest(),PseudoBlockGmresSolMgrLinearProblemFailure,
01151       "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
01152   }
01153 
01154   // Create indices for the linear systems to be solved.
01155   int startPtr = 0;
01156   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
01157   int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01158 
01159   std::vector<int> currIdx( numCurrRHS );
01160   blockSize_ = numCurrRHS;
01161   for (int i=0; i<numCurrRHS; ++i)
01162     { currIdx[i] = startPtr+i; }
01163 
01164   // Inform the linear problem of the current linear system to solve.
01165   problem_->setLSIndex( currIdx );
01166 
01168   // Parameter list
01169   Teuchos::ParameterList plist;
01170   plist.set("Num Blocks",numBlocks_);
01171 
01172   // Reset the status test.
01173   outputTest_->reset();
01174   loaDetected_ = false;
01175 
01176   // Assume convergence is achieved, then let any failed convergence set this to false.
01177   bool isConverged = true;
01178 
01180   // BlockGmres solver
01181 
01182   Teuchos::RCP<PseudoBlockGmresIter<ScalarType,MV,OP> > block_gmres_iter
01183     = Teuchos::rcp( new PseudoBlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
01184 
01185   // Enter solve() iterations
01186   {
01187 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01188     Teuchos::TimeMonitor slvtimer(*timerSolve_);
01189 #endif
01190 
01191     while ( numRHS2Solve > 0 ) {
01192 
01193       // Reset the active / converged vectors from this block
01194       std::vector<int> convRHSIdx;
01195       std::vector<int> currRHSIdx( currIdx );
01196       currRHSIdx.resize(numCurrRHS);
01197 
01198       // Set the current number of blocks with the pseudo Gmres iteration.
01199       block_gmres_iter->setNumBlocks( numBlocks_ );
01200 
01201       // Reset the number of iterations.
01202       block_gmres_iter->resetNumIters();
01203 
01204       // Reset the number of calls that the status test output knows about.
01205       outputTest_->resetNumCalls();
01206 
01207       // Get a new state struct and initialize the solver.
01208       PseudoBlockGmresIterState<ScalarType,MV> newState;
01209 
01210       // Create the first block in the current Krylov basis for each right-hand side.
01211       std::vector<int> index(1);
01212       Teuchos::RCP<MV> tmpV, R_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
01213       newState.V.resize( blockSize_ );
01214       newState.Z.resize( blockSize_ );
01215       for (int i=0; i<blockSize_; ++i) {
01216         index[0]=i;
01217         tmpV = MVT::CloneViewNonConst( *R_0, index );
01218 
01219         // Get a matrix to hold the orthonormalization coefficients.
01220         Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
01221           = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
01222 
01223         // Orthonormalize the new V_0
01224         int rank = ortho_->normalize( *tmpV, tmpZ );
01225         TEUCHOS_TEST_FOR_EXCEPTION(rank != 1, PseudoBlockGmresSolMgrOrthoFailure,
01226             "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
01227 
01228         newState.V[i] = tmpV;
01229         newState.Z[i] = tmpZ;
01230       }
01231 
01232       newState.curDim = 0;
01233       block_gmres_iter->initialize(newState);
01234       int numRestarts = 0;
01235 
01236       while(1) {
01237 
01238         // tell block_gmres_iter to iterate
01239         try {
01240           block_gmres_iter->iterate();
01241 
01243           //
01244           // check convergence first
01245           //
01247           if ( convTest_->getStatus() == Passed ) {
01248 
01249             if ( expConvTest_->getLOADetected() ) {
01250               // Oops!  There was a loss of accuracy (LOA) for one or
01251               // more right-hand sides.  That means the implicit
01252               // (a.k.a. "native") residuals claim convergence,
01253               // whereas the explicit (hence expConvTest_, i.e.,
01254               // "explicit convergence test") residuals have not
01255               // converged.
01256               //
01257               // We choose to deal with this situation by deflating
01258               // out the affected right-hand sides and moving on.
01259               loaDetected_ = true;
01260               printer_->stream(Warnings) <<
01261                 "Belos::PseudoBlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
01262               isConverged = false;
01263             }
01264 
01265             // Figure out which linear systems converged.
01266             std::vector<int> convIdx = expConvTest_->convIndices();
01267 
01268             // If the number of converged linear systems is equal to the
01269             // number of current linear systems, then we are done with this block.
01270             if (convIdx.size() == currRHSIdx.size())
01271               break;  // break from while(1){block_gmres_iter->iterate()}
01272 
01273             // Get a new state struct and initialize the solver.
01274             PseudoBlockGmresIterState<ScalarType,MV> defState;
01275 
01276             // Inform the linear problem that we are finished with this current linear system.
01277             problem_->setCurrLS();
01278 
01279             // Get the state.
01280             PseudoBlockGmresIterState<ScalarType,MV> oldState = block_gmres_iter->getState();
01281             int curDim = oldState.curDim;
01282 
01283             // Get a new state struct and reset currRHSIdx to have the right-hand sides that
01284             // are left to converge for this block.
01285             int have = 0;
01286             std::vector<int> oldRHSIdx( currRHSIdx );
01287             std::vector<int> defRHSIdx;
01288             for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
01289               bool found = false;
01290               for (unsigned int j=0; j<convIdx.size(); ++j) {
01291                 if (currRHSIdx[i] == convIdx[j]) {
01292                   found = true;
01293                   break;
01294                 }
01295               }
01296               if (found) {
01297                 defRHSIdx.push_back( i );
01298               }
01299               else {
01300                 defState.V.push_back( Teuchos::rcp_const_cast<MV>( oldState.V[i] ) );
01301                 defState.Z.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.Z[i] ) );
01302                 defState.H.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseMatrix<int,ScalarType> >( oldState.H[i] ) );
01303                 defState.sn.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.sn[i] ) );
01304                 defState.cs.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,MagnitudeType> >(oldState.cs[i] ) );
01305                 currRHSIdx[have] = currRHSIdx[i];
01306                 have++;
01307               }
01308             }
01309             defRHSIdx.resize(currRHSIdx.size()-have);
01310             currRHSIdx.resize(have);
01311 
01312             // Compute the current solution that needs to be deflated if this solver has taken any steps.
01313             if (curDim) {
01314               Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01315               Teuchos::RCP<MV> defUpdate = MVT::CloneViewNonConst( *update, defRHSIdx );
01316 
01317               // Set the deflated indices so we can update the solution.
01318               problem_->setLSIndex( convIdx );
01319 
01320               // Update the linear problem.
01321               problem_->updateSolution( defUpdate, true );
01322             }
01323 
01324             // Set the remaining indices after deflation.
01325             problem_->setLSIndex( currRHSIdx );
01326 
01327             // Set the dimension of the subspace, which is the same as the old subspace size.
01328             defState.curDim = curDim;
01329 
01330             // Initialize the solver with the deflated system.
01331             block_gmres_iter->initialize(defState);
01332           }
01334           //
01335           // check for maximum iterations
01336           //
01338           else if ( maxIterTest_->getStatus() == Passed ) {
01339             // we don't have convergence
01340             isConverged = false;
01341             break;  // break from while(1){block_gmres_iter->iterate()}
01342           }
01344           //
01345           // check for restarting, i.e. the subspace is full
01346           //
01348           else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
01349 
01350             if ( numRestarts >= maxRestarts_ ) {
01351               isConverged = false;
01352               break; // break from while(1){block_gmres_iter->iterate()}
01353             }
01354             numRestarts++;
01355 
01356             printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
01357 
01358             // Update the linear problem.
01359             Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01360             problem_->updateSolution( update, true );
01361 
01362             // Get the state.
01363             PseudoBlockGmresIterState<ScalarType,MV> oldState = block_gmres_iter->getState();
01364 
01365             // Set the new state.
01366             PseudoBlockGmresIterState<ScalarType,MV> newstate;
01367             newstate.V.resize(currRHSIdx.size());
01368             newstate.Z.resize(currRHSIdx.size());
01369 
01370             // Compute the restart vectors
01371             // NOTE: Force the linear problem to update the current residual since the solution was updated.
01372             R_0 = MVT::Clone( *(problem_->getInitPrecResVec()), currRHSIdx.size() );
01373             problem_->computeCurrPrecResVec( &*R_0 );
01374             for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
01375               index[0] = i;  // index(1) vector declared on line 891
01376 
01377               tmpV = MVT::CloneViewNonConst( *R_0, index );
01378 
01379               // Get a matrix to hold the orthonormalization coefficients.
01380               Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
01381                 = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
01382 
01383               // Orthonormalize the new V_0
01384               int rank = ortho_->normalize( *tmpV, tmpZ );
01385               TEUCHOS_TEST_FOR_EXCEPTION(rank != 1 ,PseudoBlockGmresSolMgrOrthoFailure,
01386                   "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after the restart.");
01387 
01388               newstate.V[i] = tmpV;
01389               newstate.Z[i] = tmpZ;
01390             }
01391 
01392             // Initialize the solver.
01393             newstate.curDim = 0;
01394             block_gmres_iter->initialize(newstate);
01395 
01396           } // end of restarting
01397 
01399           //
01400           // we returned from iterate(), but none of our status tests Passed.
01401           // something is wrong, and it is probably our fault.
01402           //
01404 
01405           else {
01406             TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
01407                 "Belos::PseudoBlockGmresSolMgr::solve(): Invalid return from PseudoBlockGmresIter::iterate().");
01408           }
01409         }
01410         catch (const PseudoBlockGmresIterOrthoFailure &e) {
01411 
01412           // Try to recover the most recent least-squares solution
01413           block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
01414 
01415           // Check to see if the most recent least-squares solution yielded convergence.
01416           sTest_->checkStatus( &*block_gmres_iter );
01417           if (convTest_->getStatus() != Passed)
01418             isConverged = false;
01419           break;
01420         }
01421         catch (const std::exception &e) {
01422           printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockGmresIter::iterate() at iteration "
01423                                    << block_gmres_iter->getNumIters() << std::endl
01424                                    << e.what() << std::endl;
01425           throw;
01426         }
01427       }
01428 
01429       // Compute the current solution.
01430       // Update the linear problem.
01431       if (nonnull(userConvStatusTest_)) {
01432         //std::cout << "\nnonnull(userConvStatusTest_)\n";
01433         Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01434         problem_->updateSolution( update, true );
01435       }
01436       else if (nonnull(expConvTest_->getSolution())) {
01437         //std::cout << "\nexpConvTest_->getSolution()\n";
01438         Teuchos::RCP<MV> newX = expConvTest_->getSolution();
01439         Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
01440         MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
01441       }
01442       else {
01443         //std::cout << "\nblock_gmres_iter->getCurrentUpdate()\n";
01444         Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01445         problem_->updateSolution( update, true );
01446       }
01447 
01448       // Inform the linear problem that we are finished with this block linear system.
01449       problem_->setCurrLS();
01450 
01451       // Update indices for the linear systems to be solved.
01452       startPtr += numCurrRHS;
01453       numRHS2Solve -= numCurrRHS;
01454       if ( numRHS2Solve > 0 ) {
01455         numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01456 
01457         blockSize_ = numCurrRHS;
01458         currIdx.resize( numCurrRHS  );
01459         for (int i=0; i<numCurrRHS; ++i)
01460         { currIdx[i] = startPtr+i; }
01461 
01462         // Adapt the status test quorum if we need to.
01463         if (defQuorum_ > blockSize_) {
01464           if (impConvTest_ != Teuchos::null)
01465             impConvTest_->setQuorum( blockSize_ );
01466           if (expConvTest_ != Teuchos::null)
01467             expConvTest_->setQuorum( blockSize_ );
01468         }
01469 
01470         // Set the next indices.
01471         problem_->setLSIndex( currIdx );
01472       }
01473       else {
01474         currIdx.resize( numRHS2Solve );
01475       }
01476 
01477     }// while ( numRHS2Solve > 0 )
01478 
01479   }
01480 
01481   // print final summary
01482   sTest_->print( printer_->stream(FinalSummary) );
01483 
01484   // print timing information
01485 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01486   // Calling summarize() can be expensive, so don't call unless the
01487   // user wants to print out timing details.  summarize() will do all
01488   // the work even if it's passed a "black hole" output stream.
01489   if (verbosity_ & TimingDetails)
01490     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01491 #endif
01492 
01493   // get iteration information for this solve
01494   numIters_ = maxIterTest_->getNumIters();
01495 
01496   // Save the convergence test value ("achieved tolerance") for this
01497   // solve.  For this solver, convTest_ may either be a single
01498   // residual norm test, or a combination of two residual norm tests.
01499   // In the latter case, the master convergence test convTest_ is a
01500   // SEQ combo of the implicit resp. explicit tests.  If the implicit
01501   // test never passes, then the explicit test won't ever be executed.
01502   // This manifests as expConvTest_->getTestValue()->size() < 1.  We
01503   // deal with this case by using the values returned by
01504   // impConvTest_->getTestValue().
01505   {
01506     // We'll fetch the vector of residual norms one way or the other.
01507     const std::vector<MagnitudeType>* pTestValues = NULL;
01508     if (expResTest_) {
01509       pTestValues = expConvTest_->getTestValue();
01510       if (pTestValues == NULL || pTestValues->size() < 1) {
01511         pTestValues = impConvTest_->getTestValue();
01512       }
01513     }
01514     else {
01515       // Only the implicit residual norm test is being used.
01516       pTestValues = impConvTest_->getTestValue();
01517     }
01518     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
01519       "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's "
01520       "getTestValue() method returned NULL.  Please report this bug to the "
01521       "Belos developers.");
01522     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
01523       "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's "
01524       "getTestValue() method returned a vector of length zero.  Please report "
01525       "this bug to the Belos developers.");
01526 
01527     // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
01528     // achieved tolerances for all vectors in the current solve(), or
01529     // just for the vectors from the last deflation?
01530     achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
01531   }
01532 
01533   if (!isConverged || loaDetected_) {
01534     return Unconverged; // return from PseudoBlockGmresSolMgr::solve()
01535   }
01536   return Converged; // return from PseudoBlockGmresSolMgr::solve()
01537 }
01538 
01539 
01540 template<class ScalarType, class MV, class OP>
01541 std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::description () const
01542 {
01543   std::ostringstream out;
01544 
01545   out << "\"Belos::PseudoBlockGmresSolMgr\": {";
01546   if (this->getObjectLabel () != "") {
01547     out << "Label: " << this->getObjectLabel () << ", ";
01548   }
01549   out << "Num Blocks: " << numBlocks_
01550       << ", Maximum Iterations: " << maxIters_
01551       << ", Maximum Restarts: " << maxRestarts_
01552       << ", Convergence Tolerance: " << convtol_
01553       << "}";
01554   return out.str ();
01555 }
01556 
01557 
01558 template<class ScalarType, class MV, class OP>
01559 void
01560 PseudoBlockGmresSolMgr<ScalarType, MV, OP>::
01561 describe (Teuchos::FancyOStream &out,
01562           const Teuchos::EVerbosityLevel verbLevel) const
01563 {
01564   using Teuchos::TypeNameTraits;
01565   using Teuchos::VERB_DEFAULT;
01566   using Teuchos::VERB_NONE;
01567   using Teuchos::VERB_LOW;
01568   // using Teuchos::VERB_MEDIUM;
01569   // using Teuchos::VERB_HIGH;
01570   // using Teuchos::VERB_EXTREME;
01571   using std::endl;
01572 
01573   // Set default verbosity if applicable.
01574   const Teuchos::EVerbosityLevel vl =
01575     (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
01576 
01577   if (vl != VERB_NONE) {
01578     Teuchos::OSTab tab0 (out);
01579 
01580     out << "\"Belos::PseudoBlockGmresSolMgr\":" << endl;
01581     Teuchos::OSTab tab1 (out);
01582     out << "Template parameters:" << endl;
01583     {
01584       Teuchos::OSTab tab2 (out);
01585       out << "ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
01586           << "MV: " << TypeNameTraits<MV>::name () << endl
01587           << "OP: " << TypeNameTraits<OP>::name () << endl;
01588     }
01589     if (this->getObjectLabel () != "") {
01590       out << "Label: " << this->getObjectLabel () << endl;
01591     }
01592     out << "Num Blocks: " << numBlocks_ << endl
01593         << "Maximum Iterations: " << maxIters_ << endl
01594         << "Maximum Restarts: " << maxRestarts_ << endl
01595         << "Convergence Tolerance: " << convtol_ << endl;
01596   }
01597 }
01598 
01599 } // end Belos namespace
01600 
01601 #endif /* BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines