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 #include "Teuchos_TimeMonitor.hpp"
00067 
00084 namespace Belos {
00085   
00087 
00088   
00095   class PseudoBlockGmresSolMgrLinearProblemFailure : public BelosError {public:
00096     PseudoBlockGmresSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00097     {}};
00098   
00105   class PseudoBlockGmresSolMgrOrthoFailure : public BelosError {public:
00106     PseudoBlockGmresSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00107     {}};
00108   
00109   template<class ScalarType, class MV, class OP>
00110   class PseudoBlockGmresSolMgr : public SolverManager<ScalarType,MV,OP> {
00111     
00112   private:
00113     typedef MultiVecTraits<ScalarType,MV> MVT;
00114     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00115     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00116     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00117     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00118     
00119   public:
00120     
00122 
00123 
00129     PseudoBlockGmresSolMgr();
00130     
00143     PseudoBlockGmresSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00144                             const Teuchos::RCP<Teuchos::ParameterList> &pl );
00145     
00147     virtual ~PseudoBlockGmresSolMgr() {};
00149     
00151 
00152     
00153     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00154       return *problem_;
00155     }
00156     
00159     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00160    
00163     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00164  
00170     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00171       return tuple(timerSolve_);
00172     }
00173 
00175     int getNumIters() const {
00176       return numIters_;
00177     }
00178   
00182     bool isLOADetected() const { return loaDetected_; }
00183   
00185     
00187 
00188    
00190     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00191    
00193     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00194 
00196     virtual void setUserConvStatusTest(
00197       const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &userConvStatusTest
00198       );
00199     
00201 
00203 
00204 
00208     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00210 
00212 
00213     
00231     ReturnType solve();
00232     
00234     
00237     
00239     std::string description() const;
00240     
00242     
00243   private:
00244 
00245     // Method to convert std::string to enumerated type for residual.
00246     Belos::ScaleType convertStringToScaleType( const std::string& scaleType ) {
00247       if (scaleType == "Norm of Initial Residual") {
00248         return Belos::NormOfInitRes;
00249       } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00250         return Belos::NormOfPrecInitRes;
00251       } else if (scaleType == "Norm of RHS") {
00252         return Belos::NormOfRHS;
00253       } else if (scaleType == "None") {
00254         return Belos::None;
00255       } else 
00256         TEST_FOR_EXCEPTION( true ,std::logic_error,
00257             "Belos::PseudoBlockGmresSolMgr(): Invalid residual scaling type.");
00258     }
00259 
00260     // Method for checking current status test against defined linear problem.
00261     bool checkStatusTest();
00262     
00263     // Linear problem.
00264     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00265     
00266     // Output manager.
00267     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00268     Teuchos::RCP<std::ostream> outputStream_;
00269 
00270     // Status test.
00271     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > userConvStatusTest_;
00272     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00273     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00274     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00275     Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > impConvTest_, expConvTest_;
00276     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00277 
00278     // Orthogonalization manager.
00279     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00280 
00281      // Current parameter list.
00282     Teuchos::RCP<ParameterList> params_;
00283    
00284     // Default solver values.
00285     static const MagnitudeType convtol_default_;
00286     static const MagnitudeType orthoKappa_default_;
00287     static const int maxRestarts_default_;
00288     static const int maxIters_default_;
00289     static const bool showMaxResNormOnly_default_;
00290     static const int blockSize_default_;
00291     static const int numBlocks_default_;
00292     static const int verbosity_default_;
00293     static const int outputStyle_default_;
00294     static const int outputFreq_default_;
00295     static const int defQuorum_default_;
00296     static const std::string impResScale_default_; 
00297     static const std::string expResScale_default_; 
00298     static const std::string label_default_;
00299     static const std::string orthoType_default_;
00300     static const Teuchos::RCP<std::ostream> outputStream_default_;
00301 
00302     // Current solver values.
00303     MagnitudeType convtol_, orthoKappa_;
00304     int maxRestarts_, maxIters_, numIters_;
00305     int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_, defQuorum_;
00306     bool showMaxResNormOnly_;
00307     std::string orthoType_; 
00308     std::string impResScale_, expResScale_;       
00309  
00310     // Timers.
00311     std::string label_;
00312     Teuchos::RCP<Teuchos::Time> timerSolve_;
00313 
00314     // Internal state variables.
00315     bool isSet_, isSTSet_, expResTest_;
00316     bool loaDetected_;
00317   };
00318   
00319   
00320 // Default solver values.
00321 template<class ScalarType, class MV, class OP>
00322 const typename PseudoBlockGmresSolMgr<ScalarType,MV,OP>::MagnitudeType PseudoBlockGmresSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00323 
00324 template<class ScalarType, class MV, class OP>
00325 const typename PseudoBlockGmresSolMgr<ScalarType,MV,OP>::MagnitudeType PseudoBlockGmresSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00326 
00327 template<class ScalarType, class MV, class OP>
00328 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::maxRestarts_default_ = 20;
00329 
00330 template<class ScalarType, class MV, class OP>
00331 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00332 
00333 template<class ScalarType, class MV, class OP>
00334 const bool PseudoBlockGmresSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00335 
00336 template<class ScalarType, class MV, class OP>
00337 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00338 
00339 template<class ScalarType, class MV, class OP>
00340 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::numBlocks_default_ = 300;
00341 
00342 template<class ScalarType, class MV, class OP>
00343 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00344 
00345 template<class ScalarType, class MV, class OP>
00346 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00347 
00348 template<class ScalarType, class MV, class OP>
00349 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00350 
00351 template<class ScalarType, class MV, class OP>
00352 const int PseudoBlockGmresSolMgr<ScalarType,MV,OP>::defQuorum_default_ = 1;
00353 
00354 template<class ScalarType, class MV, class OP>
00355 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of Preconditioned Initial Residual";
00356 
00357 template<class ScalarType, class MV, class OP>
00358 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of Initial Residual";
00359 
00360 template<class ScalarType, class MV, class OP>
00361 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00362 
00363 template<class ScalarType, class MV, class OP>
00364 const std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00365 
00366 template<class ScalarType, class MV, class OP>
00367 const Teuchos::RCP<std::ostream> PseudoBlockGmresSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00368 
00369 
00370 // Empty Constructor
00371 template<class ScalarType, class MV, class OP>
00372 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::PseudoBlockGmresSolMgr() :
00373   outputStream_(outputStream_default_),
00374   convtol_(convtol_default_),
00375   orthoKappa_(orthoKappa_default_),
00376   maxRestarts_(maxRestarts_default_),
00377   maxIters_(maxIters_default_),
00378   blockSize_(blockSize_default_),
00379   numBlocks_(numBlocks_default_),
00380   verbosity_(verbosity_default_),
00381   outputStyle_(outputStyle_default_),
00382   outputFreq_(outputFreq_default_),
00383   defQuorum_(defQuorum_default_),
00384   showMaxResNormOnly_(showMaxResNormOnly_default_),
00385   orthoType_(orthoType_default_),
00386   impResScale_(impResScale_default_),
00387   expResScale_(expResScale_default_),
00388   label_(label_default_),
00389   isSet_(false),
00390   isSTSet_(false),
00391   expResTest_(false),
00392   loaDetected_(false)
00393 {}
00394 
00395 // Basic Constructor
00396 template<class ScalarType, class MV, class OP>
00397 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::PseudoBlockGmresSolMgr( 
00398                                         const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00399                                         const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00400   problem_(problem),
00401   outputStream_(outputStream_default_),
00402   convtol_(convtol_default_),
00403   orthoKappa_(orthoKappa_default_),
00404   maxRestarts_(maxRestarts_default_),
00405   maxIters_(maxIters_default_),
00406   blockSize_(blockSize_default_),
00407   numBlocks_(numBlocks_default_),
00408   verbosity_(verbosity_default_),
00409   outputStyle_(outputStyle_default_),
00410   outputFreq_(outputFreq_default_),
00411   defQuorum_(defQuorum_default_),
00412   showMaxResNormOnly_(showMaxResNormOnly_default_),
00413   orthoType_(orthoType_default_),
00414   impResScale_(impResScale_default_),
00415   expResScale_(expResScale_default_),
00416   label_(label_default_),
00417   isSet_(false),
00418   isSTSet_(false),
00419   expResTest_(false),
00420   loaDetected_(false)
00421 {
00422   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00423 
00424   // If the parameter list pointer is null, then set the current parameters to the default parameter list.
00425   if (!is_null(pl)) {
00426     // Set the parameters using the list that was passed in.
00427     setParameters( pl );  
00428   }
00429 }
00430 
00431 template<class ScalarType, class MV, class OP>
00432 void PseudoBlockGmresSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00433 {
00434   // Create the internal parameter list if ones doesn't already exist.
00435   if (params_ == Teuchos::null) {
00436     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00437   }
00438   else {
00439     params->validateParameters(*getValidParameters());
00440   }
00441 
00442   // Check for maximum number of restarts
00443   if (params->isParameter("Maximum Restarts")) {
00444     maxRestarts_ = params->get("Maximum Restarts",maxRestarts_default_);
00445 
00446     // Update parameter in our list.
00447     params_->set("Maximum Restarts", maxRestarts_);
00448   }
00449 
00450   // Check for maximum number of iterations
00451   if (params->isParameter("Maximum Iterations")) {
00452     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00453 
00454     // Update parameter in our list and in status test.
00455     params_->set("Maximum Iterations", maxIters_);
00456     if (maxIterTest_!=Teuchos::null)
00457       maxIterTest_->setMaxIters( maxIters_ );
00458   }
00459 
00460   // Check for blocksize
00461   if (params->isParameter("Block Size")) {
00462     blockSize_ = params->get("Block Size",blockSize_default_);    
00463     TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00464         "Belos::PseudoBlockGmresSolMgr: \"Block Size\" must be strictly positive.");
00465 
00466     // Update parameter in our list.
00467     params_->set("Block Size", blockSize_);
00468   }
00469 
00470   // Check for the maximum number of blocks.
00471   if (params->isParameter("Num Blocks")) {
00472     numBlocks_ = params->get("Num Blocks",numBlocks_default_);
00473     TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00474         "Belos::PseudoBlockGmresSolMgr: \"Num Blocks\" must be strictly positive.");
00475 
00476     // Update parameter in our list.
00477     params_->set("Num Blocks", numBlocks_);
00478   }
00479 
00480   // Check to see if the timer label changed.
00481   if (params->isParameter("Timer Label")) {
00482     std::string tempLabel = params->get("Timer Label", label_default_);
00483 
00484     // Update parameter in our list and solver timer
00485     if (tempLabel != label_) {
00486       label_ = tempLabel;
00487       params_->set("Timer Label", label_);
00488       std::string solveLabel = label_ + ": PseudoBlockGmresSolMgr total solve time";
00489       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00490     }
00491   }
00492 
00493   // Check if the orthogonalization changed.
00494   if (params->isParameter("Orthogonalization")) {
00495     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00496     TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00497                         std::invalid_argument,
00498                         "Belos::PseudoBlockGmresSolMgr: \"Orthogonalization\" must be either \"DGKS\",\"ICGS\", or \"IMGS\".");
00499     if (tempOrthoType != orthoType_) {
00500       orthoType_ = tempOrthoType;
00501       // Create orthogonalization manager
00502       if (orthoType_=="DGKS") {
00503         if (orthoKappa_ <= 0) {
00504           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00505         }
00506         else {
00507           ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00508           Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00509         }
00510       }
00511       else if (orthoType_=="ICGS") {
00512         ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00513       } 
00514       else if (orthoType_=="IMGS") {
00515         ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00516       } 
00517     }  
00518   }
00519 
00520   // Check which orthogonalization constant to use.
00521   if (params->isParameter("Orthogonalization Constant")) {
00522     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00523 
00524     // Update parameter in our list.
00525     params_->set("Orthogonalization Constant",orthoKappa_);
00526     if (orthoType_=="DGKS") {
00527       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00528         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00529       }
00530     } 
00531   }
00532 
00533   // Check for a change in verbosity level
00534   if (params->isParameter("Verbosity")) {
00535     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00536       verbosity_ = params->get("Verbosity", verbosity_default_);
00537     } else {
00538       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00539     }
00540 
00541     // Update parameter in our list.
00542     params_->set("Verbosity", verbosity_);
00543     if (printer_ != Teuchos::null)
00544       printer_->setVerbosity(verbosity_);
00545   }
00546 
00547   // Check for a change in output style.
00548   if (params->isParameter("Output Style")) {
00549     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00550       outputStyle_ = params->get("Output Style", outputStyle_default_);
00551     } else {
00552       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00553     }
00554 
00555     // Update parameter in our list.
00556     params_->set("Output Style", verbosity_);
00557     if (outputTest_ != Teuchos::null)
00558       isSTSet_ = false;
00559   }
00560 
00561   // output stream
00562   if (params->isParameter("Output Stream")) {
00563     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00564 
00565     // Update parameter in our list.
00566     params_->set("Output Stream", outputStream_);
00567     if (printer_ != Teuchos::null)
00568       printer_->setOStream( outputStream_ );
00569   }
00570 
00571   // frequency level
00572   if (verbosity_ & Belos::StatusTestDetails) {
00573     if (params->isParameter("Output Frequency")) {
00574       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00575     }
00576 
00577     // Update parameter in out list and output status test.
00578     params_->set("Output Frequency", outputFreq_);
00579     if (outputTest_ != Teuchos::null)
00580       outputTest_->setOutputFrequency( outputFreq_ );
00581   }
00582 
00583   // Create output manager if we need to.
00584   if (printer_ == Teuchos::null) {
00585     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00586   }
00587   
00588   // Convergence
00589   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00590   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00591 
00592   // Check for convergence tolerance
00593   if (params->isParameter("Convergence Tolerance")) {
00594     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00595 
00596     // Update parameter in our list and residual tests.
00597     params_->set("Convergence Tolerance", convtol_);
00598     if (impConvTest_ != Teuchos::null)
00599       impConvTest_->setTolerance( convtol_ );
00600     if (expConvTest_ != Teuchos::null)
00601       expConvTest_->setTolerance( convtol_ );
00602   }
00603 
00604   // Check for a change in scaling, if so we need to build new residual tests.
00605   bool newImpResTest = false, newExpResTest = false; 
00606   if (params->isParameter("Implicit Residual Scaling")) {
00607     std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
00608 
00609     // Only update the scaling if it's different.
00610     if (impResScale_ != tempImpResScale) {
00611       Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
00612       impResScale_ = tempImpResScale;
00613 
00614       // Update parameter in our list and residual tests
00615       params_->set("Implicit Residual Scaling", impResScale_);
00616       if (impConvTest_ != Teuchos::null) {
00617         try {
00618           impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
00619         }
00620         catch (std::exception& e) {
00621           // Make sure the convergence test gets constructed again.
00622           isSTSet_ = false;
00623           newImpResTest = true;
00624         }
00625       }
00626     }      
00627   }
00628 
00629   if (params->isParameter("Explicit Residual Scaling")) {
00630     std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
00631 
00632     // Only update the scaling if it's different.
00633     if (expResScale_ != tempExpResScale) {
00634       Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
00635       expResScale_ = tempExpResScale;
00636 
00637       // Update parameter in our list and residual tests
00638       params_->set("Explicit Residual Scaling", expResScale_);
00639       if (expConvTest_ != Teuchos::null) {
00640         try {
00641           expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
00642         }
00643         catch (std::exception& e) {
00644           // Make sure the convergence test gets constructed again.
00645           isSTSet_ = false;
00646           newExpResTest = true;
00647         }
00648       }
00649     }      
00650   }
00651 
00652 
00653   if (params->isParameter("Show Maximum Residual Norm Only")) {
00654     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00655 
00656     // Update parameter in our list and residual tests
00657     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00658     if (impConvTest_ != Teuchos::null)
00659       impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00660     if (expConvTest_ != Teuchos::null)
00661       expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00662   }
00663 
00664   // Create status tests if we need to.
00665 
00666   // Get the deflation quorum, or number of converged systems before deflation is allowed
00667   if (params->isParameter("Deflation Quorum")) {
00668     defQuorum_ = params->get("Deflation Quorum", defQuorum_);
00669     TEST_FOR_EXCEPTION(defQuorum_ > blockSize_, std::invalid_argument,
00670         "Belos::PseudoBlockGmresSolMgr: \"Deflation Quorum\" cannot be larger than \"Block Size\".");
00671     params_->set("Deflation Quorum", defQuorum_);
00672     if (impConvTest_ != Teuchos::null)
00673       impConvTest_->setQuorum( defQuorum_ );
00674     if (expConvTest_ != Teuchos::null)
00675       expConvTest_->setQuorum( defQuorum_ );
00676   }
00677 
00678   // Create orthogonalization manager if we need to.
00679   if (ortho_ == Teuchos::null) {
00680     if (orthoType_=="DGKS") {
00681       if (orthoKappa_ <= 0) {
00682         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00683       }
00684       else {
00685         ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00686         Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00687       }
00688     }
00689     else if (orthoType_=="ICGS") {
00690       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00691     } 
00692     else if (orthoType_=="IMGS") {
00693       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00694     } 
00695     else {
00696       TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00697           "Belos::PseudoBlockGmresSolMgr(): Invalid orthogonalization type.");
00698     }  
00699   }
00700 
00701   // Create the timer if we need to.
00702   if (timerSolve_ == Teuchos::null) {
00703     std::string solveLabel = label_ + ": PseudoBlockGmresSolMgr total solve time";
00704     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00705   }
00706 
00707   // Inform the solver manager that the current parameters were set.
00708   isSet_ = true;
00709 }
00710 
00711     
00712 template<class ScalarType, class MV, class OP>
00713 void 
00714 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::setUserConvStatusTest(
00715   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &userConvStatusTest
00716   )
00717 {
00718   userConvStatusTest_ = userConvStatusTest;
00719 }
00720 
00721     
00722 template<class ScalarType, class MV, class OP>
00723 Teuchos::RCP<const Teuchos::ParameterList>
00724 PseudoBlockGmresSolMgr<ScalarType,MV,OP>::getValidParameters() const
00725 {
00726   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00727   if (is_null(validPL)) {
00728     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00729   // Set all the valid parameters and their default values.
00730     pl= Teuchos::rcp( new Teuchos::ParameterList() );
00731     pl->set("Convergence Tolerance", convtol_default_,
00732       "The relative residual tolerance that needs to be achieved by the\n"
00733       "iterative solver in order for the linera system to be declared converged.");
00734     pl->set("Maximum Restarts", maxRestarts_default_,
00735       "The maximum number of restarts allowed for each\n"
00736       "set of RHS solved.");
00737     pl->set("Maximum Iterations", maxIters_default_,
00738       "The maximum number of block iterations allowed for each\n"
00739       "set of RHS solved.");
00740     pl->set("Num Blocks", numBlocks_default_,
00741       "The maximum number of vectors allowed in the Krylov subspace\n"
00742       "for each set of RHS solved.");
00743     pl->set("Block Size", blockSize_default_,
00744       "The number of RHS solved simultaneously.");
00745     pl->set("Verbosity", verbosity_default_,
00746       "What type(s) of solver information should be outputted\n"
00747       "to the output stream.");
00748     pl->set("Output Style", outputStyle_default_,
00749       "What style is used for the solver information outputted\n"
00750       "to the output stream.");
00751     pl->set("Output Frequency", outputFreq_default_,
00752       "How often convergence information should be outputted\n"
00753       "to the output stream.");  
00754     pl->set("Deflation Quorum", defQuorum_default_,
00755       "The number of linear systems that need to converge before\n"
00756       "they are deflated.  This number should be <= block size.");
00757     pl->set("Output Stream", outputStream_default_,
00758       "A reference-counted pointer to the output stream where all\n"
00759       "solver output is sent.");
00760     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00761       "When convergence information is printed, only show the maximum\n"
00762       "relative residual norm when the block size is greater than one.");
00763     pl->set("Implicit Residual Scaling", impResScale_default_,
00764       "The type of scaling used in the implicit residual convergence test.");
00765     pl->set("Explicit Residual Scaling", expResScale_default_,
00766       "The type of scaling used in the explicit residual convergence test.");
00767     pl->set("Timer Label", label_default_,
00768       "The string to use as a prefix for the timer labels.");
00769     //  defaultParams_->set("Restart Timers", restartTimers_);
00770     pl->set("Orthogonalization", orthoType_default_,
00771       "The type of orthogonalization to use: DGKS, ICGS, IMGS.");
00772     pl->set("Orthogonalization Constant",orthoKappa_default_,
00773       "The constant used by DGKS orthogonalization to determine\n"
00774       "whether another step of classical Gram-Schmidt is necessary.");
00775     validPL = pl;
00776   }
00777   return validPL;
00778 }
00779 
00780 // Check the status test versus the defined linear problem
00781 template<class ScalarType, class MV, class OP>
00782 bool PseudoBlockGmresSolMgr<ScalarType,MV,OP>::checkStatusTest() {
00783 
00784   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00785   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestGenResNorm_t;
00786   typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP>  StatusTestImpResNorm_t;
00787 
00788   // Basic test checks maximum iterations and native residual.
00789   maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00790 
00791   // If there is a left preconditioner, we create a combined status test that checks the implicit
00792   // and then explicit residual norm to see if we have convergence.
00793   if ( !Teuchos::is_null(problem_->getLeftPrec()) ) {
00794     expResTest_ = true;
00795   }
00796 
00797   if (expResTest_) {
00798    
00799     // Implicit residual test, using the native residual to determine if convergence was achieved.
00800     Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
00801       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
00802     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00803     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00804     impConvTest_ = tmpImpConvTest;
00805 
00806     // Explicit residual test once the native residual is below the tolerance
00807     Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
00808       Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
00809     tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
00810     tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
00811     tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00812     expConvTest_ = tmpExpConvTest;
00813 
00814     // The convergence test is a combination of the "cheap" implicit test and explicit test.
00815     convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
00816   }
00817   else {
00818 
00819     // Implicit residual test, using the native residual to determine if convergence was achieved.
00820     // Use test that checks for loss of accuracy.
00821     Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
00822       Teuchos::rcp( new StatusTestImpResNorm_t( convtol_, defQuorum_ ) );
00823     tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00824     tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
00825     impConvTest_ = tmpImpConvTest;
00826     
00827     // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
00828     expConvTest_ = impConvTest_;
00829     convTest_ = impConvTest_;
00830   }
00831 
00832   if (nonnull(userConvStatusTest_) ) {
00833     // Override the overall convergence test with the users convergence test
00834     convTest_ = Teuchos::rcp(
00835       new StatusTestCombo_t( StatusTestCombo_t::SEQ, convTest_, userConvStatusTest_ ) );
00836     // NOTE: Above, you have to run the other convergence tests also because
00837     // the logic in this class depends on it.  This is very unfortunate.
00838   }
00839 
00840   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00841 
00842   // Create the status test output class.
00843   // This class manages and formats the output from the status test.
00844   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00845   outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00846   
00847   // Set the solver string for the output test
00848   std::string solverDesc = " Pseudo Block Gmres ";
00849   outputTest_->setSolverDesc( solverDesc );
00850 
00851 
00852   // The status test is now set.
00853   isSTSet_ = true;
00854 
00855   return false;
00856 }
00857 
00858 
00859 // solve()
00860 template<class ScalarType, class MV, class OP>
00861 ReturnType PseudoBlockGmresSolMgr<ScalarType,MV,OP>::solve() {
00862 
00863   // Set the current parameters if they were not set before.
00864   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00865   // then didn't set any parameters using setParameters().
00866   if (!isSet_) { setParameters( params_ ); }
00867   
00868   Teuchos::BLAS<int,ScalarType> blas;
00869   
00870   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockGmresSolMgrLinearProblemFailure,
00871                      "Belos::PseudoBlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00872 
00873   // Check if we have to create the status tests.
00874   if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
00875     TEST_FOR_EXCEPTION( checkStatusTest(),PseudoBlockGmresSolMgrLinearProblemFailure,
00876       "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
00877   }
00878 
00879   // Create indices for the linear systems to be solved.
00880   int startPtr = 0;
00881   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00882   int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
00883 
00884   std::vector<int> currIdx( numCurrRHS );
00885   blockSize_ = numCurrRHS;
00886   for (int i=0; i<numCurrRHS; ++i) 
00887     { currIdx[i] = startPtr+i; }
00888 
00889   // Inform the linear problem of the current linear system to solve.
00890   problem_->setLSIndex( currIdx );
00891 
00893   // Parameter list
00894   Teuchos::ParameterList plist;
00895   plist.set("Num Blocks",numBlocks_);
00896   
00897   // Reset the status test.  
00898   outputTest_->reset();
00899   loaDetected_ = false;
00900 
00901   // Assume convergence is achieved, then let any failed convergence set this to false.
00902   bool isConverged = true;
00903 
00905   // BlockGmres solver
00906 
00907   Teuchos::RCP<PseudoBlockGmresIter<ScalarType,MV,OP> > block_gmres_iter
00908     = Teuchos::rcp( new PseudoBlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );  
00909 
00910   // Enter solve() iterations
00911   {
00912     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00913 
00914     while ( numRHS2Solve > 0 ) {
00915 
00916       // Reset the active / converged vectors from this block
00917       std::vector<int> convRHSIdx;
00918       std::vector<int> currRHSIdx( currIdx );
00919       currRHSIdx.resize(numCurrRHS);
00920 
00921       // Set the current number of blocks with the pseudo Gmres iteration.
00922       block_gmres_iter->setNumBlocks( numBlocks_ );
00923 
00924       // Reset the number of iterations.
00925       block_gmres_iter->resetNumIters();
00926 
00927       // Reset the number of calls that the status test output knows about.
00928       outputTest_->resetNumCalls();
00929 
00930       // Get a new state struct and initialize the solver.
00931       PseudoBlockGmresIterState<ScalarType,MV> newState;
00932 
00933       // Create the first block in the current Krylov basis for each right-hand side.
00934       std::vector<int> index(1);
00935       Teuchos::RCP<MV> tmpV, R_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
00936       newState.V.resize( blockSize_ );
00937       newState.Z.resize( blockSize_ );
00938       for (int i=0; i<blockSize_; ++i) {
00939         index[0]=i;
00940         tmpV = MVT::CloneViewNonConst( *R_0, index );
00941 
00942         // Get a matrix to hold the orthonormalization coefficients.
00943         Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
00944           = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
00945 
00946         // Orthonormalize the new V_0
00947         int rank = ortho_->normalize( *tmpV, tmpZ );
00948         TEST_FOR_EXCEPTION(rank != 1, PseudoBlockGmresSolMgrOrthoFailure,
00949             "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
00950 
00951         newState.V[i] = tmpV;
00952         newState.Z[i] = tmpZ;
00953       }
00954 
00955       newState.curDim = 0;
00956       block_gmres_iter->initialize(newState);
00957       int numRestarts = 0;
00958 
00959       while(1) {
00960         
00961         // tell block_gmres_iter to iterate
00962         try {
00963           block_gmres_iter->iterate();
00964           
00966           //
00967           // check convergence first
00968           //
00970           if ( convTest_->getStatus() == Passed ) {
00971         
00972             if ( expConvTest_->getLOADetected() ) {
00973               // we don't have convergence but we will deflate out the linear systems and move on.
00974               loaDetected_ = true;
00975               isConverged = false;
00976             }
00977         
00978             // Figure out which linear systems converged.
00979             std::vector<int> convIdx = expConvTest_->convIndices();
00980         
00981             // If the number of converged linear systems is equal to the
00982             // number of current linear systems, then we are done with this block.
00983             if (convIdx.size() == currRHSIdx.size())
00984               break;  // break from while(1){block_gmres_iter->iterate()}
00985         
00986             // Get a new state struct and initialize the solver.
00987             PseudoBlockGmresIterState<ScalarType,MV> defState;
00988         
00989             // Inform the linear problem that we are finished with this current linear system.
00990             problem_->setCurrLS();
00991         
00992             // Get the state.
00993             PseudoBlockGmresIterState<ScalarType,MV> oldState = block_gmres_iter->getState();
00994             int curDim = oldState.curDim;
00995         
00996             // Get a new state struct and reset currRHSIdx to have the right-hand sides that 
00997             // are left to converge for this block.
00998             int have = 0;
00999             std::vector<int> oldRHSIdx( currRHSIdx );
01000             std::vector<int> defRHSIdx;
01001             for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
01002               bool found = false;
01003               for (unsigned int j=0; j<convIdx.size(); ++j) {
01004                 if (currRHSIdx[i] == convIdx[j]) {
01005                   found = true;
01006                   break;
01007                 }
01008               }
01009               if (found) {
01010                 defRHSIdx.push_back( i );
01011               }
01012               else {
01013                 defState.V.push_back( Teuchos::rcp_const_cast<MV>( oldState.V[i] ) );
01014                 defState.Z.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.Z[i] ) );
01015                 defState.H.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseMatrix<int,ScalarType> >( oldState.H[i] ) );
01016                 defState.sn.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.sn[i] ) );
01017                 defState.cs.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,MagnitudeType> >(oldState.cs[i] ) );
01018                 currRHSIdx[have] = currRHSIdx[i];
01019                 have++;
01020               }
01021             }
01022             defRHSIdx.resize(currRHSIdx.size()-have);
01023             currRHSIdx.resize(have);
01024         
01025             // Compute the current solution that needs to be deflated if this solver has taken any steps.
01026             if (curDim) {
01027               Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01028               Teuchos::RCP<MV> defUpdate = MVT::CloneViewNonConst( *update, defRHSIdx );
01029               
01030               // Set the deflated indices so we can update the solution.
01031               problem_->setLSIndex( convIdx );
01032               
01033               // Update the linear problem.
01034               problem_->updateSolution( defUpdate, true );
01035             }
01036             
01037             // Set the remaining indices after deflation.
01038             problem_->setLSIndex( currRHSIdx );
01039             
01040             // Set the dimension of the subspace, which is the same as the old subspace size.
01041             defState.curDim = curDim;
01042             
01043             // Initialize the solver with the deflated system.
01044             block_gmres_iter->initialize(defState);
01045           }
01047           //
01048           // check for maximum iterations
01049           //
01051           else if ( maxIterTest_->getStatus() == Passed ) {
01052             // we don't have convergence
01053             isConverged = false;
01054             break;  // break from while(1){block_gmres_iter->iterate()}
01055           }
01057           //
01058           // check for restarting, i.e. the subspace is full
01059           //
01061           else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
01062         
01063             if ( numRestarts >= maxRestarts_ ) {
01064               isConverged = false;
01065               break; // break from while(1){block_gmres_iter->iterate()}
01066             }
01067             numRestarts++;
01068         
01069             printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
01070             
01071             // Update the linear problem.
01072             Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01073             problem_->updateSolution( update, true );
01074             
01075             // Get the state.
01076             PseudoBlockGmresIterState<ScalarType,MV> oldState = block_gmres_iter->getState();
01077             
01078             // Set the new state.
01079             PseudoBlockGmresIterState<ScalarType,MV> newstate;
01080             newstate.V.resize(currRHSIdx.size());
01081             newstate.Z.resize(currRHSIdx.size());
01082         
01083             // Compute the restart vectors
01084             // NOTE: Force the linear problem to update the current residual since the solution was updated.
01085             R_0 = MVT::Clone( *(problem_->getInitPrecResVec()), currRHSIdx.size() );
01086             problem_->computeCurrPrecResVec( &*R_0 );
01087             for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
01088               index[0] = i;  // index(1) vector declared on line 891
01089         
01090               tmpV = MVT::CloneViewNonConst( *R_0, index );
01091         
01092               // Get a matrix to hold the orthonormalization coefficients.
01093               Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
01094                 = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
01095               
01096               // Orthonormalize the new V_0
01097               int rank = ortho_->normalize( *tmpV, tmpZ );
01098               TEST_FOR_EXCEPTION(rank != 1 ,PseudoBlockGmresSolMgrOrthoFailure,
01099                   "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after the restart.");
01100               
01101               newstate.V[i] = tmpV;
01102               newstate.Z[i] = tmpZ;
01103             }
01104         
01105             // Initialize the solver.
01106             newstate.curDim = 0;
01107             block_gmres_iter->initialize(newstate);
01108         
01109           } // end of restarting
01110         
01112           //
01113           // we returned from iterate(), but none of our status tests Passed.
01114           // something is wrong, and it is probably our fault.
01115           //
01117         
01118           else {
01119             TEST_FOR_EXCEPTION(true,std::logic_error,
01120                 "Belos::PseudoBlockGmresSolMgr::solve(): Invalid return from PseudoBlockGmresIter::iterate().");
01121           }
01122         }
01123         catch (const PseudoBlockGmresIterOrthoFailure &e) {
01124         
01125           // Try to recover the most recent least-squares solution
01126           block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
01127         
01128           // Check to see if the most recent least-squares solution yielded convergence.
01129           sTest_->checkStatus( &*block_gmres_iter );
01130           if (convTest_->getStatus() != Passed)
01131             isConverged = false;
01132           break;
01133         }
01134         catch (const std::exception &e) {
01135           printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockGmresIter::iterate() at iteration " 
01136                                    << block_gmres_iter->getNumIters() << std::endl 
01137                                    << e.what() << std::endl;
01138           throw;
01139         }
01140       }
01141       
01142       // Compute the current solution.
01143       // Update the linear problem.
01144       if (nonnull(userConvStatusTest_)) {
01145         //std::cout << "\nnonnull(userConvStatusTest_)\n";
01146         Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01147         problem_->updateSolution( update, true );
01148       }
01149       else if (nonnull(expConvTest_->getSolution())) {
01150         //std::cout << "\nexpConvTest_->getSolution()\n";
01151         Teuchos::RCP<MV> newX = expConvTest_->getSolution();
01152         Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
01153         MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
01154       }
01155       else {
01156         //std::cout << "\nblock_gmres_iter->getCurrentUpdate()\n";
01157         Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
01158         problem_->updateSolution( update, true );
01159       }
01160 
01161       // Inform the linear problem that we are finished with this block linear system.
01162       problem_->setCurrLS();
01163       
01164       // Update indices for the linear systems to be solved.
01165       startPtr += numCurrRHS;
01166       numRHS2Solve -= numCurrRHS;
01167       if ( numRHS2Solve > 0 ) {
01168         numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01169 
01170         blockSize_ = numCurrRHS;
01171         currIdx.resize( numCurrRHS  );
01172         for (int i=0; i<numCurrRHS; ++i) 
01173         { currIdx[i] = startPtr+i; }
01174 
01175         // Adapt the status test quorum if we need to.
01176         if (defQuorum_ > blockSize_) {
01177           if (impConvTest_ != Teuchos::null)
01178             impConvTest_->setQuorum( blockSize_ );
01179           if (expConvTest_ != Teuchos::null)
01180             expConvTest_->setQuorum( blockSize_ );
01181         }
01182 
01183         // Set the next indices.
01184         problem_->setLSIndex( currIdx );
01185       }
01186       else {
01187         currIdx.resize( numRHS2Solve );
01188       }
01189 
01190     }// while ( numRHS2Solve > 0 )
01191 
01192   }
01193 
01194   // print final summary
01195   sTest_->print( printer_->stream(FinalSummary) );
01196  
01197   // print timing information
01198   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01199  
01200   // get iteration information for this solve
01201   numIters_ = maxIterTest_->getNumIters();
01202  
01203   if (!isConverged || loaDetected_) {
01204     return Unconverged; // return from PseudoBlockGmresSolMgr::solve() 
01205   }
01206   return Converged; // return from PseudoBlockGmresSolMgr::solve() 
01207 }
01208 
01209 //  This method requires the solver manager to return a std::string that describes itself.
01210 template<class ScalarType, class MV, class OP>
01211 std::string PseudoBlockGmresSolMgr<ScalarType,MV,OP>::description() const
01212 {
01213   std::ostringstream oss;
01214   oss << "Belos::PseudoBlockGmresSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01215   oss << "{";
01216   oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" <<blockSize_;
01217   oss << ", Num Blocks="<<numBlocks_<< ", Max Restarts=" << maxRestarts_;
01218   oss << "}";
01219   return oss.str();
01220 }
01221   
01222 } // end Belos namespace
01223 
01224 #endif /* BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines