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

Generated on Tue Jul 13 09:27:02 2010 for Belos by  doxygen 1.4.7