BelosBlockCGSolMgr.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_BLOCK_CG_SOLMGR_HPP
00030 #define BELOS_BLOCK_CG_SOLMGR_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 
00039 #include "BelosLinearProblem.hpp"
00040 #include "BelosSolverManager.hpp"
00041 
00042 #include "BelosCGIter.hpp"
00043 #include "BelosBlockCGIter.hpp"
00044 #include "BelosDGKSOrthoManager.hpp"
00045 #include "BelosICGSOrthoManager.hpp"
00046 #include "BelosIMGSOrthoManager.hpp"
00047 #include "BelosStatusTestMaxIters.hpp"
00048 #include "BelosStatusTestGenResNorm.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 
00069 namespace Belos {
00070   
00072 
00073   
00080   class BlockCGSolMgrLinearProblemFailure : public BelosError {public:
00081     BlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00082     {}};
00083   
00090   class BlockCGSolMgrOrthoFailure : public BelosError {public:
00091     BlockCGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00092     {}};
00093   
00094   template<class ScalarType, class MV, class OP>
00095   class BlockCGSolMgr : public SolverManager<ScalarType,MV,OP> {
00096     
00097   private:
00098     typedef MultiVecTraits<ScalarType,MV> MVT;
00099     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00100     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00101     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00102     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00103     
00104   public:
00105     
00107 
00108 
00114      BlockCGSolMgr();
00115 
00143     BlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00144        const Teuchos::RCP<Teuchos::ParameterList> &pl );
00145     
00147     virtual ~BlockCGSolMgr() {};
00149     
00151 
00152     
00153     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00154       return *problem_;
00155     }
00156 
00159     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00160     
00163     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00164     
00170     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00171       return tuple(timerSolve_);
00172     }
00173 
00175     int getNumIters() const {
00176       return numIters_;
00177     }
00178 
00181     bool isLOADetected() const { return false; }
00182  
00184     
00186 
00187     
00188     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00189     
00190     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00191     
00193     
00195 
00196     
00214     ReturnType solve();
00215     
00217     
00220     
00222     std::string description() const;
00223     
00225     
00226   private:
00227 
00228     // Linear problem.
00229     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00230     
00231     // Output manager.
00232     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00233     Teuchos::RCP<std::ostream> outputStream_;
00234 
00235     // Status test.
00236     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00237     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00238     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00239     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00240 
00241     // Orthogonalization manager.
00242     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00243     
00244     // Current parameter list.
00245     Teuchos::RCP<ParameterList> params_;
00246     
00247     // Default solver values.
00248     static const MagnitudeType convtol_default_;
00249     static const MagnitudeType orthoKappa_default_;
00250     static const int maxIters_default_;
00251     static const bool adaptiveBlockSize_default_;
00252     static const bool showMaxResNormOnly_default_;
00253     static const int blockSize_default_;
00254     static const int verbosity_default_;
00255     static const int outputFreq_default_;
00256     static const std::string label_default_;
00257     static const std::string orthoType_default_;
00258     static const Teuchos::RCP<std::ostream> outputStream_default_;
00259 
00260     // Current solver values.
00261     MagnitudeType convtol_, orthoKappa_;
00262     int maxIters_, numIters_;
00263     int blockSize_, verbosity_, outputFreq_;
00264     bool adaptiveBlockSize_, showMaxResNormOnly_;
00265     std::string orthoType_; 
00266     
00267     // Timers.
00268     std::string label_;
00269     Teuchos::RCP<Teuchos::Time> timerSolve_;
00270 
00271     // Internal state variables.
00272     bool isSet_;
00273   };
00274 
00275 
00276 // Default solver values.
00277 template<class ScalarType, class MV, class OP>
00278 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType BlockCGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00279 
00280 template<class ScalarType, class MV, class OP>
00281 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType BlockCGSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00282 
00283 template<class ScalarType, class MV, class OP>
00284 const int BlockCGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00285 
00286 template<class ScalarType, class MV, class OP>
00287 const bool BlockCGSolMgr<ScalarType,MV,OP>::adaptiveBlockSize_default_ = true;
00288 
00289 template<class ScalarType, class MV, class OP>
00290 const bool BlockCGSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00291 
00292 template<class ScalarType, class MV, class OP>
00293 const int BlockCGSolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00294 
00295 template<class ScalarType, class MV, class OP>
00296 const int BlockCGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00297 
00298 template<class ScalarType, class MV, class OP>
00299 const int BlockCGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00300 
00301 template<class ScalarType, class MV, class OP>
00302 const std::string BlockCGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00303 
00304 template<class ScalarType, class MV, class OP>
00305 const std::string BlockCGSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00306 
00307 template<class ScalarType, class MV, class OP>
00308 const Teuchos::RCP<std::ostream> BlockCGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00309 
00310 
00311 // Empty Constructor
00312 template<class ScalarType, class MV, class OP>
00313 BlockCGSolMgr<ScalarType,MV,OP>::BlockCGSolMgr() :
00314   outputStream_(outputStream_default_),
00315   convtol_(convtol_default_),
00316   orthoKappa_(orthoKappa_default_),
00317   maxIters_(maxIters_default_),
00318   blockSize_(blockSize_default_),
00319   verbosity_(verbosity_default_),
00320   outputFreq_(outputFreq_default_),
00321   adaptiveBlockSize_(adaptiveBlockSize_default_),
00322   showMaxResNormOnly_(showMaxResNormOnly_default_),
00323   orthoType_(orthoType_default_),
00324   label_(label_default_),
00325   isSet_(false)
00326 {}
00327 
00328 
00329 // Basic Constructor
00330 template<class ScalarType, class MV, class OP>
00331 BlockCGSolMgr<ScalarType,MV,OP>::BlockCGSolMgr( 
00332                  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00333                  const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00334   problem_(problem),
00335   outputStream_(outputStream_default_),
00336   convtol_(convtol_default_),
00337   orthoKappa_(orthoKappa_default_),
00338   maxIters_(maxIters_default_),
00339   blockSize_(blockSize_default_),
00340   verbosity_(verbosity_default_),
00341   outputFreq_(outputFreq_default_),
00342   adaptiveBlockSize_(adaptiveBlockSize_default_),
00343   showMaxResNormOnly_(showMaxResNormOnly_default_),
00344   orthoType_(orthoType_default_),
00345   label_(label_default_),
00346   isSet_(false)
00347 {
00348   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00349 
00350   // If the parameter list pointer is null, then set the current parameters to the default parameter list.
00351   if ( !is_null(pl) ) {
00352     setParameters( pl );  
00353   }
00354 }
00355 
00356 template<class ScalarType, class MV, class OP>
00357 void BlockCGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00358 {
00359   // Create the internal parameter list if ones doesn't already exist.
00360   if (params_ == Teuchos::null) {
00361     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00362   }
00363   else {
00364     params->validateParameters(*getValidParameters());
00365   }
00366 
00367   // Check for maximum number of iterations
00368   if (params->isParameter("Maximum Iterations")) {
00369     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00370 
00371     // Update parameter in our list and in status test.
00372     params_->set("Maximum Iterations", maxIters_);
00373     if (maxIterTest_!=Teuchos::null)
00374       maxIterTest_->setMaxIters( maxIters_ );
00375   }
00376 
00377   // Check for blocksize
00378   if (params->isParameter("Block Size")) {
00379     blockSize_ = params->get("Block Size",blockSize_default_);    
00380     TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00381            "Belos::BlockCGSolMgr: \"Block Size\" must be strictly positive.");
00382 
00383     // Update parameter in our list.
00384     params_->set("Block Size", blockSize_);
00385   }
00386 
00387   // Check if the blocksize should be adaptive
00388   if (params->isParameter("Adapative Block Size")) {
00389     adaptiveBlockSize_ = params->get("Adaptive Block Size",adaptiveBlockSize_default_);
00390     
00391     // Update parameter in our list.
00392     params_->set("Adaptive Block Size", adaptiveBlockSize_);
00393   }
00394 
00395   // Check to see if the timer label changed.
00396   if (params->isParameter("Timer Label")) {
00397     std::string tempLabel = params->get("Timer Label", label_default_);
00398 
00399     // Update parameter in our list and solver timer
00400     if (tempLabel != label_) {
00401       label_ = tempLabel;
00402       params_->set("Timer Label", label_);
00403       std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
00404       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00405     }
00406   }
00407 
00408   // Check if the orthogonalization changed.
00409   if (params->isParameter("Orthogonalization")) {
00410     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00411     TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00412                         std::invalid_argument,
00413       "Belos::BlockCGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00414     if (tempOrthoType != orthoType_) {
00415       orthoType_ = tempOrthoType;
00416       // Create orthogonalization manager
00417       if (orthoType_=="DGKS") {
00418   if (orthoKappa_ <= 0) {
00419     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00420   }
00421   else {
00422     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00423     Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00424   }
00425       }
00426       else if (orthoType_=="ICGS") {
00427   ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00428       } 
00429       else if (orthoType_=="IMGS") {
00430   ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00431       } 
00432     }  
00433   }
00434 
00435   // Check which orthogonalization constant to use.
00436   if (params->isParameter("Orthogonalization Constant")) {
00437     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00438 
00439     // Update parameter in our list.
00440     params_->set("Orthogonalization Constant",orthoKappa_);
00441     if (orthoType_=="DGKS") {
00442       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00443   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00444       }
00445     } 
00446   }
00447 
00448   // Check for a change in verbosity level
00449   if (params->isParameter("Verbosity")) {
00450     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00451       verbosity_ = params->get("Verbosity", verbosity_default_);
00452     } else {
00453       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00454     }
00455 
00456     // Update parameter in our list.
00457     params_->set("Verbosity", verbosity_);
00458     if (printer_ != Teuchos::null)
00459       printer_->setVerbosity(verbosity_);
00460   }
00461 
00462   // output stream
00463   if (params->isParameter("Output Stream")) {
00464     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00465 
00466     // Update parameter in our list.
00467     params_->set("Output Stream", outputStream_);
00468     if (printer_ != Teuchos::null)
00469       printer_->setOStream( outputStream_ );
00470   }
00471 
00472   // frequency level
00473   if (verbosity_ & Belos::StatusTestDetails) {
00474     if (params->isParameter("Output Frequency")) {
00475       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00476     }
00477 
00478     // Update parameter in out list and output status test.
00479     params_->set("Output Frequency", outputFreq_);
00480     if (outputTest_ != Teuchos::null)
00481       outputTest_->setOutputFrequency( outputFreq_ );
00482   }
00483 
00484   // Create output manager if we need to.
00485   if (printer_ == Teuchos::null) {
00486     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00487   }  
00488   
00489   // Convergence
00490   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00491   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00492 
00493   // Check for convergence tolerance
00494   if (params->isParameter("Convergence Tolerance")) {
00495     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00496 
00497     // Update parameter in our list and residual tests.
00498     params_->set("Convergence Tolerance", convtol_);
00499     if (convTest_ != Teuchos::null)
00500       convTest_->setTolerance( convtol_ );
00501   }
00502   
00503   if (params->isParameter("Show Maximum Residual Norm Only")) {
00504     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00505 
00506     // Update parameter in our list and residual tests
00507     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00508     if (convTest_ != Teuchos::null)
00509       convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00510   }
00511 
00512   // Create status tests if we need to.
00513 
00514   // Basic test checks maximum iterations and native residual.
00515   if (maxIterTest_ == Teuchos::null)
00516     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00517   
00518   // Implicit residual test, using the native residual to determine if convergence was achieved.
00519   if (convTest_ == Teuchos::null)
00520     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
00521   
00522   if (sTest_ == Teuchos::null)
00523     sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00524   
00525   if (outputTest_ == Teuchos::null) {
00526     if (outputFreq_ > 0) {
00527       outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_, 
00528                     sTest_, 
00529                     outputFreq_, 
00530                     Passed+Failed+Undefined ) ); 
00531     }
00532     else {
00533       outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_, 
00534                     sTest_, 1 ) );
00535     }
00536   }
00537 
00538   // Create orthogonalization manager if we need to.
00539   if (ortho_ == Teuchos::null) {
00540     if (orthoType_=="DGKS") {
00541       if (orthoKappa_ <= 0) {
00542   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00543       }
00544       else {
00545   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00546   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00547       }
00548     }
00549     else if (orthoType_=="ICGS") {
00550       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00551     } 
00552     else if (orthoType_=="IMGS") {
00553       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00554     } 
00555     else {
00556       TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00557        "Belos::BlockCGSolMgr(): Invalid orthogonalization type.");
00558     }  
00559   }
00560 
00561   // Create the timer if we need to.
00562   if (timerSolve_ == Teuchos::null) {
00563     std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
00564     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00565   }
00566 
00567   // Inform the solver manager that the current parameters were set.
00568   isSet_ = true;
00569 }
00570 
00571     
00572 template<class ScalarType, class MV, class OP>
00573 Teuchos::RCP<const Teuchos::ParameterList> 
00574 BlockCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00575 {
00576   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00577   
00578   // Set all the valid parameters and their default values.
00579   if(is_null(validPL)) {
00580     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00581     pl->set("Convergence Tolerance", convtol_default_,
00582       "The relative residual tolerance that needs to be achieved by the\n"
00583       "iterative solver in order for the linear system to be declared converged.");
00584     pl->set("Maximum Iterations", maxIters_default_,
00585       "The maximum number of block iterations allowed for each\n"
00586       "set of RHS solved.");
00587     pl->set("Block Size", blockSize_default_,
00588       "The number of vectors in each block.");
00589     pl->set("Adaptive Block Size", adaptiveBlockSize_default_,
00590       "Whether the solver manager should adapt to the block size\n"
00591       "based on the number of RHS to solve.");
00592     pl->set("Verbosity", verbosity_default_,
00593       "What type(s) of solver information should be outputted\n"
00594       "to the output stream.");
00595     pl->set("Output Frequency", outputFreq_default_,
00596       "How often convergence information should be outputted\n"
00597       "to the output stream.");  
00598     pl->set("Output Stream", outputStream_default_,
00599       "A reference-counted pointer to the output stream where all\n"
00600       "solver output is sent.");
00601     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00602       "When convergence information is printed, only show the maximum\n"
00603       "relative residual norm when the block size is greater than one.");
00604     pl->set("Timer Label", label_default_,
00605       "The string to use as a prefix for the timer labels.");
00606     //  pl->set("Restart Timers", restartTimers_);
00607     pl->set("Orthogonalization", orthoType_default_,
00608       "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
00609     pl->set("Orthogonalization Constant",orthoKappa_default_,
00610       "The constant used by DGKS orthogonalization to determine\n"
00611       "whether another step of classical Gram-Schmidt is necessary.");
00612     validPL = pl;
00613   }
00614   return validPL;
00615 }
00616 
00617   
00618 // solve()
00619 template<class ScalarType, class MV, class OP>
00620 ReturnType BlockCGSolMgr<ScalarType,MV,OP>::solve() {
00621 
00622   // Set the current parameters if they were not set before.
00623   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00624   // then didn't set any parameters using setParameters().
00625   if (!isSet_) {
00626     setParameters(Teuchos::parameterList(*getValidParameters()));
00627   }
00628 
00629   Teuchos::BLAS<int,ScalarType> blas;
00630   Teuchos::LAPACK<int,ScalarType> lapack;
00631   
00632   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),BlockCGSolMgrLinearProblemFailure,
00633                      "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00634 
00635   // Create indices for the linear systems to be solved.
00636   int startPtr = 0;
00637   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00638   int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
00639 
00640   std::vector<int> currIdx, currIdx2;
00641   //  If an adaptive block size is allowed then only the linear systems that need to be solved are solved.
00642   //  Otherwise, the index set is generated that informs the linear problem that some linear systems are augmented.
00643   if ( adaptiveBlockSize_ ) {
00644     blockSize_ = numCurrRHS;
00645     currIdx.resize( numCurrRHS  );
00646     currIdx2.resize( numCurrRHS  );
00647     for (int i=0; i<numCurrRHS; ++i) 
00648       { currIdx[i] = startPtr+i; currIdx2[i]=i; }
00649     
00650   }
00651   else {
00652     currIdx.resize( blockSize_ );
00653     currIdx2.resize( blockSize_ );
00654     for (int i=0; i<numCurrRHS; ++i) 
00655       { currIdx[i] = startPtr+i; currIdx2[i]=i; }
00656     for (int i=numCurrRHS; i<blockSize_; ++i) 
00657       { currIdx[i] = -1; currIdx2[i] = i; }
00658   }
00659 
00660   // Inform the linear problem of the current linear system to solve.
00661   problem_->setLSIndex( currIdx );
00662 
00664   // Parameter list
00665   Teuchos::ParameterList plist;
00666   plist.set("Block Size",blockSize_);
00667   
00668   // Reset the status test.  
00669   outputTest_->reset();
00670 
00671   // Assume convergence is achieved, then let any failed convergence set this to false.
00672   bool isConverged = true;  
00673 
00675   // BlockCG solver
00676 
00677   Teuchos::RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
00678   if (blockSize_ == 1)
00679     block_cg_iter = Teuchos::rcp( new CGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
00680   else
00681     block_cg_iter = Teuchos::rcp( new BlockCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
00682   
00683 
00684   // Enter solve() iterations
00685   {
00686     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00687 
00688     while ( numRHS2Solve > 0 ) {
00689       //
00690       // Reset the active / converged vectors from this block
00691       std::vector<int> convRHSIdx;
00692       std::vector<int> currRHSIdx( currIdx );
00693       currRHSIdx.resize(numCurrRHS);
00694 
00695       // Reset the number of iterations.
00696       block_cg_iter->resetNumIters();
00697 
00698       // Reset the number of calls that the status test output knows about.
00699       outputTest_->resetNumCalls();
00700 
00701       // Get the current residual for this block of linear systems.
00702       Teuchos::RCP<MV> R_0 = MVT::CloneView( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
00703 
00704       // Set the new state and initialize the solver.
00705       CGIterationState<ScalarType,MV> newstate;
00706       newstate.R = R_0;
00707       block_cg_iter->initializeCG(newstate);
00708 
00709       while(1) {
00710   
00711   // tell block_cg_iter to iterate
00712   try {
00713     block_cg_iter->iterate();
00714     
00716     //
00717     // check convergence first
00718     //
00720     if ( convTest_->getStatus() == Passed ) {
00721       // We have convergence of at least one linear system.
00722 
00723       // Figure out which linear systems converged.
00724       std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
00725 
00726       // If the number of converged linear systems is equal to the
00727             // number of current linear systems, then we are done with this block.
00728       if (convIdx.size() == currRHSIdx.size())
00729         break;  // break from while(1){block_cg_iter->iterate()}
00730 
00731       // Inform the linear problem that we are finished with this current linear system.
00732       problem_->setCurrLS();
00733 
00734       // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
00735       int have = 0;
00736             std::vector<int> unconvIdx(currRHSIdx.size());
00737       for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
00738         bool found = false;
00739         for (unsigned int j=0; j<convIdx.size(); ++j) {
00740     if (currRHSIdx[i] == convIdx[j]) {
00741       found = true;
00742       break;
00743     }
00744         }
00745         if (!found) {
00746                 currIdx2[have] = currIdx2[i];
00747     currRHSIdx[have++] = currRHSIdx[i];
00748         }
00749               else {
00750               } 
00751       }
00752       currRHSIdx.resize(have);
00753       currIdx2.resize(have);
00754 
00755       // Set the remaining indices after deflation.
00756       problem_->setLSIndex( currRHSIdx );
00757 
00758       // Get the current residual std::vector.
00759       std::vector<MagnitudeType> norms;
00760             R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
00761       for (int i=0; i<have; ++i) { currIdx2[i] = i; }
00762 
00763       // Set the new blocksize for the solver.
00764       block_cg_iter->setBlockSize( have );
00765 
00766       // Set the new state and initialize the solver.
00767       CGIterationState<ScalarType,MV> defstate;
00768       defstate.R = R_0;
00769       block_cg_iter->initializeCG(defstate);
00770     }
00772     //
00773     // check for maximum iterations
00774     //
00776     else if ( maxIterTest_->getStatus() == Passed ) {
00777       // we don't have convergence
00778       isConverged = false;
00779       break;  // break from while(1){block_cg_iter->iterate()}
00780     }
00781 
00783     //
00784     // we returned from iterate(), but none of our status tests Passed.
00785     // something is wrong, and it is probably our fault.
00786     //
00788 
00789     else {
00790       TEST_FOR_EXCEPTION(true,std::logic_error,
00791              "Belos::BlockCGSolMgr::solve(): Invalid return from CGIteration::iterate().");
00792     }
00793   }
00794   catch (std::exception e) {
00795     printer_->stream(Errors) << "Error! Caught std::exception in CGIteration::iterate() at iteration " 
00796            << block_cg_iter->getNumIters() << std::endl 
00797            << e.what() << std::endl;
00798     throw;
00799   }
00800       }
00801       
00802       // Inform the linear problem that we are finished with this block linear system.
00803       problem_->setCurrLS();
00804       
00805       // Update indices for the linear systems to be solved.
00806       startPtr += numCurrRHS;
00807       numRHS2Solve -= numCurrRHS;
00808       if ( numRHS2Solve > 0 ) {
00809   numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
00810 
00811   if ( adaptiveBlockSize_ ) {
00812           blockSize_ = numCurrRHS;
00813     currIdx.resize( numCurrRHS  );
00814     currIdx2.resize( numCurrRHS  );
00815     for (int i=0; i<numCurrRHS; ++i) 
00816       { currIdx[i] = startPtr+i; currIdx2[i] = i; }   
00817   }
00818   else {
00819     currIdx.resize( blockSize_ );
00820     currIdx2.resize( blockSize_ );
00821     for (int i=0; i<numCurrRHS; ++i) 
00822       { currIdx[i] = startPtr+i; currIdx2[i] = i; }
00823     for (int i=numCurrRHS; i<blockSize_; ++i) 
00824       { currIdx[i] = -1; currIdx2[i] = i; }
00825   }
00826   // Set the next indices.
00827   problem_->setLSIndex( currIdx );
00828 
00829   // Set the new blocksize for the solver.
00830   block_cg_iter->setBlockSize( blockSize_ );  
00831       }
00832       else {
00833         currIdx.resize( numRHS2Solve );
00834       }
00835       
00836     }// while ( numRHS2Solve > 0 )
00837     
00838   }
00839 
00840   // print final summary
00841   sTest_->print( printer_->stream(FinalSummary) );
00842  
00843   // print timing information
00844   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
00845  
00846   // get iteration information for this solve
00847   numIters_ = maxIterTest_->getNumIters();
00848  
00849   if (!isConverged) {
00850     return Unconverged; // return from BlockCGSolMgr::solve() 
00851   }
00852   return Converged; // return from BlockCGSolMgr::solve() 
00853 }
00854 
00855 //  This method requires the solver manager to return a std::string that describes itself.
00856 template<class ScalarType, class MV, class OP>
00857 std::string BlockCGSolMgr<ScalarType,MV,OP>::description() const
00858 {
00859   std::ostringstream oss;
00860   oss << "Belos::BlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00861   oss << "{";
00862   oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" << blockSize_;
00863   oss << "}";
00864   return oss.str();
00865 }
00866   
00867 } // end Belos namespace
00868 
00869 #endif /* BELOS_BLOCK_CG_SOLMGR_HPP */

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