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

Generated on Wed May 12 21:45:51 2010 for Belos by  doxygen 1.4.7