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 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00063 #include "Teuchos_TimeMonitor.hpp"
00064 #endif
00065 
00082 namespace Belos {
00083   
00085 
00086   
00093   class PseudoBlockCGSolMgrLinearProblemFailure : public BelosError {public:
00094     PseudoBlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00095     {}};
00096   
00103   class PseudoBlockCGSolMgrOrthoFailure : public BelosError {public:
00104     PseudoBlockCGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00105     {}};
00106   
00107   template<class ScalarType, class MV, class OP>
00108   class PseudoBlockCGSolMgr : public SolverManager<ScalarType,MV,OP> {
00109     
00110   private:
00111     typedef MultiVecTraits<ScalarType,MV> MVT;
00112     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00113     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00114     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00115     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00116     
00117   public:
00118     
00120 
00121 
00127     PseudoBlockCGSolMgr();
00128     
00138     PseudoBlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00139                    const Teuchos::RCP<Teuchos::ParameterList> &pl );
00140     
00142     virtual ~PseudoBlockCGSolMgr() {};
00144     
00146 
00147     
00148     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00149       return *problem_;
00150     }
00151     
00154     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00155    
00158     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00159  
00165     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00166       return Teuchos::tuple(timerSolve_);
00167     }
00168 
00170     int getNumIters() const {
00171       return numIters_;
00172     }
00173   
00177     bool isLOADetected() const { return false; }
00178   
00180     
00182 
00183    
00185     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00186    
00188     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00189     
00191 
00193 
00194 
00198     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00200 
00202 
00203     
00221     ReturnType solve();
00222     
00224     
00227     
00229     std::string description() const;
00230     
00232     
00233   private:
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<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00246     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00247 
00248     // Orthogonalization manager.
00249     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00250 
00251     // Current parameter list.
00252     Teuchos::RCP<Teuchos::ParameterList> params_;
00253 
00259     mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
00260    
00261     // Default solver values.
00262     static const MagnitudeType convtol_default_;
00263     static const int maxIters_default_;
00264     static const bool assertPositiveDefiniteness_default_;
00265     static const bool showMaxResNormOnly_default_;
00266     static const int verbosity_default_;
00267     static const int outputStyle_default_;
00268     static const int outputFreq_default_;
00269     static const int defQuorum_default_;
00270     static const std::string resScale_default_; 
00271     static const std::string label_default_;
00272     static const Teuchos::RCP<std::ostream> outputStream_default_;
00273 
00274     // Current solver values.
00275     MagnitudeType convtol_;
00276     int maxIters_, numIters_;
00277     int verbosity_, outputStyle_, outputFreq_, defQuorum_;
00278     bool assertPositiveDefiniteness_, showMaxResNormOnly_;
00279     std::string resScale_;       
00280  
00281     // Timers.
00282     std::string label_;
00283     Teuchos::RCP<Teuchos::Time> timerSolve_;
00284 
00285     // Internal state variables.
00286     bool isSet_;
00287   };
00288   
00289   
00290 // Default solver values.
00291 template<class ScalarType, class MV, class OP>
00292 const typename PseudoBlockCGSolMgr<ScalarType,MV,OP>::MagnitudeType PseudoBlockCGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00293 
00294 template<class ScalarType, class MV, class OP>
00295 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00296 
00297 template<class ScalarType, class MV, class OP>
00298 const bool PseudoBlockCGSolMgr<ScalarType,MV,OP>::assertPositiveDefiniteness_default_ = true;
00299 
00300 template<class ScalarType, class MV, class OP>
00301 const bool PseudoBlockCGSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00302 
00303 template<class ScalarType, class MV, class OP>
00304 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00305 
00306 template<class ScalarType, class MV, class OP>
00307 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00308 
00309 template<class ScalarType, class MV, class OP>
00310 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00311 
00312 template<class ScalarType, class MV, class OP>
00313 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::defQuorum_default_ = 1;
00314 
00315 template<class ScalarType, class MV, class OP>
00316 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::resScale_default_ = "Norm of Initial Residual";
00317 
00318 template<class ScalarType, class MV, class OP>
00319 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00320 
00321 template<class ScalarType, class MV, class OP>
00322 const Teuchos::RCP<std::ostream> PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00323 
00324 
00325 // Empty Constructor
00326 template<class ScalarType, class MV, class OP>
00327 PseudoBlockCGSolMgr<ScalarType,MV,OP>::PseudoBlockCGSolMgr() :
00328   outputStream_(outputStream_default_),
00329   convtol_(convtol_default_),
00330   maxIters_(maxIters_default_),
00331   verbosity_(verbosity_default_),
00332   outputStyle_(outputStyle_default_),
00333   outputFreq_(outputFreq_default_),
00334   defQuorum_(defQuorum_default_),
00335   assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
00336   showMaxResNormOnly_(showMaxResNormOnly_default_),
00337   resScale_(resScale_default_),
00338   label_(label_default_),
00339   isSet_(false)
00340 {}
00341 
00342 // Basic Constructor
00343 template<class ScalarType, class MV, class OP>
00344 PseudoBlockCGSolMgr<ScalarType,MV,OP>::PseudoBlockCGSolMgr( 
00345                const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00346                const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00347   problem_(problem),
00348   outputStream_(outputStream_default_),
00349   convtol_(convtol_default_),
00350   maxIters_(maxIters_default_),
00351   verbosity_(verbosity_default_),
00352   outputStyle_(outputStyle_default_),
00353   outputFreq_(outputFreq_default_),
00354   defQuorum_(defQuorum_default_),
00355   assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
00356   showMaxResNormOnly_(showMaxResNormOnly_default_),
00357   resScale_(resScale_default_),
00358   label_(label_default_),
00359   isSet_(false)
00360 {
00361   TEUCHOS_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   using Teuchos::ParameterList;
00374   using Teuchos::parameterList;
00375   using Teuchos::RCP;
00376 
00377   RCP<const ParameterList> defaultParams = getValidParameters();
00378 
00379   // Create the internal parameter list if one doesn't already exist.
00380   if (params_.is_null()) {
00381     params_ = parameterList (*defaultParams);
00382   } else {
00383     params->validateParameters (*defaultParams);
00384   }
00385 
00386   // Check for maximum number of iterations
00387   if (params->isParameter("Maximum Iterations")) {
00388     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00389 
00390     // Update parameter in our list and in status test.
00391     params_->set("Maximum Iterations", maxIters_);
00392     if (maxIterTest_!=Teuchos::null)
00393       maxIterTest_->setMaxIters( maxIters_ );
00394   }
00395 
00396   // Check if positive definiteness assertions are to be performed
00397   if (params->isParameter("Assert Positive Definiteness")) {
00398     assertPositiveDefiniteness_ = params->get("Assert Positive Definiteness",assertPositiveDefiniteness_default_);
00399 
00400     // Update parameter in our list.
00401     params_->set("Assert Positive Definiteness", assertPositiveDefiniteness_);
00402   }
00403 
00404   // Check to see if the timer label changed.
00405   if (params->isParameter("Timer Label")) {
00406     std::string tempLabel = params->get("Timer Label", label_default_);
00407 
00408     // Update parameter in our list and solver timer
00409     if (tempLabel != label_) {
00410       label_ = tempLabel;
00411       params_->set("Timer Label", label_);
00412       std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00413 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00414       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00415 #endif
00416     }
00417   }
00418 
00419   // Check for a change in verbosity level
00420   if (params->isParameter("Verbosity")) {
00421     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00422       verbosity_ = params->get("Verbosity", verbosity_default_);
00423     } else {
00424       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00425     }
00426 
00427     // Update parameter in our list.
00428     params_->set("Verbosity", verbosity_);
00429     if (printer_ != Teuchos::null)
00430       printer_->setVerbosity(verbosity_);
00431   }
00432 
00433   // Check for a change in output style
00434   if (params->isParameter("Output Style")) {
00435     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00436       outputStyle_ = params->get("Output Style", outputStyle_default_);
00437     } else {
00438       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00439     }
00440 
00441     // Reconstruct the convergence test if the explicit residual test is not being used.
00442     params_->set("Output Style", outputStyle_);
00443     outputTest_ = Teuchos::null;
00444   }
00445 
00446   // output stream
00447   if (params->isParameter("Output Stream")) {
00448     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00449 
00450     // Update parameter in our list.
00451     params_->set("Output Stream", outputStream_);
00452     if (printer_ != Teuchos::null)
00453       printer_->setOStream( outputStream_ );
00454   }
00455 
00456   // frequency level
00457   if (verbosity_ & Belos::StatusTestDetails) {
00458     if (params->isParameter("Output Frequency")) {
00459       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00460     }
00461 
00462     // Update parameter in out list and output status test.
00463     params_->set("Output Frequency", outputFreq_);
00464     if (outputTest_ != Teuchos::null)
00465       outputTest_->setOutputFrequency( outputFreq_ );
00466   }
00467 
00468   // Create output manager if we need to.
00469   if (printer_ == Teuchos::null) {
00470     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00471   }
00472   
00473   // Convergence
00474   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00475   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00476 
00477   // Check for convergence tolerance
00478   if (params->isParameter("Convergence Tolerance")) {
00479     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00480 
00481     // Update parameter in our list and residual tests.
00482     params_->set("Convergence Tolerance", convtol_);
00483     if (convTest_ != Teuchos::null)
00484       convTest_->setTolerance( convtol_ );
00485   }
00486 
00487   if (params->isParameter("Show Maximum Residual Norm Only")) {
00488     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00489 
00490     // Update parameter in our list and residual tests
00491     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00492     if (convTest_ != Teuchos::null)
00493       convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00494   }
00495 
00496   // Check for a change in scaling, if so we need to build new residual tests.
00497   bool newResTest = false;
00498   {
00499     // "Residual Scaling" is the old parameter name; "Implicit
00500     // Residual Scaling" is the new name.  We support both options for
00501     // backwards compatibility.
00502     std::string tempResScale = resScale_;
00503     bool implicitResidualScalingName = false;
00504     if (params->isParameter ("Residual Scaling")) {
00505       tempResScale = params->get<std::string> ("Residual Scaling");
00506     } 
00507     else if (params->isParameter ("Implicit Residual Scaling")) {
00508       tempResScale = params->get<std::string> ("Implicit Residual Scaling");
00509       implicitResidualScalingName = true;
00510     }
00511 
00512     // Only update the scaling if it's different.
00513     if (resScale_ != tempResScale) {
00514       Belos::ScaleType resScaleType = convertStringToScaleType( tempResScale );
00515       resScale_ = tempResScale;
00516 
00517       // Update parameter in our list and residual tests, using the
00518       // given parameter name.
00519       if (implicitResidualScalingName) {
00520   params_->set ("Implicit Residual Scaling", resScale_);
00521       } 
00522       else {
00523   params_->set ("Residual Scaling", resScale_); 
00524       }
00525 
00526       if (! convTest_.is_null()) {
00527         try {
00528           convTest_->defineScaleForm( resScaleType, Belos::TwoNorm );
00529         }
00530         catch (std::exception& e) {
00531     // Make sure the convergence test gets constructed again.
00532     newResTest = true;
00533         }
00534       }
00535     }      
00536   }
00537 
00538   // Get the deflation quorum, or number of converged systems before deflation is allowed
00539   if (params->isParameter("Deflation Quorum")) {
00540     defQuorum_ = params->get("Deflation Quorum", defQuorum_);
00541     params_->set("Deflation Quorum", defQuorum_);
00542     if (convTest_ != Teuchos::null)
00543       convTest_->setQuorum( defQuorum_ );
00544   }
00545 
00546   // Create status tests if we need to.
00547 
00548   // Basic test checks maximum iterations and native residual.
00549   if (maxIterTest_ == Teuchos::null)
00550     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00551   
00552   // Implicit residual test, using the native residual to determine if convergence was achieved.
00553   if (convTest_ == Teuchos::null || newResTest) {
00554     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, defQuorum_, showMaxResNormOnly_ ) );
00555     convTest_->defineScaleForm( convertStringToScaleType( resScale_ ), Belos::TwoNorm );
00556   } 
00557 
00558   if (sTest_ == Teuchos::null || newResTest)
00559     sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00560   
00561   if (outputTest_ == Teuchos::null || newResTest) {
00562 
00563     // Create the status test output class.
00564     // This class manages and formats the output from the status test.
00565     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00566     outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00567 
00568     // Set the solver string for the output test
00569     std::string solverDesc = " Pseudo Block CG ";
00570     outputTest_->setSolverDesc( solverDesc );
00571 
00572   }
00573 
00574   // Create the timer if we need to.
00575   if (timerSolve_ == Teuchos::null) {
00576     std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00577 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00578     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00579 #endif
00580   }
00581 
00582   // Inform the solver manager that the current parameters were set.
00583   isSet_ = true;
00584 }
00585 
00586     
00587 template<class ScalarType, class MV, class OP>
00588 Teuchos::RCP<const Teuchos::ParameterList>
00589 PseudoBlockCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00590 {
00591   using Teuchos::ParameterList;
00592   using Teuchos::parameterList;
00593   using Teuchos::RCP;
00594 
00595   if (validParams_.is_null()) {
00596     // Set all the valid parameters and their default values.
00597     RCP<ParameterList> pl = parameterList ();
00598     pl->set("Convergence Tolerance", convtol_default_,
00599       "The relative residual tolerance that needs to be achieved by the\n"
00600       "iterative solver in order for the linera system to be declared converged.");
00601     pl->set("Maximum Iterations", maxIters_default_,
00602       "The maximum number of block iterations allowed for each\n"
00603       "set of RHS solved.");
00604     pl->set("Assert Positive Definiteness", assertPositiveDefiniteness_default_,
00605       "Whether or not to assert that the linear operator\n"
00606       "and the preconditioner are indeed positive definite.");
00607     pl->set("Verbosity", verbosity_default_,
00608       "What type(s) of solver information should be outputted\n"
00609       "to the output stream.");
00610     pl->set("Output Style", outputStyle_default_,
00611       "What style is used for the solver information outputted\n"
00612       "to the output stream.");
00613     pl->set("Output Frequency", outputFreq_default_,
00614       "How often convergence information should be outputted\n"
00615       "to the output stream.");  
00616     pl->set("Deflation Quorum", defQuorum_default_,
00617       "The number of linear systems that need to converge before\n"
00618       "they are deflated.  This number should be <= block size.");
00619     pl->set("Output Stream", outputStream_default_,
00620       "A reference-counted pointer to the output stream where all\n"
00621       "solver output is sent.");
00622     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00623       "When convergence information is printed, only show the maximum\n"
00624       "relative residual norm when the block size is greater than one.");
00625     pl->set("Implicit Residual Scaling", resScale_default_,
00626       "The type of scaling used in the residual convergence test.");
00627     // We leave the old name as a valid parameter for backwards
00628     // compatibility (so that validateParametersAndSetDefaults()
00629     // doesn't raise an exception if it encounters "Residual
00630     // Scaling").  The new name was added for compatibility with other
00631     // solvers, none of which use "Residual Scaling".
00632     pl->set("Residual Scaling", resScale_default_,
00633       "The type of scaling used in the residual convergence test.  This "
00634       "name is deprecated; the new name is \"Implicit Residual Scaling\".");
00635     pl->set("Timer Label", label_default_,
00636       "The string to use as a prefix for the timer labels.");
00637     //  defaultParams_->set("Restart Timers", restartTimers_);
00638     validParams_ = pl;
00639   }
00640   return validParams_;
00641 }
00642 
00643 
00644 // solve()
00645 template<class ScalarType, class MV, class OP>
00646 ReturnType PseudoBlockCGSolMgr<ScalarType,MV,OP>::solve() {
00647 
00648   // Set the current parameters if they were not set before.
00649   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00650   // then didn't set any parameters using setParameters().
00651   if (!isSet_) { setParameters( params_ ); }
00652   
00653   Teuchos::BLAS<int,ScalarType> blas;
00654   
00655   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockCGSolMgrLinearProblemFailure,
00656                      "Belos::PseudoBlockCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00657 
00658   // Create indices for the linear systems to be solved.
00659   int startPtr = 0;
00660   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00661   int numCurrRHS = numRHS2Solve;
00662 
00663   std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
00664   for (int i=0; i<numRHS2Solve; ++i) {
00665     currIdx[i] = startPtr+i; 
00666     currIdx2[i]=i;
00667   }
00668 
00669   // Inform the linear problem of the current linear system to solve.
00670   problem_->setLSIndex( currIdx );
00671 
00673   // Parameter list
00674   Teuchos::ParameterList plist;
00675 
00676   plist.set("Assert Positive Definiteness",assertPositiveDefiniteness_);
00677 
00678   // Reset the status test.  
00679   outputTest_->reset();
00680 
00681   // Assume convergence is achieved, then let any failed convergence set this to false.
00682   bool isConverged = true;  
00683 
00685   // Pseudo-Block CG solver
00686 
00687   Teuchos::RCP<PseudoBlockCGIter<ScalarType,MV,OP> > block_cg_iter
00688     = Teuchos::rcp( new PseudoBlockCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );  
00689 
00690   // Enter solve() iterations
00691   {
00692 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00693     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00694 #endif
00695 
00696     while ( numRHS2Solve > 0 ) {
00697 
00698       // Reset the active / converged vectors from this block
00699       std::vector<int> convRHSIdx;
00700       std::vector<int> currRHSIdx( currIdx );
00701       currRHSIdx.resize(numCurrRHS);
00702 
00703       // Reset the number of iterations.
00704       block_cg_iter->resetNumIters();
00705 
00706       // Reset the number of calls that the status test output knows about.
00707       outputTest_->resetNumCalls();
00708 
00709       // Get the current residual for this block of linear systems.
00710       Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
00711 
00712       // Get a new state struct and initialize the solver.
00713       CGIterationState<ScalarType,MV> newState;
00714       newState.R = R_0;
00715       block_cg_iter->initializeCG(newState);
00716 
00717       while(1) {
00718   
00719   // tell block_gmres_iter to iterate
00720   try {
00721     block_cg_iter->iterate();
00722     
00724     //
00725     // check convergence first
00726     //
00728     if ( convTest_->getStatus() == Passed ) {
00729       
00730       // Figure out which linear systems converged.
00731       std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
00732 
00733       // If the number of converged linear systems is equal to the
00734             // number of current linear systems, then we are done with this block.
00735       if (convIdx.size() == currRHSIdx.size())
00736         break;  // break from while(1){block_cg_iter->iterate()}
00737 
00738       // Inform the linear problem that we are finished with this current linear system.
00739       problem_->setCurrLS();
00740 
00741       // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
00742       int have = 0;
00743             std::vector<int> unconvIdx(currRHSIdx.size());
00744       for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
00745         bool found = false;
00746         for (unsigned int j=0; j<convIdx.size(); ++j) {
00747     if (currRHSIdx[i] == convIdx[j]) {
00748       found = true;
00749       break;
00750     }
00751         }
00752         if (!found) {
00753                 currIdx2[have] = currIdx2[i];
00754     currRHSIdx[have++] = currRHSIdx[i];
00755         }
00756       }
00757       currRHSIdx.resize(have);
00758       currIdx2.resize(have);
00759       
00760       // Set the remaining indices after deflation.
00761       problem_->setLSIndex( currRHSIdx );
00762       
00763       // Get the current residual vector.
00764       std::vector<MagnitudeType> norms;
00765             R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
00766       for (int i=0; i<have; ++i) { currIdx2[i] = i; }
00767 
00768       // Set the new state and initialize the solver.
00769       CGIterationState<ScalarType,MV> defstate;
00770       defstate.R = R_0;
00771       block_cg_iter->initializeCG(defstate);
00772     }
00773 
00775     //
00776     // check for maximum iterations
00777     //
00779     else if ( maxIterTest_->getStatus() == Passed ) {
00780       // we don't have convergence
00781       isConverged = false;
00782       break;  // break from while(1){block_cg_iter->iterate()}
00783     }
00784 
00786     //
00787     // we returned from iterate(), but none of our status tests Passed.
00788     // something is wrong, and it is probably our fault.
00789     //
00791 
00792     else {
00793       TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00794              "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
00795     }
00796   } 
00797   catch (const std::exception &e) {
00798     printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration " 
00799            << block_cg_iter->getNumIters() << std::endl 
00800            << e.what() << std::endl;
00801     throw;
00802   }
00803       }
00804       
00805       // Inform the linear problem that we are finished with this block linear system.
00806       problem_->setCurrLS();
00807       
00808       // Update indices for the linear systems to be solved.
00809       startPtr += numCurrRHS;
00810       numRHS2Solve -= numCurrRHS;
00811 
00812       if ( numRHS2Solve > 0 ) { 
00813 
00814   numCurrRHS = numRHS2Solve;  
00815   currIdx.resize( numCurrRHS );
00816   currIdx2.resize( numCurrRHS );
00817   for (int i=0; i<numCurrRHS; ++i) 
00818     { currIdx[i] = startPtr+i; currIdx2[i] = i; }
00819   
00820   // Set the next indices.
00821   problem_->setLSIndex( currIdx );  
00822       }
00823       else {
00824   currIdx.resize( numRHS2Solve );
00825       }
00826       
00827     }// while ( numRHS2Solve > 0 )
00828     
00829   }
00830                      
00831   // print final summary
00832   sTest_->print( printer_->stream(FinalSummary) );
00833  
00834   // print timing information
00835 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00836   // Calling summarize() can be expensive, so don't call unless the
00837   // user wants to print out timing details.  summarize() will do all
00838   // the work even if it's passed a "black hole" output stream.
00839   if (verbosity_ & TimingDetails)
00840     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
00841 #endif
00842  
00843   // get iteration information for this solve
00844   numIters_ = maxIterTest_->getNumIters();
00845  
00846   if (!isConverged ) {
00847     return Unconverged; // return from PseudoBlockCGSolMgr::solve() 
00848   }
00849   return Converged; // return from PseudoBlockCGSolMgr::solve() 
00850 }
00851 
00852 //  This method requires the solver manager to return a std::string that describes itself.
00853 template<class ScalarType, class MV, class OP>
00854 std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::description() const
00855 {
00856   std::ostringstream oss;
00857   oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00858   oss << "{";
00859   oss << "}";
00860   return oss.str();
00861 }
00862   
00863 } // end Belos namespace
00864 
00865 #endif /* BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines