Belos Package Browser (Single Doxygen Collection) Development
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 #include "BelosStatusTestMaxIters.hpp"
00060 #include "BelosStatusTestGenResNorm.hpp"
00061 #include "BelosStatusTestImpResNorm.hpp"
00062 #include "BelosStatusTestCombo.hpp"
00063 #include "BelosStatusTestOutputFactory.hpp"
00064 #include "BelosOutputManager.hpp"
00065 #include "Teuchos_BLAS.hpp"
00066 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00067 #include "Teuchos_TimeMonitor.hpp"
00068 #endif
00069 
00077 namespace Belos {
00078   
00080 
00081   
00088   class PseudoBlockGmresSolMgrLinearProblemFailure : public BelosError {public:
00089     PseudoBlockGmresSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00090     {}};
00091   
00098   class PseudoBlockGmresSolMgrOrthoFailure : public BelosError {public:
00099     PseudoBlockGmresSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00100     {}};
00101 
00125   template<class ScalarType, class MV, class OP>
00126   class PseudoBlockGmresSolMgr : public SolverManager<ScalarType,MV,OP> {
00127     
00128   private:
00129     typedef MultiVecTraits<ScalarType,MV> MVT;
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     PseudoBlockGmresSolMgr();
00146     
00159     PseudoBlockGmresSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00160                             const Teuchos::RCP<Teuchos::ParameterList> &pl );
00161     
00163     virtual ~PseudoBlockGmresSolMgr() {};
00165     
00167 
00168     
00169     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00170       return *problem_;
00171     }
00172     
00175     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00176    
00179     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00180  
00186     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00187       return Teuchos::tuple(timerSolve_);
00188     }
00189 
00200     MagnitudeType achievedTol() const {
00201       return achievedTol_;
00202     }
00203 
00205     int getNumIters() const {
00206       return numIters_;
00207     }
00208   
00265     bool isLOADetected() const { return loaDetected_; }
00266   
00268     
00270 
00271    
00273     void setProblem (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem) {
00274       problem_ = problem; 
00275     }
00276    
00278     void setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params);
00279 
00286     virtual void setUserConvStatusTest(
00287       const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &userConvStatusTest
00288       );
00289     
00291 
00293 
00294 
00298     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00300 
00302 
00303     
00321     ReturnType solve();
00322     
00324     
00327     
00329     std::string description() const;
00330     
00332     
00333   private:
00334 
00349     bool checkStatusTest();
00350     
00352     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00353     
00354     // Output manager.
00355     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00356     Teuchos::RCP<std::ostream> outputStream_;
00357 
00358     // Status tests.
00359     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > userConvStatusTest_;
00360     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00361     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00362     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00363     Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > impConvTest_, expConvTest_;
00364     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00365 
00366     // Orthogonalization manager.
00367     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00368 
00369      // Current parameter list.
00370     Teuchos::RCP<Teuchos::ParameterList> params_;
00371    
00372     // Default solver values.
00373     static const MagnitudeType convtol_default_;
00374     static const MagnitudeType orthoKappa_default_;
00375     static const int maxRestarts_default_;
00376     static const int maxIters_default_;
00377     static const bool showMaxResNormOnly_default_;
00378     static const int blockSize_default_;
00379     static const int numBlocks_default_;
00380     static const int verbosity_default_;
00381     static const int outputStyle_default_;
00382     static const int outputFreq_default_;
00383     static const int defQuorum_default_;
00384     static const std::string impResScale_default_; 
00385     static const std::string expResScale_default_; 
00386     static const std::string label_default_;
00387     static const std::string orthoType_default_;
00388     static const Teuchos::RCP<std::ostream> outputStream_default_;
00389 
00390     // Current solver values.
00391     MagnitudeType convtol_, orthoKappa_, achievedTol_;
00392     int maxRestarts_, maxIters_, numIters_;
00393     int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_, defQuorum_;
00394     bool showMaxResNormOnly_;
00395     std::string orthoType_; 
00396     std::string impResScale_, expResScale_;       
00397  
00398     // Timers.
00399     std::string label_;
00400     Teuchos::RCP<Teuchos::Time> timerSolve_;
00401 
00402     // Internal state variables.
00403     bool isSet_, isSTSet_, expResTest_;
00404     bool loaDetected_;
00405   };
00406   
00407   
00408 // Default solver values.
00409 template<class ScalarType, class MV, class OP>
00410 const typename PseudoBlockGmresSolMgr<ScalarType,MV,OP>::MagnitudeType PseudoBlockGmresSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00411 
00412 template<class ScalarType, class MV, class OP>
00413 const typename PseudoBlockGmresSolMgr<ScalarType,MV,OP>::MagnitudeType PseudoBlockGmresSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00414 
00415 template<class ScalarType, class MV, class OP>
00416 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::maxRestarts_default_ = 20;
00417 
00418 template<class ScalarType, class MV, class OP>
00419 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00420 
00421 template<class ScalarType, class MV, class OP>
00422 const bool PseudoBlockGmresSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00423 
00424 template<class ScalarType, class MV, class OP>
00425 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00426 
00427 template<class ScalarType, class MV, class OP>
00428 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::numBlocks_default_ = 300;
00429 
00430 template<class ScalarType, class MV, class OP>
00431 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00432 
00433 template<class ScalarType, class MV, class OP>
00434 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00435 
00436 template<class ScalarType, class MV, class OP>
00437 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00438 
00439 template<class ScalarType, class MV, class OP>
00440 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::defQuorum_default_ = 1;
00441 
00442 template<class ScalarType, class MV, class OP>
00443 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of Preconditioned Initial Residual";
00444 
00445 template<class ScalarType, class MV, class OP>
00446 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of Initial Residual";
00447 
00448 template<class ScalarType, class MV, class OP>
00449 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00450 
00451 template<class ScalarType, class MV, class OP>
00452 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00453 
00454 template<class ScalarType, class MV, class OP>
00455 const Teuchos::RCP<std::ostream> PseudoBlockGmresSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00456 
00457 
00458 // Empty Constructor
00459 template<class ScalarType, class MV, class OP>
00460 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::PseudoBlockGmresSolMgr() :
00461   outputStream_(outputStream_default_),
00462   convtol_(convtol_default_),
00463   orthoKappa_(orthoKappa_default_),
00464   achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
00465   maxRestarts_(maxRestarts_default_),
00466   maxIters_(maxIters_default_),
00467   numIters_(0),
00468   blockSize_(blockSize_default_),
00469   numBlocks_(numBlocks_default_),
00470   verbosity_(verbosity_default_),
00471   outputStyle_(outputStyle_default_),
00472   outputFreq_(outputFreq_default_),
00473   defQuorum_(defQuorum_default_),
00474   showMaxResNormOnly_(showMaxResNormOnly_default_),
00475   orthoType_(orthoType_default_),
00476   impResScale_(impResScale_default_),
00477   expResScale_(expResScale_default_),
00478   label_(label_default_),
00479   isSet_(false),
00480   isSTSet_(false),
00481   expResTest_(false),
00482   loaDetected_(false)
00483 {}
00484 
00485 // Basic Constructor
00486 template<class ScalarType, class MV, class OP>
00487 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::
00488 PseudoBlockGmresSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00489                         const Teuchos::RCP<Teuchos::ParameterList> &pl) : 
00490   problem_(problem),
00491   outputStream_(outputStream_default_),
00492   convtol_(convtol_default_),
00493   orthoKappa_(orthoKappa_default_),
00494   achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
00495   maxRestarts_(maxRestarts_default_),
00496   maxIters_(maxIters_default_),
00497   numIters_(0),
00498   blockSize_(blockSize_default_),
00499   numBlocks_(numBlocks_default_),
00500   verbosity_(verbosity_default_),
00501   outputStyle_(outputStyle_default_),
00502   outputFreq_(outputFreq_default_),
00503   defQuorum_(defQuorum_default_),
00504   showMaxResNormOnly_(showMaxResNormOnly_default_),
00505   orthoType_(orthoType_default_),
00506   impResScale_(impResScale_default_),
00507   expResScale_(expResScale_default_),
00508   label_(label_default_),
00509   isSet_(false),
00510   isSTSet_(false),
00511   expResTest_(false),
00512   loaDetected_(false)
00513 {
00514   TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00515 
00516   // If the parameter list pointer is null, then set the current parameters to the default parameter list.
00517   if (!is_null(pl)) {
00518     // Set the parameters using the list that was passed in.
00519     setParameters( pl );  
00520   }
00521 }
00522 
00523 template<class ScalarType, class MV, class OP>
00524 void PseudoBlockGmresSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00525 {
00526   // Create the internal parameter list if ones doesn't already exist.
00527   if (params_ == Teuchos::null) {
00528     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00529   }
00530   else {
00531     params->validateParameters(*getValidParameters());
00532   }
00533 
00534   // Check for maximum number of restarts
00535   if (params->isParameter("Maximum Restarts")) {
00536     maxRestarts_ = params->get("Maximum Restarts",maxRestarts_default_);
00537 
00538     // Update parameter in our list.
00539     params_->set("Maximum Restarts", maxRestarts_);
00540   }
00541 
00542   // Check for maximum number of iterations
00543   if (params->isParameter("Maximum Iterations")) {
00544     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00545 
00546     // Update parameter in our list and in status test.
00547     params_->set("Maximum Iterations", maxIters_);
00548     if (maxIterTest_!=Teuchos::null)
00549       maxIterTest_->setMaxIters( maxIters_ );
00550   }
00551 
00552   // Check for blocksize
00553   if (params->isParameter("Block Size")) {
00554     blockSize_ = params->get("Block Size",blockSize_default_);    
00555     TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00556         "Belos::PseudoBlockGmresSolMgr: \"Block Size\" must be strictly positive.");
00557 
00558     // Update parameter in our list.
00559     params_->set("Block Size", blockSize_);
00560   }
00561 
00562   // Check for the maximum number of blocks.
00563   if (params->isParameter("Num Blocks")) {
00564     numBlocks_ = params->get("Num Blocks",numBlocks_default_);
00565     TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00566         "Belos::PseudoBlockGmresSolMgr: \"Num Blocks\" must be strictly positive.");
00567 
00568     // Update parameter in our list.
00569     params_->set("Num Blocks", numBlocks_);
00570   }
00571 
00572   // Check to see if the timer label changed.
00573   if (params->isParameter("Timer Label")) {
00574     std::string tempLabel = params->get("Timer Label", label_default_);
00575 
00576     // Update parameter in our list and solver timer
00577     if (tempLabel != label_) {
00578       label_ = tempLabel;
00579       params_->set("Timer Label", label_);
00580       std::string solveLabel = label_ + ": PseudoBlockGmresSolMgr total solve time";
00581 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00582       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00583 #endif
00584     }
00585   }
00586 
00587   // Check if the orthogonalization changed.
00588   if (params->isParameter("Orthogonalization")) {
00589     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00590     TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00591                         std::invalid_argument,
00592                         "Belos::PseudoBlockGmresSolMgr: \"Orthogonalization\" must be either \"DGKS\",\"ICGS\", or \"IMGS\".");
00593     if (tempOrthoType != orthoType_) {
00594       orthoType_ = tempOrthoType;
00595       // Create orthogonalization manager
00596       if (orthoType_=="DGKS") {
00597         if (orthoKappa_ <= 0) {
00598           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00599         }
00600         else {
00601           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00602           Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00603         }
00604       }
00605       else if (orthoType_=="ICGS") {
00606         ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00607       } 
00608       else if (orthoType_=="IMGS") {
00609         ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00610       } 
00611     }  
00612   }
00613 
00614   // Check which orthogonalization constant to use.
00615   if (params->isParameter("Orthogonalization Constant")) {
00616     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00617 
00618     // Update parameter in our list.
00619     params_->set("Orthogonalization Constant",orthoKappa_);
00620     if (orthoType_=="DGKS") {
00621       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00622         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00623       }
00624     } 
00625   }
00626 
00627   // Check for a change in verbosity level
00628   if (params->isParameter("Verbosity")) {
00629     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00630       verbosity_ = params->get("Verbosity", verbosity_default_);
00631     } else {
00632       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00633     }
00634 
00635     // Update parameter in our list.
00636     params_->set("Verbosity", verbosity_);
00637     if (printer_ != Teuchos::null)
00638       printer_->setVerbosity(verbosity_);
00639   }
00640 
00641   // Check for a change in output style.
00642   if (params->isParameter("Output Style")) {
00643     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00644       outputStyle_ = params->get("Output Style", outputStyle_default_);
00645     } else {
00646       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00647     }
00648 
00649     // Update parameter in our list.
00650     params_->set("Output Style", verbosity_);
00651     if (outputTest_ != Teuchos::null)
00652       isSTSet_ = false;
00653   }
00654 
00655   // output stream
00656   if (params->isParameter("Output Stream")) {
00657     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00658 
00659     // Update parameter in our list.
00660     params_->set("Output Stream", outputStream_);
00661     if (printer_ != Teuchos::null)
00662       printer_->setOStream( outputStream_ );
00663   }
00664 
00665   // frequency level
00666   if (verbosity_ & Belos::StatusTestDetails) {
00667     if (params->isParameter("Output Frequency")) {
00668       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00669     }
00670 
00671     // Update parameter in out list and output status test.
00672     params_->set("Output Frequency", outputFreq_);
00673     if (outputTest_ != Teuchos::null)
00674       outputTest_->setOutputFrequency( outputFreq_ );
00675   }
00676 
00677   // Create output manager if we need to.
00678   if (printer_ == Teuchos::null) {
00679     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00680   }
00681   
00682   // Convergence
00683   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00684   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00685 
00686   // Check for convergence tolerance
00687   if (params->isParameter("Convergence Tolerance")) {
00688     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00689 
00690     // Update parameter in our list and residual tests.
00691     params_->set("Convergence Tolerance", convtol_);
00692     if (impConvTest_ != Teuchos::null)
00693       impConvTest_->setTolerance( convtol_ );
00694     if (expConvTest_ != Teuchos::null)
00695       expConvTest_->setTolerance( convtol_ );
00696   }
00697 
00698   // Check for a change in scaling, if so we need to build new residual tests.
00699   bool newImpResTest = false, newExpResTest = false; 
00700   if (params->isParameter("Implicit Residual Scaling")) {
00701     std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
00702 
00703     // Only update the scaling if it's different.
00704     if (impResScale_ != tempImpResScale) {
00705       Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
00706       impResScale_ = tempImpResScale;
00707 
00708       // Update parameter in our list and residual tests
00709       params_->set("Implicit Residual Scaling", impResScale_);
00710       if (impConvTest_ != Teuchos::null) {
00711         try {
00712           impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
00713         }
00714         catch (std::exception& e) {
00715           // Make sure the convergence test gets constructed again.
00716           isSTSet_ = false;
00717           newImpResTest = true;
00718         }
00719       }
00720     }      
00721   }
00722 
00723   if (params->isParameter("Explicit Residual Scaling")) {
00724     std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
00725 
00726     // Only update the scaling if it's different.
00727     if (expResScale_ != tempExpResScale) {
00728       Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
00729       expResScale_ = tempExpResScale;
00730 
00731       // Update parameter in our list and residual tests
00732       params_->set("Explicit Residual Scaling", expResScale_);
00733       if (expConvTest_ != Teuchos::null) {
00734         try {
00735           expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
00736         }
00737         catch (std::exception& e) {
00738           // Make sure the convergence test gets constructed again.
00739           isSTSet_ = false;
00740           newExpResTest = true;
00741         }
00742       }
00743     }      
00744   }
00745 
00746 
00747   if (params->isParameter("Show Maximum Residual Norm Only")) {
00748     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00749 
00750     // Update parameter in our list and residual tests
00751     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00752     if (impConvTest_ != Teuchos::null)
00753       impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00754     if (expConvTest_ != Teuchos::null)
00755       expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00756   }
00757 
00758   // Create status tests if we need to.
00759 
00760   // Get the deflation quorum, or number of converged systems before deflation is allowed
00761   if (params->isParameter("Deflation Quorum")) {
00762     defQuorum_ = params->get("Deflation Quorum", defQuorum_);
00763     TEUCHOS_TEST_FOR_EXCEPTION(defQuorum_ > blockSize_, std::invalid_argument,
00764         "Belos::PseudoBlockGmresSolMgr: \"Deflation Quorum\" cannot be larger than \"Block Size\".");
00765     params_->set("Deflation Quorum", defQuorum_);
00766     if (impConvTest_ != Teuchos::null)
00767       impConvTest_->setQuorum( defQuorum_ );
00768     if (expConvTest_ != Teuchos::null)
00769       expConvTest_->setQuorum( defQuorum_ );
00770   }
00771 
00772   // Create orthogonalization manager if we need to.
00773   if (ortho_ == Teuchos::null) {
00774     if (orthoType_=="DGKS") {
00775       if (orthoKappa_ <= 0) {
00776         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00777       }
00778       else {
00779         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00780         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00781       }
00782     }
00783     else if (orthoType_=="ICGS") {
00784       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00785     } 
00786     else if (orthoType_=="IMGS") {
00787       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00788     } 
00789     else {
00790       TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00791           "Belos::PseudoBlockGmresSolMgr(): Invalid orthogonalization type.");
00792     }  
00793   }
00794 
00795   // Create the timer if we need to.
00796   if (timerSolve_ == Teuchos::null) {
00797     std::string solveLabel = label_ + ": PseudoBlockGmresSolMgr total solve time";
00798 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00799     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00800 #endif
00801   }
00802 
00803   // Inform the solver manager that the current parameters were set.
00804   isSet_ = true;
00805 }
00806 
00807     
00808 template<class ScalarType, class MV, class OP>
00809 void 
00810 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::setUserConvStatusTest(
00811   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &userConvStatusTest
00812   )
00813 {
00814   userConvStatusTest_ = userConvStatusTest;
00815 }
00816 
00817     
00818 template<class ScalarType, class MV, class OP>
00819 Teuchos::RCP<const Teuchos::ParameterList>
00820 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::getValidParameters() const
00821 {
00822   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00823   if (is_null(validPL)) {
00824     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00825   // Set all the valid parameters and their default values.
00826     pl= Teuchos::rcp( new Teuchos::ParameterList() );
00827     pl->set("Convergence Tolerance", convtol_default_,
00828       "The relative residual tolerance that needs to be achieved by the\n"
00829       "iterative solver in order for the linear system to be declared converged.");
00830     pl->set("Maximum Restarts", maxRestarts_default_,
00831       "The maximum number of restarts allowed for each\n"
00832       "set of RHS solved.");
00833     pl->set("Maximum Iterations", maxIters_default_,
00834       "The maximum number of block iterations allowed for each\n"
00835       "set of RHS solved.");
00836     pl->set("Num Blocks", numBlocks_default_,
00837       "The maximum number of vectors allowed in the Krylov subspace\n"
00838       "for each set of RHS solved.");
00839     pl->set("Block Size", blockSize_default_,
00840       "The number of RHS solved simultaneously.");
00841     pl->set("Verbosity", verbosity_default_,
00842       "What type(s) of solver information should be outputted\n"
00843       "to the output stream.");
00844     pl->set("Output Style", outputStyle_default_,
00845       "What style is used for the solver information outputted\n"
00846       "to the output stream.");
00847     pl->set("Output Frequency", outputFreq_default_,
00848       "How often convergence information should be outputted\n"
00849       "to the output stream.");  
00850     pl->set("Deflation Quorum", defQuorum_default_,
00851       "The number of linear systems that need to converge before\n"
00852       "they are deflated.  This number should be <= block size.");
00853     pl->set("Output Stream", outputStream_default_,
00854       "A reference-counted pointer to the output stream where all\n"
00855       "solver output is sent.");
00856     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00857       "When convergence information is printed, only show the maximum\n"
00858       "relative residual norm when the block size is greater than one.");
00859     pl->set("Implicit Residual Scaling", impResScale_default_,
00860       "The type of scaling used in the implicit residual convergence test.");
00861     pl->set("Explicit Residual Scaling", expResScale_default_,
00862       "The type of scaling used in the explicit residual convergence test.");
00863     pl->set("Timer Label", label_default_,
00864       "The string to use as a prefix for the timer labels.");
00865     //  defaultParams_->set("Restart Timers", restartTimers_);
00866     pl->set("Orthogonalization", orthoType_default_,
00867       "The type of orthogonalization to use: DGKS, ICGS, IMGS.");
00868     pl->set("Orthogonalization Constant",orthoKappa_default_,
00869       "The constant used by DGKS orthogonalization to determine\n"
00870       "whether another step of classical Gram-Schmidt is necessary.");
00871     validPL = pl;
00872   }
00873   return validPL;
00874 }
00875 
00876 // Check the status test versus the defined linear problem
00877 template<class ScalarType, class MV, class OP>
00878 bool PseudoBlockGmresSolMgr<ScalarType,MV,OP>::checkStatusTest() {
00879 
00880   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00881   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestGenResNorm_t;
00882   typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP>  StatusTestImpResNorm_t;
00883 
00884   // Basic test checks maximum iterations and native residual.
00885   maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00886 
00887   // If there is a left preconditioner, we create a combined status test that checks the implicit
00888   // and then explicit residual norm to see if we have convergence.
00889   if ( !Teuchos::is_null(problem_->getLeftPrec()) ) {
00890     expResTest_ = true;
00891   }
00892 
00893   if (expResTest_) {
00894    
00895     // Implicit residual test, using the native residual to determine if convergence was achieved.
00896     Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
00897       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
00898     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00899     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00900     impConvTest_ = tmpImpConvTest;
00901 
00902     // Explicit residual test once the native residual is below the tolerance
00903     Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
00904       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
00905     tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
00906     tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
00907     tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00908     expConvTest_ = tmpExpConvTest;
00909 
00910     // The convergence test is a combination of the "cheap" implicit test and explicit test.
00911     convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
00912   }
00913   else {
00914 
00915     // Implicit residual test, using the native residual to determine if convergence was achieved.
00916     // Use test that checks for loss of accuracy.
00917     Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
00918       Teuchos::rcp( new StatusTestImpResNorm_t( convtol_, defQuorum_ ) );
00919     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00920     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00921     impConvTest_ = tmpImpConvTest;
00922     
00923     // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
00924     expConvTest_ = impConvTest_;
00925     convTest_ = impConvTest_;
00926   }
00927 
00928   if (nonnull(userConvStatusTest_) ) {
00929     // Override the overall convergence test with the users convergence test
00930     convTest_ = Teuchos::rcp(
00931       new StatusTestCombo_t( StatusTestCombo_t::SEQ, convTest_, userConvStatusTest_ ) );
00932     // NOTE: Above, you have to run the other convergence tests also because
00933     // the logic in this class depends on it.  This is very unfortunate.
00934   }
00935 
00936   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00937 
00938   // Create the status test output class.
00939   // This class manages and formats the output from the status test.
00940   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00941   outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00942   
00943   // Set the solver string for the output test
00944   std::string solverDesc = " Pseudo Block Gmres ";
00945   outputTest_->setSolverDesc( solverDesc );
00946 
00947 
00948   // The status test is now set.
00949   isSTSet_ = true;
00950 
00951   return false;
00952 }
00953 
00954 
00955 // solve()
00956 template<class ScalarType, class MV, class OP>
00957 ReturnType PseudoBlockGmresSolMgr<ScalarType,MV,OP>::solve() {
00958 
00959   // Set the current parameters if they were not set before.
00960   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00961   // then didn't set any parameters using setParameters().
00962   if (!isSet_) { setParameters( params_ ); }
00963   
00964   Teuchos::BLAS<int,ScalarType> blas;
00965   
00966   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockGmresSolMgrLinearProblemFailure,
00967                      "Belos::PseudoBlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00968 
00969   // Check if we have to create the status tests.
00970   if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
00971     TEUCHOS_TEST_FOR_EXCEPTION( checkStatusTest(),PseudoBlockGmresSolMgrLinearProblemFailure,
00972       "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
00973   }
00974 
00975   // Create indices for the linear systems to be solved.
00976   int startPtr = 0;
00977   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00978   int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
00979 
00980   std::vector<int> currIdx( numCurrRHS );
00981   blockSize_ = numCurrRHS;
00982   for (int i=0; i<numCurrRHS; ++i) 
00983     { currIdx[i] = startPtr+i; }
00984 
00985   // Inform the linear problem of the current linear system to solve.
00986   problem_->setLSIndex( currIdx );
00987 
00989   // Parameter list
00990   Teuchos::ParameterList plist;
00991   plist.set("Num Blocks",numBlocks_);
00992   
00993   // Reset the status test.  
00994   outputTest_->reset();
00995   loaDetected_ = false;
00996 
00997   // Assume convergence is achieved, then let any failed convergence set this to false.
00998   bool isConverged = true;
00999 
01001   // BlockGmres solver
01002 
01003   Teuchos::RCP<PseudoBlockGmresIter<ScalarType,MV,OP> > block_gmres_iter
01004     = Teuchos::rcp( new PseudoBlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );  
01005 
01006   // Enter solve() iterations
01007   {
01008 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01009     Teuchos::TimeMonitor slvtimer(*timerSolve_);
01010 #endif
01011 
01012     while ( numRHS2Solve > 0 ) {
01013 
01014       // Reset the active / converged vectors from this block
01015       std::vector<int> convRHSIdx;
01016       std::vector<int> currRHSIdx( currIdx );
01017       currRHSIdx.resize(numCurrRHS);
01018 
01019       // Set the current number of blocks with the pseudo Gmres iteration.
01020       block_gmres_iter->setNumBlocks( numBlocks_ );
01021 
01022       // Reset the number of iterations.
01023       block_gmres_iter->resetNumIters();
01024 
01025       // Reset the number of calls that the status test output knows about.
01026       outputTest_->resetNumCalls();
01027 
01028       // Get a new state struct and initialize the solver.
01029       PseudoBlockGmresIterState<ScalarType,MV> newState;
01030 
01031       // Create the first block in the current Krylov basis for each right-hand side.
01032       std::vector<int> index(1);
01033       Teuchos::RCP<MV> tmpV, R_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
01034       newState.V.resize( blockSize_ );
01035       newState.Z.resize( blockSize_ );
01036       for (int i=0; i<blockSize_; ++i) {
01037         index[0]=i;
01038         tmpV = MVT::CloneViewNonConst( *R_0, index );
01039 
01040         // Get a matrix to hold the orthonormalization coefficients.
01041         Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
01042           = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
01043 
01044         // Orthonormalize the new V_0
01045         int rank = ortho_->normalize( *tmpV, tmpZ );
01046         TEUCHOS_TEST_FOR_EXCEPTION(rank != 1, PseudoBlockGmresSolMgrOrthoFailure,
01047             "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
01048 
01049         newState.V[i] = tmpV;
01050         newState.Z[i] = tmpZ;
01051       }
01052 
01053       newState.curDim = 0;
01054       block_gmres_iter->initialize(newState);
01055       int numRestarts = 0;
01056 
01057       while(1) {
01058         
01059         // tell block_gmres_iter to iterate
01060         try {
01061           block_gmres_iter->iterate();
01062           
01064           //
01065           // check convergence first
01066           //
01068           if ( convTest_->getStatus() == Passed ) {
01069         
01070             if ( expConvTest_->getLOADetected() ) {
01071         // Oops!  There was a loss of accuracy (LOA) for one or
01072         // more right-hand sides.  That means the implicit
01073         // (a.k.a. "native") residuals claim convergence,
01074         // whereas the explicit (hence expConvTest_, i.e.,
01075         // "explicit convergence test") residuals have not
01076         // converged.
01077         //
01078         // We choose to deal with this situation by deflating
01079         // out the affected right-hand sides and moving on.
01080               loaDetected_ = true;
01081               printer_->stream(Warnings) <<
01082                 "Belos::PseudoBlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
01083               isConverged = false;
01084             }
01085         
01086             // Figure out which linear systems converged.
01087             std::vector<int> convIdx = expConvTest_->convIndices();
01088         
01089             // If the number of converged linear systems is equal to the
01090             // number of current linear systems, then we are done with this block.
01091             if (convIdx.size() == currRHSIdx.size())
01092               break;  // break from while(1){block_gmres_iter->iterate()}
01093         
01094             // Get a new state struct and initialize the solver.
01095             PseudoBlockGmresIterState<ScalarType,MV> defState;
01096         
01097             // Inform the linear problem that we are finished with this current linear system.
01098             problem_->setCurrLS();
01099         
01100             // Get the state.
01101             PseudoBlockGmresIterState<ScalarType,MV> oldState = block_gmres_iter->getState();
01102             int curDim = oldState.curDim;
01103         
01104             // Get a new state struct and reset currRHSIdx to have the right-hand sides that 
01105             // are left to converge for this block.
01106             int have = 0;
01107             std::vector<int> oldRHSIdx( currRHSIdx );
01108             std::vector<int> defRHSIdx;
01109             for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
01110               bool found = false;
01111               for (unsigned int j=0; j<convIdx.size(); ++j) {
01112                 if (currRHSIdx[i] == convIdx[j]) {
01113                   found = true;
01114                   break;
01115                 }
01116               }
01117               if (found) {
01118                 defRHSIdx.push_back( i );
01119               }
01120               else {
01121                 defState.V.push_back( Teuchos::rcp_const_cast<MV>( oldState.V[i] ) );
01122                 defState.Z.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.Z[i] ) );
01123                 defState.H.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseMatrix<int,ScalarType> >( oldState.H[i] ) );
01124                 defState.sn.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.sn[i] ) );
01125                 defState.cs.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,MagnitudeType> >(oldState.cs[i] ) );
01126                 currRHSIdx[have] = currRHSIdx[i];
01127                 have++;
01128               }
01129             }
01130             defRHSIdx.resize(currRHSIdx.size()-have);
01131             currRHSIdx.resize(have);
01132         
01133             // Compute the current solution that needs to be deflated if this solver has taken any steps.
01134             if (curDim) {
01135               Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01136               Teuchos::RCP<MV> defUpdate = MVT::CloneViewNonConst( *update, defRHSIdx );
01137               
01138               // Set the deflated indices so we can update the solution.
01139               problem_->setLSIndex( convIdx );
01140               
01141               // Update the linear problem.
01142               problem_->updateSolution( defUpdate, true );
01143             }
01144             
01145             // Set the remaining indices after deflation.
01146             problem_->setLSIndex( currRHSIdx );
01147             
01148             // Set the dimension of the subspace, which is the same as the old subspace size.
01149             defState.curDim = curDim;
01150             
01151             // Initialize the solver with the deflated system.
01152             block_gmres_iter->initialize(defState);
01153           }
01155           //
01156           // check for maximum iterations
01157           //
01159           else if ( maxIterTest_->getStatus() == Passed ) {
01160             // we don't have convergence
01161             isConverged = false;
01162             break;  // break from while(1){block_gmres_iter->iterate()}
01163           }
01165           //
01166           // check for restarting, i.e. the subspace is full
01167           //
01169           else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
01170         
01171             if ( numRestarts >= maxRestarts_ ) {
01172               isConverged = false;
01173               break; // break from while(1){block_gmres_iter->iterate()}
01174             }
01175             numRestarts++;
01176         
01177             printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
01178             
01179             // Update the linear problem.
01180             Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01181             problem_->updateSolution( update, true );
01182             
01183             // Get the state.
01184             PseudoBlockGmresIterState<ScalarType,MV> oldState = block_gmres_iter->getState();
01185             
01186             // Set the new state.
01187             PseudoBlockGmresIterState<ScalarType,MV> newstate;
01188             newstate.V.resize(currRHSIdx.size());
01189             newstate.Z.resize(currRHSIdx.size());
01190         
01191             // Compute the restart vectors
01192             // NOTE: Force the linear problem to update the current residual since the solution was updated.
01193             R_0 = MVT::Clone( *(problem_->getInitPrecResVec()), currRHSIdx.size() );
01194             problem_->computeCurrPrecResVec( &*R_0 );
01195             for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
01196               index[0] = i;  // index(1) vector declared on line 891
01197         
01198               tmpV = MVT::CloneViewNonConst( *R_0, index );
01199         
01200               // Get a matrix to hold the orthonormalization coefficients.
01201               Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
01202                 = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
01203               
01204               // Orthonormalize the new V_0
01205               int rank = ortho_->normalize( *tmpV, tmpZ );
01206               TEUCHOS_TEST_FOR_EXCEPTION(rank != 1 ,PseudoBlockGmresSolMgrOrthoFailure,
01207                   "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after the restart.");
01208               
01209               newstate.V[i] = tmpV;
01210               newstate.Z[i] = tmpZ;
01211             }
01212         
01213             // Initialize the solver.
01214             newstate.curDim = 0;
01215             block_gmres_iter->initialize(newstate);
01216         
01217           } // end of restarting
01218         
01220           //
01221           // we returned from iterate(), but none of our status tests Passed.
01222           // something is wrong, and it is probably our fault.
01223           //
01225         
01226           else {
01227             TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
01228                 "Belos::PseudoBlockGmresSolMgr::solve(): Invalid return from PseudoBlockGmresIter::iterate().");
01229           }
01230         }
01231         catch (const PseudoBlockGmresIterOrthoFailure &e) {
01232         
01233           // Try to recover the most recent least-squares solution
01234           block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
01235         
01236           // Check to see if the most recent least-squares solution yielded convergence.
01237           sTest_->checkStatus( &*block_gmres_iter );
01238           if (convTest_->getStatus() != Passed)
01239             isConverged = false;
01240           break;
01241         }
01242         catch (const std::exception &e) {
01243           printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockGmresIter::iterate() at iteration " 
01244                                    << block_gmres_iter->getNumIters() << std::endl 
01245                                    << e.what() << std::endl;
01246           throw;
01247         }
01248       }
01249       
01250       // Compute the current solution.
01251       // Update the linear problem.
01252       if (nonnull(userConvStatusTest_)) {
01253         //std::cout << "\nnonnull(userConvStatusTest_)\n";
01254         Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01255         problem_->updateSolution( update, true );
01256       }
01257       else if (nonnull(expConvTest_->getSolution())) {
01258         //std::cout << "\nexpConvTest_->getSolution()\n";
01259         Teuchos::RCP<MV> newX = expConvTest_->getSolution();
01260         Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
01261         MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
01262       }
01263       else {
01264         //std::cout << "\nblock_gmres_iter->getCurrentUpdate()\n";
01265         Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01266         problem_->updateSolution( update, true );
01267       }
01268 
01269       // Inform the linear problem that we are finished with this block linear system.
01270       problem_->setCurrLS();
01271       
01272       // Update indices for the linear systems to be solved.
01273       startPtr += numCurrRHS;
01274       numRHS2Solve -= numCurrRHS;
01275       if ( numRHS2Solve > 0 ) {
01276         numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01277 
01278         blockSize_ = numCurrRHS;
01279         currIdx.resize( numCurrRHS  );
01280         for (int i=0; i<numCurrRHS; ++i) 
01281         { currIdx[i] = startPtr+i; }
01282 
01283         // Adapt the status test quorum if we need to.
01284         if (defQuorum_ > blockSize_) {
01285           if (impConvTest_ != Teuchos::null)
01286             impConvTest_->setQuorum( blockSize_ );
01287           if (expConvTest_ != Teuchos::null)
01288             expConvTest_->setQuorum( blockSize_ );
01289         }
01290 
01291         // Set the next indices.
01292         problem_->setLSIndex( currIdx );
01293       }
01294       else {
01295         currIdx.resize( numRHS2Solve );
01296       }
01297 
01298     }// while ( numRHS2Solve > 0 )
01299 
01300   }
01301 
01302   // print final summary
01303   sTest_->print( printer_->stream(FinalSummary) );
01304  
01305   // print timing information
01306 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01307   // Calling summarize() can be expensive, so don't call unless the
01308   // user wants to print out timing details.  summarize() will do all
01309   // the work even if it's passed a "black hole" output stream.
01310   if (verbosity_ & TimingDetails)
01311     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01312 #endif
01313  
01314   // get iteration information for this solve
01315   numIters_ = maxIterTest_->getNumIters();
01316 
01317   // Save the convergence test value ("achieved tolerance") for this
01318   // solve.  For this solver, convTest_ may either be a single
01319   // residual norm test, or a combination of two residual norm tests.
01320   // In the latter case, the master convergence test convTest_ is a
01321   // SEQ combo of the implicit resp. explicit tests.  If the implicit
01322   // test never passes, then the explicit test won't ever be executed.
01323   // This manifests as expConvTest_->getTestValue()->size() < 1.  We
01324   // deal with this case by using the values returned by
01325   // impConvTest_->getTestValue().
01326   {
01327     // We'll fetch the vector of residual norms one way or the other.
01328     const std::vector<MagnitudeType>* pTestValues = NULL;
01329     if (expResTest_) {
01330       pTestValues = expConvTest_->getTestValue();
01331       if (pTestValues == NULL || pTestValues->size() < 1) {
01332   pTestValues = impConvTest_->getTestValue();
01333       }
01334     } 
01335     else {
01336       // Only the implicit residual norm test is being used.
01337       pTestValues = impConvTest_->getTestValue();
01338     }
01339     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
01340       "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's "
01341       "getTestValue() method returned NULL.  Please report this bug to the "
01342       "Belos developers.");
01343     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
01344       "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's "
01345       "getTestValue() method returned a vector of length zero.  Please report "
01346       "this bug to the Belos developers.");
01347 
01348     // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
01349     // achieved tolerances for all vectors in the current solve(), or
01350     // just for the vectors from the last deflation?
01351     achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
01352   }
01353  
01354   if (!isConverged || loaDetected_) {
01355     return Unconverged; // return from PseudoBlockGmresSolMgr::solve() 
01356   }
01357   return Converged; // return from PseudoBlockGmresSolMgr::solve() 
01358 }
01359 
01360 //  This method requires the solver manager to return a std::string that describes itself.
01361 template<class ScalarType, class MV, class OP>
01362 std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::description() const
01363 {
01364   std::ostringstream oss;
01365   oss << "Belos::PseudoBlockGmresSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01366   oss << "{";
01367   oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" <<blockSize_;
01368   oss << ", Num Blocks="<<numBlocks_<< ", Max Restarts=" << maxRestarts_;
01369   oss << "}";
01370   return oss.str();
01371 }
01372   
01373 } // end Belos namespace
01374 
01375 #endif /* BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines