BelosPseudoBlockCGSolMgr.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_CG_SOLMGR_HPP
00030 #define BELOS_PSEUDO_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 "BelosPseudoBlockCGIter.hpp"
00043 #include "BelosStatusTestMaxIters.hpp"
00044 #include "BelosStatusTestGenResNorm.hpp"
00045 #include "BelosStatusTestCombo.hpp"
00046 #include "BelosStatusTestOutput.hpp"
00047 #include "BelosOutputManager.hpp"
00048 #include "Teuchos_BLAS.hpp"
00049 #include "Teuchos_TimeMonitor.hpp"
00050 
00066 namespace Belos {
00067   
00069 
00070   
00077   class PseudoBlockCGSolMgrLinearProblemFailure : public BelosError {public:
00078     PseudoBlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00079     {}};
00080   
00087   class PseudoBlockCGSolMgrOrthoFailure : public BelosError {public:
00088     PseudoBlockCGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00089     {}};
00090   
00091   template<class ScalarType, class MV, class OP>
00092   class PseudoBlockCGSolMgr : public SolverManager<ScalarType,MV,OP> {
00093     
00094   private:
00095     typedef MultiVecTraits<ScalarType,MV> MVT;
00096     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00097     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00098     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00099     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00100     
00101   public:
00102     
00104 
00105 
00111     PseudoBlockCGSolMgr();
00112     
00121     PseudoBlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00122                    const Teuchos::RCP<Teuchos::ParameterList> &pl );
00123     
00125     virtual ~PseudoBlockCGSolMgr() {};
00127     
00129 
00130     
00131     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00132       return *problem_;
00133     }
00134     
00137     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00138    
00141     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00142  
00148     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00149       return tuple(timerSolve_);
00150     }
00151 
00153     int getNumIters() const {
00154       return numIters_;
00155     }
00156   
00160     bool isLOADetected() const { return false; }
00161   
00163     
00165 
00166    
00168     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00169    
00171     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00172     
00174 
00176 
00177 
00181     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00183 
00185 
00186     
00204     ReturnType solve();
00205     
00207     
00210     
00212     std::string description() const;
00213     
00215     
00216   private:
00217 
00218     // Method to convert std::string to enumerated type for residual.
00219     Belos::ScaleType convertStringToScaleType( const std::string& scaleType ) {
00220       if (scaleType == "Norm of Initial Residual") {
00221         return Belos::NormOfInitRes;
00222       } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00223         return Belos::NormOfPrecInitRes;
00224       } else if (scaleType == "Norm of RHS") {
00225         return Belos::NormOfRHS;
00226       } else if (scaleType == "None") {
00227         return Belos::None;
00228       } else 
00229         TEST_FOR_EXCEPTION( true ,std::logic_error,
00230           "Belos::PseudoBlockCGSolMgr(): Invalid residual scaling type.");
00231     }
00232 
00233     // Linear problem.
00234     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00235     
00236     // Output manager.
00237     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00238     Teuchos::RCP<std::ostream> outputStream_;
00239 
00240     // Status test.
00241     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00242     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00243     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00244     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00245 
00246     // Orthogonalization manager.
00247     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00248 
00249      // Current parameter list.
00250     Teuchos::RCP<ParameterList> params_;
00251    
00252     // Default solver values.
00253     static const MagnitudeType convtol_default_;
00254     static const int maxIters_default_;
00255     static const bool showMaxResNormOnly_default_;
00256     static const int verbosity_default_;
00257     static const int outputFreq_default_;
00258     static const int defQuorum_default_;
00259     static const std::string resScale_default_; 
00260     static const std::string label_default_;
00261     static const Teuchos::RCP<std::ostream> outputStream_default_;
00262 
00263     // Current solver values.
00264     MagnitudeType convtol_;
00265     int maxIters_, numIters_;
00266     int verbosity_, outputFreq_, defQuorum_;
00267     bool showMaxResNormOnly_;
00268     std::string resScale_;       
00269  
00270     // Timers.
00271     std::string label_;
00272     Teuchos::RCP<Teuchos::Time> timerSolve_;
00273 
00274     // Internal state variables.
00275     bool isSet_, isSTSet_;
00276   };
00277   
00278   
00279 // Default solver values.
00280 template<class ScalarType, class MV, class OP>
00281 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType PseudoBlockCGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00282 
00283 template<class ScalarType, class MV, class OP>
00284 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00285 
00286 template<class ScalarType, class MV, class OP>
00287 const bool PseudoBlockCGSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00288 
00289 template<class ScalarType, class MV, class OP>
00290 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00291 
00292 template<class ScalarType, class MV, class OP>
00293 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00294 
00295 template<class ScalarType, class MV, class OP>
00296 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::defQuorum_default_ = 1;
00297 
00298 template<class ScalarType, class MV, class OP>
00299 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::resScale_default_ = "Norm of Initial Residual";
00300 
00301 template<class ScalarType, class MV, class OP>
00302 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00303 
00304 template<class ScalarType, class MV, class OP>
00305 const Teuchos::RCP<std::ostream> PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00306 
00307 
00308 // Empty Constructor
00309 template<class ScalarType, class MV, class OP>
00310 PseudoBlockCGSolMgr<ScalarType,MV,OP>::PseudoBlockCGSolMgr() :
00311   outputStream_(outputStream_default_),
00312   convtol_(convtol_default_),
00313   maxIters_(maxIters_default_),
00314   verbosity_(verbosity_default_),
00315   outputFreq_(outputFreq_default_),
00316   defQuorum_(defQuorum_default_),
00317   showMaxResNormOnly_(showMaxResNormOnly_default_),
00318   resScale_(resScale_default_),
00319   label_(label_default_),
00320   isSet_(false),
00321   isSTSet_(false)
00322 {}
00323 
00324 // Basic Constructor
00325 template<class ScalarType, class MV, class OP>
00326 PseudoBlockCGSolMgr<ScalarType,MV,OP>::PseudoBlockCGSolMgr( 
00327                const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00328                const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00329   problem_(problem),
00330   outputStream_(outputStream_default_),
00331   convtol_(convtol_default_),
00332   maxIters_(maxIters_default_),
00333   verbosity_(verbosity_default_),
00334   outputFreq_(outputFreq_default_),
00335   defQuorum_(defQuorum_default_),
00336   showMaxResNormOnly_(showMaxResNormOnly_default_),
00337   resScale_(resScale_default_),
00338   label_(label_default_),
00339   isSet_(false),
00340   isSTSet_(false)
00341 {
00342   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00343 
00344   // If the parameter list pointer is null, then set the current parameters to the default parameter list.
00345   if (!is_null(pl)) {
00346     // Set the parameters using the list that was passed in.
00347     setParameters( pl );  
00348   }
00349 }
00350 
00351 template<class ScalarType, class MV, class OP>
00352 void PseudoBlockCGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00353 {
00354   // Create the internal parameter list if ones doesn't already exist.
00355   if (params_ == Teuchos::null) {
00356     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00357   }
00358   else {
00359     params->validateParameters(*getValidParameters());
00360   }
00361 
00362   // Check for maximum number of iterations
00363   if (params->isParameter("Maximum Iterations")) {
00364     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00365 
00366     // Update parameter in our list and in status test.
00367     params_->set("Maximum Iterations", maxIters_);
00368     if (maxIterTest_!=Teuchos::null)
00369       maxIterTest_->setMaxIters( maxIters_ );
00370   }
00371 
00372   // Check to see if the timer label changed.
00373   if (params->isParameter("Timer Label")) {
00374     std::string tempLabel = params->get("Timer Label", label_default_);
00375 
00376     // Update parameter in our list and solver timer
00377     if (tempLabel != label_) {
00378       label_ = tempLabel;
00379       params_->set("Timer Label", label_);
00380       std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00381       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00382     }
00383   }
00384 
00385   // Check for a change in verbosity level
00386   if (params->isParameter("Verbosity")) {
00387     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00388       verbosity_ = params->get("Verbosity", verbosity_default_);
00389     } else {
00390       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00391     }
00392 
00393     // Update parameter in our list.
00394     params_->set("Verbosity", verbosity_);
00395     if (printer_ != Teuchos::null)
00396       printer_->setVerbosity(verbosity_);
00397   }
00398 
00399   // output stream
00400   if (params->isParameter("Output Stream")) {
00401     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00402 
00403     // Update parameter in our list.
00404     params_->set("Output Stream", outputStream_);
00405     if (printer_ != Teuchos::null)
00406       printer_->setOStream( outputStream_ );
00407   }
00408 
00409   // frequency level
00410   if (verbosity_ & Belos::StatusTestDetails) {
00411     if (params->isParameter("Output Frequency")) {
00412       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00413     }
00414 
00415     // Update parameter in out list and output status test.
00416     params_->set("Output Frequency", outputFreq_);
00417     if (outputTest_ != Teuchos::null)
00418       outputTest_->setOutputFrequency( outputFreq_ );
00419   }
00420 
00421   // Create output manager if we need to.
00422   if (printer_ == Teuchos::null) {
00423     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00424   }
00425   
00426   // Convergence
00427   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00428   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00429 
00430   // Check for convergence tolerance
00431   if (params->isParameter("Convergence Tolerance")) {
00432     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00433 
00434     // Update parameter in our list and residual tests.
00435     params_->set("Convergence Tolerance", convtol_);
00436     if (convTest_ != Teuchos::null)
00437       convTest_->setTolerance( convtol_ );
00438   }
00439 
00440   if (params->isParameter("Show Maximum Residual Norm Only")) {
00441     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00442 
00443     // Update parameter in our list and residual tests
00444     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00445     if (convTest_ != Teuchos::null)
00446       convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00447   }
00448 
00449   // Check for a change in scaling, if so we need to build new residual tests.
00450   bool newResTest = false;
00451   if (params->isParameter("Residual Scaling")) {
00452     std::string tempResScale = Teuchos::getParameter<std::string>( *params, "Residual Scaling" );
00453 
00454     // Only update the scaling if it's different.
00455     if (resScale_ != tempResScale) {
00456       Belos::ScaleType resScaleType = convertStringToScaleType( tempResScale );
00457       resScale_ = tempResScale;
00458 
00459       // Update parameter in our list and residual tests
00460       params_->set("Residual Scaling", resScale_);
00461       if (convTest_ != Teuchos::null) {
00462         try {
00463           convTest_->defineScaleForm( resScaleType, Belos::TwoNorm );
00464         }
00465         catch (std::exception& e) {
00466     // Make sure the convergence test gets constructed again.
00467     newResTest = true;
00468         }
00469       }
00470     }      
00471   }
00472 
00473   // Get the deflation quorum, or number of converged systems before deflation is allowed
00474   if (params->isParameter("Deflation Quorum")) {
00475     defQuorum_ = params->get("Deflation Quorum", defQuorum_);
00476     params_->set("Deflation Quorum", defQuorum_);
00477     if (convTest_ != Teuchos::null)
00478       convTest_->setQuorum( defQuorum_ );
00479   }
00480 
00481   // Create status tests if we need to.
00482 
00483   // Basic test checks maximum iterations and native residual.
00484   if (maxIterTest_ == Teuchos::null)
00485     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00486   
00487   // Implicit residual test, using the native residual to determine if convergence was achieved.
00488   if (convTest_ == Teuchos::null || newResTest) {
00489     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, defQuorum_, showMaxResNormOnly_ ) );
00490     convTest_->defineScaleForm( convertStringToScaleType( resScale_ ), Belos::TwoNorm );
00491   } 
00492 
00493   if (sTest_ == Teuchos::null || newResTest)
00494     sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00495   
00496   if (outputTest_ == Teuchos::null || newResTest) {
00497     if (outputFreq_ > 0) {
00498       outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_, 
00499                     sTest_, 
00500                     outputFreq_, 
00501                     Passed+Failed+Undefined ) ); 
00502     }
00503     else {
00504       outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_, 
00505                     sTest_, 1 ) );
00506     }
00507   }
00508 
00509   // Create the timer if we need to.
00510   if (timerSolve_ == Teuchos::null) {
00511     std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00512     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00513   }
00514 
00515   // Inform the solver manager that the current parameters were set.
00516   isSet_ = true;
00517 }
00518 
00519     
00520 template<class ScalarType, class MV, class OP>
00521 Teuchos::RCP<const Teuchos::ParameterList>
00522 PseudoBlockCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00523 {
00524   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00525   if (is_null(validPL)) {
00526     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00527     // Set all the valid parameters and their default values.
00528     pl= Teuchos::rcp( new Teuchos::ParameterList() );
00529     pl->set("Convergence Tolerance", convtol_default_,
00530       "The relative residual tolerance that needs to be achieved by the\n"
00531       "iterative solver in order for the linera system to be declared converged.");
00532     pl->set("Maximum Iterations", maxIters_default_,
00533       "The maximum number of block iterations allowed for each\n"
00534       "set of RHS solved.");
00535     pl->set("Verbosity", verbosity_default_,
00536       "What type(s) of solver information should be outputted\n"
00537       "to the output stream.");
00538     pl->set("Output Frequency", outputFreq_default_,
00539       "How often convergence information should be outputted\n"
00540       "to the output stream.");  
00541     pl->set("Deflation Quorum", defQuorum_default_,
00542       "The number of linear systems that need to converge before\n"
00543       "they are deflated.  This number should be <= block size.");
00544     pl->set("Output Stream", outputStream_default_,
00545       "A reference-counted pointer to the output stream where all\n"
00546       "solver output is sent.");
00547     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00548       "When convergence information is printed, only show the maximum\n"
00549       "relative residual norm when the block size is greater than one.");
00550     pl->set("Residual Scaling", resScale_default_,
00551       "The type of scaling used in the residual convergence test.");
00552     pl->set("Timer Label", label_default_,
00553       "The string to use as a prefix for the timer labels.");
00554     //  defaultParams_->set("Restart Timers", restartTimers_);
00555     validPL = pl;
00556   }
00557   return validPL;
00558 }
00559 
00560 
00561 // solve()
00562 template<class ScalarType, class MV, class OP>
00563 ReturnType PseudoBlockCGSolMgr<ScalarType,MV,OP>::solve() {
00564 
00565   // Set the current parameters if they were not set before.
00566   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00567   // then didn't set any parameters using setParameters().
00568   if (!isSet_) { setParameters( params_ ); }
00569   
00570   Teuchos::BLAS<int,ScalarType> blas;
00571   
00572   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockCGSolMgrLinearProblemFailure,
00573                      "Belos::PseudoBlockCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00574 
00575   // Create indices for the linear systems to be solved.
00576   int startPtr = 0;
00577   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00578   int numCurrRHS = numRHS2Solve;
00579 
00580   std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
00581   for (int i=0; i<numRHS2Solve; ++i) {
00582     currIdx[i] = startPtr+i; 
00583     currIdx2[i]=i;
00584   }
00585 
00586   // Inform the linear problem of the current linear system to solve.
00587   problem_->setLSIndex( currIdx );
00588 
00590   // Parameter list
00591   Teuchos::ParameterList plist;
00592 
00593   // Reset the status test.  
00594   outputTest_->reset();
00595 
00596   // Assume convergence is achieved, then let any failed convergence set this to false.
00597   bool isConverged = true;  
00598 
00600   // Pseudo-Block CG solver
00601 
00602   Teuchos::RCP<PseudoBlockCGIter<ScalarType,MV,OP> > block_cg_iter
00603     = Teuchos::rcp( new PseudoBlockCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );  
00604 
00605   // Enter solve() iterations
00606   {
00607     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00608 
00609     while ( numRHS2Solve > 0 ) {
00610 
00611       // Reset the active / converged vectors from this block
00612       std::vector<int> convRHSIdx;
00613       std::vector<int> currRHSIdx( currIdx );
00614       currRHSIdx.resize(numCurrRHS);
00615 
00616       // Reset the number of iterations.
00617       block_cg_iter->resetNumIters();
00618 
00619       // Reset the number of calls that the status test output knows about.
00620       outputTest_->resetNumCalls();
00621 
00622       // Get the current residual for this block of linear systems.
00623       Teuchos::RCP<MV> R_0 = MVT::CloneView( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
00624 
00625       // Get a new state struct and initialize the solver.
00626       CGIterationState<ScalarType,MV> newState;
00627       newState.R = R_0;
00628       block_cg_iter->initializeCG(newState);
00629 
00630       while(1) {
00631   
00632   // tell block_gmres_iter to iterate
00633   try {
00634     block_cg_iter->iterate();
00635     
00637     //
00638     // check convergence first
00639     //
00641     if ( convTest_->getStatus() == Passed ) {
00642       
00643       // Figure out which linear systems converged.
00644       std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
00645 
00646       // If the number of converged linear systems is equal to the
00647             // number of current linear systems, then we are done with this block.
00648       if (convIdx.size() == currRHSIdx.size())
00649         break;  // break from while(1){block_cg_iter->iterate()}
00650 
00651       // Inform the linear problem that we are finished with this current linear system.
00652       problem_->setCurrLS();
00653 
00654       // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
00655       int have = 0;
00656             std::vector<int> unconvIdx(currRHSIdx.size());
00657       for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
00658         bool found = false;
00659         for (unsigned int j=0; j<convIdx.size(); ++j) {
00660     if (currRHSIdx[i] == convIdx[j]) {
00661       found = true;
00662       break;
00663     }
00664         }
00665         if (!found) {
00666                 currIdx2[have] = currIdx2[i];
00667     currRHSIdx[have++] = currRHSIdx[i];
00668         }
00669       }
00670       currRHSIdx.resize(have);
00671       currIdx2.resize(have);
00672       
00673       // Set the remaining indices after deflation.
00674       problem_->setLSIndex( currRHSIdx );
00675       
00676       // Get the current residual vector.
00677       std::vector<MagnitudeType> norms;
00678             R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
00679       for (int i=0; i<have; ++i) { currIdx2[i] = i; }
00680 
00681       // Set the new state and initialize the solver.
00682       CGIterationState<ScalarType,MV> defstate;
00683       defstate.R = R_0;
00684       block_cg_iter->initializeCG(defstate);
00685     }
00686 
00688     //
00689     // check for maximum iterations
00690     //
00692     else if ( maxIterTest_->getStatus() == Passed ) {
00693       // we don't have convergence
00694       isConverged = false;
00695       break;  // break from while(1){block_cg_iter->iterate()}
00696     }
00697 
00699     //
00700     // we returned from iterate(), but none of our status tests Passed.
00701     // something is wrong, and it is probably our fault.
00702     //
00704 
00705     else {
00706       TEST_FOR_EXCEPTION(true,std::logic_error,
00707              "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
00708     }
00709   } 
00710   catch (const std::exception &e) {
00711     printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration " 
00712            << block_cg_iter->getNumIters() << std::endl 
00713            << e.what() << std::endl;
00714     throw;
00715   }
00716       }
00717       
00718       // Inform the linear problem that we are finished with this block linear system.
00719       problem_->setCurrLS();
00720       
00721       // Update indices for the linear systems to be solved.
00722       startPtr += numCurrRHS;
00723       numRHS2Solve -= numCurrRHS;
00724 
00725       if ( numRHS2Solve > 0 ) { 
00726 
00727   numCurrRHS = numRHS2Solve;  
00728   currIdx.resize( numCurrRHS );
00729   currIdx2.resize( numCurrRHS );
00730   for (int i=0; i<numCurrRHS; ++i) 
00731     { currIdx[i] = startPtr+i; currIdx2[i] = i; }
00732   
00733   // Set the next indices.
00734   problem_->setLSIndex( currIdx );  
00735       }
00736       else {
00737   currIdx.resize( numRHS2Solve );
00738       }
00739       
00740     }// while ( numRHS2Solve > 0 )
00741     
00742   }
00743                      
00744   // print final summary
00745   sTest_->print( printer_->stream(FinalSummary) );
00746  
00747   // print timing information
00748   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
00749  
00750   // get iteration information for this solve
00751   numIters_ = maxIterTest_->getNumIters();
00752  
00753   if (!isConverged ) {
00754     return Unconverged; // return from PseudoBlockCGSolMgr::solve() 
00755   }
00756   return Converged; // return from PseudoBlockCGSolMgr::solve() 
00757 }
00758 
00759 //  This method requires the solver manager to return a std::string that describes itself.
00760 template<class ScalarType, class MV, class OP>
00761 std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::description() const
00762 {
00763   std::ostringstream oss;
00764   oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00765   oss << "{";
00766   oss << "}";
00767   return oss.str();
00768 }
00769   
00770 } // end Belos namespace
00771 
00772 #endif /* BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP */

Generated on Wed May 12 21:30:08 2010 for Belos by  doxygen 1.4.7