Belos Package Browser (Single Doxygen Collection) Development
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 the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
00043 #define BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 
00052 #include "BelosLinearProblem.hpp"
00053 #include "BelosSolverManager.hpp"
00054 
00055 #include "BelosPseudoBlockCGIter.hpp"
00056 #include "BelosStatusTestMaxIters.hpp"
00057 #include "BelosStatusTestGenResNorm.hpp"
00058 #include "BelosStatusTestCombo.hpp"
00059 #include "BelosStatusTestOutputFactory.hpp"
00060 #include "BelosOutputManager.hpp"
00061 #include "Teuchos_BLAS.hpp"
00062 #include "Teuchos_TimeMonitor.hpp"
00063 
00080 namespace Belos {
00081   
00083 
00084   
00091   class PseudoBlockCGSolMgrLinearProblemFailure : public BelosError {public:
00092     PseudoBlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00093     {}};
00094   
00101   class PseudoBlockCGSolMgrOrthoFailure : public BelosError {public:
00102     PseudoBlockCGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00103     {}};
00104   
00105   template<class ScalarType, class MV, class OP>
00106   class PseudoBlockCGSolMgr : public SolverManager<ScalarType,MV,OP> {
00107     
00108   private:
00109     typedef MultiVecTraits<ScalarType,MV> MVT;
00110     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00111     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00112     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00113     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00114     
00115   public:
00116     
00118 
00119 
00125     PseudoBlockCGSolMgr();
00126     
00136     PseudoBlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00137                    const Teuchos::RCP<Teuchos::ParameterList> &pl );
00138     
00140     virtual ~PseudoBlockCGSolMgr() {};
00142     
00144 
00145     
00146     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00147       return *problem_;
00148     }
00149     
00152     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00153    
00156     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00157  
00163     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00164       return tuple(timerSolve_);
00165     }
00166 
00168     int getNumIters() const {
00169       return numIters_;
00170     }
00171   
00175     bool isLOADetected() const { return false; }
00176   
00178     
00180 
00181    
00183     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00184    
00186     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00187     
00189 
00191 
00192 
00196     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00198 
00200 
00201     
00219     ReturnType solve();
00220     
00222     
00225     
00227     std::string description() const;
00228     
00230     
00231   private:
00232 
00233     // Method to convert std::string to enumerated type for residual.
00234     Belos::ScaleType convertStringToScaleType( const std::string& scaleType ) {
00235       if (scaleType == "Norm of Initial Residual") {
00236         return Belos::NormOfInitRes;
00237       } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00238         return Belos::NormOfPrecInitRes;
00239       } else if (scaleType == "Norm of RHS") {
00240         return Belos::NormOfRHS;
00241       } else if (scaleType == "None") {
00242         return Belos::None;
00243       } else 
00244         TEST_FOR_EXCEPTION( true ,std::logic_error,
00245           "Belos::PseudoBlockCGSolMgr(): Invalid residual scaling type.");
00246     }
00247 
00248     // Linear problem.
00249     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00250     
00251     // Output manager.
00252     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00253     Teuchos::RCP<std::ostream> outputStream_;
00254 
00255     // Status test.
00256     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00257     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00258     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00259     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00260 
00261     // Orthogonalization manager.
00262     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00263 
00264      // Current parameter list.
00265     Teuchos::RCP<ParameterList> params_;
00266    
00267     // Default solver values.
00268     static const MagnitudeType convtol_default_;
00269     static const int maxIters_default_;
00270     static const bool showMaxResNormOnly_default_;
00271     static const int verbosity_default_;
00272     static const int outputStyle_default_;
00273     static const int outputFreq_default_;
00274     static const int defQuorum_default_;
00275     static const std::string resScale_default_; 
00276     static const std::string label_default_;
00277     static const Teuchos::RCP<std::ostream> outputStream_default_;
00278 
00279     // Current solver values.
00280     MagnitudeType convtol_;
00281     int maxIters_, numIters_;
00282     int verbosity_, outputStyle_, outputFreq_, defQuorum_;
00283     bool showMaxResNormOnly_;
00284     std::string resScale_;       
00285  
00286     // Timers.
00287     std::string label_;
00288     Teuchos::RCP<Teuchos::Time> timerSolve_;
00289 
00290     // Internal state variables.
00291     bool isSet_;
00292   };
00293   
00294   
00295 // Default solver values.
00296 template<class ScalarType, class MV, class OP>
00297 const typename PseudoBlockCGSolMgr<ScalarType,MV,OP>::MagnitudeType PseudoBlockCGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00298 
00299 template<class ScalarType, class MV, class OP>
00300 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00301 
00302 template<class ScalarType, class MV, class OP>
00303 const bool PseudoBlockCGSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00304 
00305 template<class ScalarType, class MV, class OP>
00306 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00307 
00308 template<class ScalarType, class MV, class OP>
00309 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00310 
00311 template<class ScalarType, class MV, class OP>
00312 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00313 
00314 template<class ScalarType, class MV, class OP>
00315 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::defQuorum_default_ = 1;
00316 
00317 template<class ScalarType, class MV, class OP>
00318 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::resScale_default_ = "Norm of Initial Residual";
00319 
00320 template<class ScalarType, class MV, class OP>
00321 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00322 
00323 template<class ScalarType, class MV, class OP>
00324 const Teuchos::RCP<std::ostream> PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00325 
00326 
00327 // Empty Constructor
00328 template<class ScalarType, class MV, class OP>
00329 PseudoBlockCGSolMgr<ScalarType,MV,OP>::PseudoBlockCGSolMgr() :
00330   outputStream_(outputStream_default_),
00331   convtol_(convtol_default_),
00332   maxIters_(maxIters_default_),
00333   verbosity_(verbosity_default_),
00334   outputStyle_(outputStyle_default_),
00335   outputFreq_(outputFreq_default_),
00336   defQuorum_(defQuorum_default_),
00337   showMaxResNormOnly_(showMaxResNormOnly_default_),
00338   resScale_(resScale_default_),
00339   label_(label_default_),
00340   isSet_(false)
00341 {}
00342 
00343 // Basic Constructor
00344 template<class ScalarType, class MV, class OP>
00345 PseudoBlockCGSolMgr<ScalarType,MV,OP>::PseudoBlockCGSolMgr( 
00346                const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00347                const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00348   problem_(problem),
00349   outputStream_(outputStream_default_),
00350   convtol_(convtol_default_),
00351   maxIters_(maxIters_default_),
00352   verbosity_(verbosity_default_),
00353   outputStyle_(outputStyle_default_),
00354   outputFreq_(outputFreq_default_),
00355   defQuorum_(defQuorum_default_),
00356   showMaxResNormOnly_(showMaxResNormOnly_default_),
00357   resScale_(resScale_default_),
00358   label_(label_default_),
00359   isSet_(false)
00360 {
00361   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00362 
00363   // If the parameter list pointer is null, then set the current parameters to the default parameter list.
00364   if (!is_null(pl)) {
00365     // Set the parameters using the list that was passed in.
00366     setParameters( pl );  
00367   }
00368 }
00369 
00370 template<class ScalarType, class MV, class OP>
00371 void PseudoBlockCGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00372 {
00373   // Create the internal parameter list if ones doesn't already exist.
00374   if (params_ == Teuchos::null) {
00375     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00376   }
00377   else {
00378     params->validateParameters(*getValidParameters());
00379   }
00380 
00381   // Check for maximum number of iterations
00382   if (params->isParameter("Maximum Iterations")) {
00383     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00384 
00385     // Update parameter in our list and in status test.
00386     params_->set("Maximum Iterations", maxIters_);
00387     if (maxIterTest_!=Teuchos::null)
00388       maxIterTest_->setMaxIters( maxIters_ );
00389   }
00390 
00391   // Check to see if the timer label changed.
00392   if (params->isParameter("Timer Label")) {
00393     std::string tempLabel = params->get("Timer Label", label_default_);
00394 
00395     // Update parameter in our list and solver timer
00396     if (tempLabel != label_) {
00397       label_ = tempLabel;
00398       params_->set("Timer Label", label_);
00399       std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00400       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00401     }
00402   }
00403 
00404   // Check for a change in verbosity level
00405   if (params->isParameter("Verbosity")) {
00406     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00407       verbosity_ = params->get("Verbosity", verbosity_default_);
00408     } else {
00409       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00410     }
00411 
00412     // Update parameter in our list.
00413     params_->set("Verbosity", verbosity_);
00414     if (printer_ != Teuchos::null)
00415       printer_->setVerbosity(verbosity_);
00416   }
00417 
00418   // Check for a change in output style
00419   if (params->isParameter("Output Style")) {
00420     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00421       outputStyle_ = params->get("Output Style", outputStyle_default_);
00422     } else {
00423       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00424     }
00425 
00426     // Reconstruct the convergence test if the explicit residual test is not being used.
00427     params_->set("Output Style", outputStyle_);
00428     outputTest_ = Teuchos::null;
00429   }
00430 
00431   // output stream
00432   if (params->isParameter("Output Stream")) {
00433     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00434 
00435     // Update parameter in our list.
00436     params_->set("Output Stream", outputStream_);
00437     if (printer_ != Teuchos::null)
00438       printer_->setOStream( outputStream_ );
00439   }
00440 
00441   // frequency level
00442   if (verbosity_ & Belos::StatusTestDetails) {
00443     if (params->isParameter("Output Frequency")) {
00444       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00445     }
00446 
00447     // Update parameter in out list and output status test.
00448     params_->set("Output Frequency", outputFreq_);
00449     if (outputTest_ != Teuchos::null)
00450       outputTest_->setOutputFrequency( outputFreq_ );
00451   }
00452 
00453   // Create output manager if we need to.
00454   if (printer_ == Teuchos::null) {
00455     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00456   }
00457   
00458   // Convergence
00459   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00460   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00461 
00462   // Check for convergence tolerance
00463   if (params->isParameter("Convergence Tolerance")) {
00464     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00465 
00466     // Update parameter in our list and residual tests.
00467     params_->set("Convergence Tolerance", convtol_);
00468     if (convTest_ != Teuchos::null)
00469       convTest_->setTolerance( convtol_ );
00470   }
00471 
00472   if (params->isParameter("Show Maximum Residual Norm Only")) {
00473     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00474 
00475     // Update parameter in our list and residual tests
00476     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00477     if (convTest_ != Teuchos::null)
00478       convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00479   }
00480 
00481   // Check for a change in scaling, if so we need to build new residual tests.
00482   bool newResTest = false;
00483   if (params->isParameter("Residual Scaling")) {
00484     std::string tempResScale = Teuchos::getParameter<std::string>( *params, "Residual Scaling" );
00485 
00486     // Only update the scaling if it's different.
00487     if (resScale_ != tempResScale) {
00488       Belos::ScaleType resScaleType = convertStringToScaleType( tempResScale );
00489       resScale_ = tempResScale;
00490 
00491       // Update parameter in our list and residual tests
00492       params_->set("Residual Scaling", resScale_);
00493       if (convTest_ != Teuchos::null) {
00494         try {
00495           convTest_->defineScaleForm( resScaleType, Belos::TwoNorm );
00496         }
00497         catch (std::exception& e) {
00498     // Make sure the convergence test gets constructed again.
00499     newResTest = true;
00500         }
00501       }
00502     }      
00503   }
00504 
00505   // Get the deflation quorum, or number of converged systems before deflation is allowed
00506   if (params->isParameter("Deflation Quorum")) {
00507     defQuorum_ = params->get("Deflation Quorum", defQuorum_);
00508     params_->set("Deflation Quorum", defQuorum_);
00509     if (convTest_ != Teuchos::null)
00510       convTest_->setQuorum( defQuorum_ );
00511   }
00512 
00513   // Create status tests if we need to.
00514 
00515   // Basic test checks maximum iterations and native residual.
00516   if (maxIterTest_ == Teuchos::null)
00517     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00518   
00519   // Implicit residual test, using the native residual to determine if convergence was achieved.
00520   if (convTest_ == Teuchos::null || newResTest) {
00521     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, defQuorum_, showMaxResNormOnly_ ) );
00522     convTest_->defineScaleForm( convertStringToScaleType( resScale_ ), Belos::TwoNorm );
00523   } 
00524 
00525   if (sTest_ == Teuchos::null || newResTest)
00526     sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00527   
00528   if (outputTest_ == Teuchos::null || newResTest) {
00529 
00530     // Create the status test output class.
00531     // This class manages and formats the output from the status test.
00532     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00533     outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00534 
00535     // Set the solver string for the output test
00536     std::string solverDesc = " Pseudo Block CG ";
00537     outputTest_->setSolverDesc( solverDesc );
00538 
00539   }
00540 
00541   // Create the timer if we need to.
00542   if (timerSolve_ == Teuchos::null) {
00543     std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00544     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00545   }
00546 
00547   // Inform the solver manager that the current parameters were set.
00548   isSet_ = true;
00549 }
00550 
00551     
00552 template<class ScalarType, class MV, class OP>
00553 Teuchos::RCP<const Teuchos::ParameterList>
00554 PseudoBlockCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00555 {
00556   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00557   if (is_null(validPL)) {
00558     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00559     // Set all the valid parameters and their default values.
00560     pl= Teuchos::rcp( new Teuchos::ParameterList() );
00561     pl->set("Convergence Tolerance", convtol_default_,
00562       "The relative residual tolerance that needs to be achieved by the\n"
00563       "iterative solver in order for the linera system to be declared converged.");
00564     pl->set("Maximum Iterations", maxIters_default_,
00565       "The maximum number of block iterations allowed for each\n"
00566       "set of RHS solved.");
00567     pl->set("Verbosity", verbosity_default_,
00568       "What type(s) of solver information should be outputted\n"
00569       "to the output stream.");
00570     pl->set("Output Style", outputStyle_default_,
00571       "What style is used for the solver information outputted\n"
00572       "to the output stream.");
00573     pl->set("Output Frequency", outputFreq_default_,
00574       "How often convergence information should be outputted\n"
00575       "to the output stream.");  
00576     pl->set("Deflation Quorum", defQuorum_default_,
00577       "The number of linear systems that need to converge before\n"
00578       "they are deflated.  This number should be <= block size.");
00579     pl->set("Output Stream", outputStream_default_,
00580       "A reference-counted pointer to the output stream where all\n"
00581       "solver output is sent.");
00582     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00583       "When convergence information is printed, only show the maximum\n"
00584       "relative residual norm when the block size is greater than one.");
00585     pl->set("Residual Scaling", resScale_default_,
00586       "The type of scaling used in the residual convergence test.");
00587     pl->set("Timer Label", label_default_,
00588       "The string to use as a prefix for the timer labels.");
00589     //  defaultParams_->set("Restart Timers", restartTimers_);
00590     validPL = pl;
00591   }
00592   return validPL;
00593 }
00594 
00595 
00596 // solve()
00597 template<class ScalarType, class MV, class OP>
00598 ReturnType PseudoBlockCGSolMgr<ScalarType,MV,OP>::solve() {
00599 
00600   // Set the current parameters if they were not set before.
00601   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00602   // then didn't set any parameters using setParameters().
00603   if (!isSet_) { setParameters( params_ ); }
00604   
00605   Teuchos::BLAS<int,ScalarType> blas;
00606   
00607   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockCGSolMgrLinearProblemFailure,
00608                      "Belos::PseudoBlockCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00609 
00610   // Create indices for the linear systems to be solved.
00611   int startPtr = 0;
00612   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00613   int numCurrRHS = numRHS2Solve;
00614 
00615   std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
00616   for (int i=0; i<numRHS2Solve; ++i) {
00617     currIdx[i] = startPtr+i; 
00618     currIdx2[i]=i;
00619   }
00620 
00621   // Inform the linear problem of the current linear system to solve.
00622   problem_->setLSIndex( currIdx );
00623 
00625   // Parameter list
00626   Teuchos::ParameterList plist;
00627 
00628   // Reset the status test.  
00629   outputTest_->reset();
00630 
00631   // Assume convergence is achieved, then let any failed convergence set this to false.
00632   bool isConverged = true;  
00633 
00635   // Pseudo-Block CG solver
00636 
00637   Teuchos::RCP<PseudoBlockCGIter<ScalarType,MV,OP> > block_cg_iter
00638     = Teuchos::rcp( new PseudoBlockCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );  
00639 
00640   // Enter solve() iterations
00641   {
00642     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00643 
00644     while ( numRHS2Solve > 0 ) {
00645 
00646       // Reset the active / converged vectors from this block
00647       std::vector<int> convRHSIdx;
00648       std::vector<int> currRHSIdx( currIdx );
00649       currRHSIdx.resize(numCurrRHS);
00650 
00651       // Reset the number of iterations.
00652       block_cg_iter->resetNumIters();
00653 
00654       // Reset the number of calls that the status test output knows about.
00655       outputTest_->resetNumCalls();
00656 
00657       // Get the current residual for this block of linear systems.
00658       Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
00659 
00660       // Get a new state struct and initialize the solver.
00661       CGIterationState<ScalarType,MV> newState;
00662       newState.R = R_0;
00663       block_cg_iter->initializeCG(newState);
00664 
00665       while(1) {
00666   
00667   // tell block_gmres_iter to iterate
00668   try {
00669     block_cg_iter->iterate();
00670     
00672     //
00673     // check convergence first
00674     //
00676     if ( convTest_->getStatus() == Passed ) {
00677       
00678       // Figure out which linear systems converged.
00679       std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
00680 
00681       // If the number of converged linear systems is equal to the
00682             // number of current linear systems, then we are done with this block.
00683       if (convIdx.size() == currRHSIdx.size())
00684         break;  // break from while(1){block_cg_iter->iterate()}
00685 
00686       // Inform the linear problem that we are finished with this current linear system.
00687       problem_->setCurrLS();
00688 
00689       // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
00690       int have = 0;
00691             std::vector<int> unconvIdx(currRHSIdx.size());
00692       for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
00693         bool found = false;
00694         for (unsigned int j=0; j<convIdx.size(); ++j) {
00695     if (currRHSIdx[i] == convIdx[j]) {
00696       found = true;
00697       break;
00698     }
00699         }
00700         if (!found) {
00701                 currIdx2[have] = currIdx2[i];
00702     currRHSIdx[have++] = currRHSIdx[i];
00703         }
00704       }
00705       currRHSIdx.resize(have);
00706       currIdx2.resize(have);
00707       
00708       // Set the remaining indices after deflation.
00709       problem_->setLSIndex( currRHSIdx );
00710       
00711       // Get the current residual vector.
00712       std::vector<MagnitudeType> norms;
00713             R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
00714       for (int i=0; i<have; ++i) { currIdx2[i] = i; }
00715 
00716       // Set the new state and initialize the solver.
00717       CGIterationState<ScalarType,MV> defstate;
00718       defstate.R = R_0;
00719       block_cg_iter->initializeCG(defstate);
00720     }
00721 
00723     //
00724     // check for maximum iterations
00725     //
00727     else if ( maxIterTest_->getStatus() == Passed ) {
00728       // we don't have convergence
00729       isConverged = false;
00730       break;  // break from while(1){block_cg_iter->iterate()}
00731     }
00732 
00734     //
00735     // we returned from iterate(), but none of our status tests Passed.
00736     // something is wrong, and it is probably our fault.
00737     //
00739 
00740     else {
00741       TEST_FOR_EXCEPTION(true,std::logic_error,
00742              "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
00743     }
00744   } 
00745   catch (const std::exception &e) {
00746     printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration " 
00747            << block_cg_iter->getNumIters() << std::endl 
00748            << e.what() << std::endl;
00749     throw;
00750   }
00751       }
00752       
00753       // Inform the linear problem that we are finished with this block linear system.
00754       problem_->setCurrLS();
00755       
00756       // Update indices for the linear systems to be solved.
00757       startPtr += numCurrRHS;
00758       numRHS2Solve -= numCurrRHS;
00759 
00760       if ( numRHS2Solve > 0 ) { 
00761 
00762   numCurrRHS = numRHS2Solve;  
00763   currIdx.resize( numCurrRHS );
00764   currIdx2.resize( numCurrRHS );
00765   for (int i=0; i<numCurrRHS; ++i) 
00766     { currIdx[i] = startPtr+i; currIdx2[i] = i; }
00767   
00768   // Set the next indices.
00769   problem_->setLSIndex( currIdx );  
00770       }
00771       else {
00772   currIdx.resize( numRHS2Solve );
00773       }
00774       
00775     }// while ( numRHS2Solve > 0 )
00776     
00777   }
00778                      
00779   // print final summary
00780   sTest_->print( printer_->stream(FinalSummary) );
00781  
00782   // print timing information
00783   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
00784  
00785   // get iteration information for this solve
00786   numIters_ = maxIterTest_->getNumIters();
00787  
00788   if (!isConverged ) {
00789     return Unconverged; // return from PseudoBlockCGSolMgr::solve() 
00790   }
00791   return Converged; // return from PseudoBlockCGSolMgr::solve() 
00792 }
00793 
00794 //  This method requires the solver manager to return a std::string that describes itself.
00795 template<class ScalarType, class MV, class OP>
00796 std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::description() const
00797 {
00798   std::ostringstream oss;
00799   oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00800   oss << "{";
00801   oss << "}";
00802   return oss.str();
00803 }
00804   
00805 } // end Belos namespace
00806 
00807 #endif /* BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines