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   numIters_(0),
00332   verbosity_(verbosity_default_),
00333   outputStyle_(outputStyle_default_),
00334   outputFreq_(outputFreq_default_),
00335   defQuorum_(defQuorum_default_),
00336   assertPositiveDefiniteness_(assertPositiveDefiniteness_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>::
00346 PseudoBlockCGSolMgr (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   numIters_(0),
00353   verbosity_(verbosity_default_),
00354   outputStyle_(outputStyle_default_),
00355   outputFreq_(outputFreq_default_),
00356   defQuorum_(defQuorum_default_),
00357   assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
00358   showMaxResNormOnly_(showMaxResNormOnly_default_),
00359   resScale_(resScale_default_),
00360   label_(label_default_),
00361   isSet_(false)
00362 {
00363   TEUCHOS_TEST_FOR_EXCEPTION(
00364     problem_.is_null (), std::invalid_argument,
00365     "Belos::PseudoBlockCGSolMgr two-argument constructor: "
00366     "'problem' is null.  You must supply a non-null Belos::LinearProblem "
00367     "instance when calling this constructor.");
00368 
00369   if (! pl.is_null ()) {
00370     // Set the parameters using the list that was passed in.
00371     setParameters (pl);
00372   }
00373 }
00374 
00375 template<class ScalarType, class MV, class OP>
00376 void PseudoBlockCGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00377 {
00378   using Teuchos::ParameterList;
00379   using Teuchos::parameterList;
00380   using Teuchos::RCP;
00381 
00382   RCP<const ParameterList> defaultParams = getValidParameters();
00383 
00384   // Create the internal parameter list if one doesn't already exist.
00385   if (params_.is_null()) {
00386     params_ = parameterList (*defaultParams);
00387   } else {
00388     params->validateParameters (*defaultParams);
00389   }
00390 
00391   // Check for maximum number of iterations
00392   if (params->isParameter("Maximum Iterations")) {
00393     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00394 
00395     // Update parameter in our list and in status test.
00396     params_->set("Maximum Iterations", maxIters_);
00397     if (maxIterTest_!=Teuchos::null)
00398       maxIterTest_->setMaxIters( maxIters_ );
00399   }
00400 
00401   // Check if positive definiteness assertions are to be performed
00402   if (params->isParameter("Assert Positive Definiteness")) {
00403     assertPositiveDefiniteness_ = params->get("Assert Positive Definiteness",assertPositiveDefiniteness_default_);
00404 
00405     // Update parameter in our list.
00406     params_->set("Assert Positive Definiteness", assertPositiveDefiniteness_);
00407   }
00408 
00409   // Check to see if the timer label changed.
00410   if (params->isParameter("Timer Label")) {
00411     std::string tempLabel = params->get("Timer Label", label_default_);
00412 
00413     // Update parameter in our list and solver timer
00414     if (tempLabel != label_) {
00415       label_ = tempLabel;
00416       params_->set("Timer Label", label_);
00417       std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00418 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00419       timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
00420 #endif
00421       if (ortho_ != Teuchos::null) {
00422         ortho_->setLabel( label_ );
00423       }
00424     }
00425   }
00426 
00427   // Check for a change in verbosity level
00428   if (params->isParameter("Verbosity")) {
00429     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00430       verbosity_ = params->get("Verbosity", verbosity_default_);
00431     } else {
00432       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00433     }
00434 
00435     // Update parameter in our list.
00436     params_->set("Verbosity", verbosity_);
00437     if (printer_ != Teuchos::null)
00438       printer_->setVerbosity(verbosity_);
00439   }
00440 
00441   // Check for a change in output style
00442   if (params->isParameter("Output Style")) {
00443     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00444       outputStyle_ = params->get("Output Style", outputStyle_default_);
00445     } else {
00446       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00447     }
00448 
00449     // Reconstruct the convergence test if the explicit residual test is not being used.
00450     params_->set("Output Style", outputStyle_);
00451     outputTest_ = Teuchos::null;
00452   }
00453 
00454   // output stream
00455   if (params->isParameter("Output Stream")) {
00456     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00457 
00458     // Update parameter in our list.
00459     params_->set("Output Stream", outputStream_);
00460     if (printer_ != Teuchos::null)
00461       printer_->setOStream( outputStream_ );
00462   }
00463 
00464   // frequency level
00465   if (verbosity_ & Belos::StatusTestDetails) {
00466     if (params->isParameter("Output Frequency")) {
00467       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00468     }
00469 
00470     // Update parameter in out list and output status test.
00471     params_->set("Output Frequency", outputFreq_);
00472     if (outputTest_ != Teuchos::null)
00473       outputTest_->setOutputFrequency( outputFreq_ );
00474   }
00475 
00476   // Create output manager if we need to.
00477   if (printer_ == Teuchos::null) {
00478     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00479   }
00480 
00481   // Convergence
00482   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00483   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00484 
00485   // Check for convergence tolerance
00486   if (params->isParameter("Convergence Tolerance")) {
00487     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00488 
00489     // Update parameter in our list and residual tests.
00490     params_->set("Convergence Tolerance", convtol_);
00491     if (convTest_ != Teuchos::null)
00492       convTest_->setTolerance( convtol_ );
00493   }
00494 
00495   if (params->isParameter("Show Maximum Residual Norm Only")) {
00496     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00497 
00498     // Update parameter in our list and residual tests
00499     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00500     if (convTest_ != Teuchos::null)
00501       convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00502   }
00503 
00504   // Check for a change in scaling, if so we need to build new residual tests.
00505   bool newResTest = false;
00506   {
00507     // "Residual Scaling" is the old parameter name; "Implicit
00508     // Residual Scaling" is the new name.  We support both options for
00509     // backwards compatibility.
00510     std::string tempResScale = resScale_;
00511     bool implicitResidualScalingName = false;
00512     if (params->isParameter ("Residual Scaling")) {
00513       tempResScale = params->get<std::string> ("Residual Scaling");
00514     }
00515     else if (params->isParameter ("Implicit Residual Scaling")) {
00516       tempResScale = params->get<std::string> ("Implicit Residual Scaling");
00517       implicitResidualScalingName = true;
00518     }
00519 
00520     // Only update the scaling if it's different.
00521     if (resScale_ != tempResScale) {
00522       Belos::ScaleType resScaleType = convertStringToScaleType( tempResScale );
00523       resScale_ = tempResScale;
00524 
00525       // Update parameter in our list and residual tests, using the
00526       // given parameter name.
00527       if (implicitResidualScalingName) {
00528         params_->set ("Implicit Residual Scaling", resScale_);
00529       }
00530       else {
00531         params_->set ("Residual Scaling", resScale_);
00532       }
00533 
00534       if (! convTest_.is_null()) {
00535         try {
00536           convTest_->defineScaleForm( resScaleType, Belos::TwoNorm );
00537         }
00538         catch (std::exception& e) {
00539           // Make sure the convergence test gets constructed again.
00540           newResTest = true;
00541         }
00542       }
00543     }
00544   }
00545 
00546   // Get the deflation quorum, or number of converged systems before deflation is allowed
00547   if (params->isParameter("Deflation Quorum")) {
00548     defQuorum_ = params->get("Deflation Quorum", defQuorum_);
00549     params_->set("Deflation Quorum", defQuorum_);
00550     if (convTest_ != Teuchos::null)
00551       convTest_->setQuorum( defQuorum_ );
00552   }
00553 
00554   // Create status tests if we need to.
00555 
00556   // Basic test checks maximum iterations and native residual.
00557   if (maxIterTest_ == Teuchos::null)
00558     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00559 
00560   // Implicit residual test, using the native residual to determine if convergence was achieved.
00561   if (convTest_ == Teuchos::null || newResTest) {
00562     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, defQuorum_, showMaxResNormOnly_ ) );
00563     convTest_->defineScaleForm( convertStringToScaleType( resScale_ ), Belos::TwoNorm );
00564   }
00565 
00566   if (sTest_ == Teuchos::null || newResTest)
00567     sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00568 
00569   if (outputTest_ == Teuchos::null || newResTest) {
00570 
00571     // Create the status test output class.
00572     // This class manages and formats the output from the status test.
00573     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00574     outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00575 
00576     // Set the solver string for the output test
00577     std::string solverDesc = " Pseudo Block CG ";
00578     outputTest_->setSolverDesc( solverDesc );
00579 
00580   }
00581 
00582   // Create the timer if we need to.
00583   if (timerSolve_ == Teuchos::null) {
00584     std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00585 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00586     timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
00587 #endif
00588   }
00589 
00590   // Inform the solver manager that the current parameters were set.
00591   isSet_ = true;
00592 }
00593 
00594 
00595 template<class ScalarType, class MV, class OP>
00596 Teuchos::RCP<const Teuchos::ParameterList>
00597 PseudoBlockCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00598 {
00599   using Teuchos::ParameterList;
00600   using Teuchos::parameterList;
00601   using Teuchos::RCP;
00602 
00603   if (validParams_.is_null()) {
00604     // Set all the valid parameters and their default values.
00605     RCP<ParameterList> pl = parameterList ();
00606     pl->set("Convergence Tolerance", convtol_default_,
00607       "The relative residual tolerance that needs to be achieved by the\n"
00608       "iterative solver in order for the linera system to be declared converged.");
00609     pl->set("Maximum Iterations", maxIters_default_,
00610       "The maximum number of block iterations allowed for each\n"
00611       "set of RHS solved.");
00612     pl->set("Assert Positive Definiteness", assertPositiveDefiniteness_default_,
00613       "Whether or not to assert that the linear operator\n"
00614       "and the preconditioner are indeed positive definite.");
00615     pl->set("Verbosity", verbosity_default_,
00616       "What type(s) of solver information should be outputted\n"
00617       "to the output stream.");
00618     pl->set("Output Style", outputStyle_default_,
00619       "What style is used for the solver information outputted\n"
00620       "to the output stream.");
00621     pl->set("Output Frequency", outputFreq_default_,
00622       "How often convergence information should be outputted\n"
00623       "to the output stream.");
00624     pl->set("Deflation Quorum", defQuorum_default_,
00625       "The number of linear systems that need to converge before\n"
00626       "they are deflated.  This number should be <= block size.");
00627     pl->set("Output Stream", outputStream_default_,
00628       "A reference-counted pointer to the output stream where all\n"
00629       "solver output is sent.");
00630     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00631       "When convergence information is printed, only show the maximum\n"
00632       "relative residual norm when the block size is greater than one.");
00633     pl->set("Implicit Residual Scaling", resScale_default_,
00634       "The type of scaling used in the residual convergence test.");
00635     // We leave the old name as a valid parameter for backwards
00636     // compatibility (so that validateParametersAndSetDefaults()
00637     // doesn't raise an exception if it encounters "Residual
00638     // Scaling").  The new name was added for compatibility with other
00639     // solvers, none of which use "Residual Scaling".
00640     pl->set("Residual Scaling", resScale_default_,
00641             "The type of scaling used in the residual convergence test.  This "
00642             "name is deprecated; the new name is \"Implicit Residual Scaling\".");
00643     pl->set("Timer Label", label_default_,
00644       "The string to use as a prefix for the timer labels.");
00645     //  defaultParams_->set("Restart Timers", restartTimers_);
00646     validParams_ = pl;
00647   }
00648   return validParams_;
00649 }
00650 
00651 
00652 // solve()
00653 template<class ScalarType, class MV, class OP>
00654 ReturnType PseudoBlockCGSolMgr<ScalarType,MV,OP>::solve() {
00655 
00656   // Set the current parameters if they were not set before.
00657   // NOTE:  This may occur if the user generated the solver manager with the default constructor and
00658   // then didn't set any parameters using setParameters().
00659   if (!isSet_) { setParameters( params_ ); }
00660 
00661   Teuchos::BLAS<int,ScalarType> blas;
00662 
00663   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockCGSolMgrLinearProblemFailure,
00664                      "Belos::PseudoBlockCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00665 
00666   // Create indices for the linear systems to be solved.
00667   int startPtr = 0;
00668   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00669   int numCurrRHS = numRHS2Solve;
00670 
00671   std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
00672   for (int i=0; i<numRHS2Solve; ++i) {
00673     currIdx[i] = startPtr+i;
00674     currIdx2[i]=i;
00675   }
00676 
00677   // Inform the linear problem of the current linear system to solve.
00678   problem_->setLSIndex( currIdx );
00679 
00681   // Parameter list
00682   Teuchos::ParameterList plist;
00683 
00684   plist.set("Assert Positive Definiteness",assertPositiveDefiniteness_);
00685 
00686   // Reset the status test.
00687   outputTest_->reset();
00688 
00689   // Assume convergence is achieved, then let any failed convergence set this to false.
00690   bool isConverged = true;
00691 
00693   // Pseudo-Block CG solver
00694 
00695   Teuchos::RCP<PseudoBlockCGIter<ScalarType,MV,OP> > block_cg_iter
00696     = Teuchos::rcp( new PseudoBlockCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
00697 
00698   // Enter solve() iterations
00699   {
00700 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00701     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00702 #endif
00703 
00704     while ( numRHS2Solve > 0 ) {
00705 
00706       // Reset the active / converged vectors from this block
00707       std::vector<int> convRHSIdx;
00708       std::vector<int> currRHSIdx( currIdx );
00709       currRHSIdx.resize(numCurrRHS);
00710 
00711       // Reset the number of iterations.
00712       block_cg_iter->resetNumIters();
00713 
00714       // Reset the number of calls that the status test output knows about.
00715       outputTest_->resetNumCalls();
00716 
00717       // Get the current residual for this block of linear systems.
00718       Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
00719 
00720       // Get a new state struct and initialize the solver.
00721       CGIterationState<ScalarType,MV> newState;
00722       newState.R = R_0;
00723       block_cg_iter->initializeCG(newState);
00724 
00725       while(1) {
00726 
00727         // tell block_gmres_iter to iterate
00728         try {
00729           block_cg_iter->iterate();
00730 
00732           //
00733           // check convergence first
00734           //
00736           if ( convTest_->getStatus() == Passed ) {
00737 
00738             // Figure out which linear systems converged.
00739             std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
00740 
00741             // If the number of converged linear systems is equal to the
00742             // number of current linear systems, then we are done with this block.
00743             if (convIdx.size() == currRHSIdx.size())
00744               break;  // break from while(1){block_cg_iter->iterate()}
00745 
00746             // Inform the linear problem that we are finished with this current linear system.
00747             problem_->setCurrLS();
00748 
00749             // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
00750             int have = 0;
00751             std::vector<int> unconvIdx(currRHSIdx.size());
00752             for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
00753               bool found = false;
00754               for (unsigned int j=0; j<convIdx.size(); ++j) {
00755                 if (currRHSIdx[i] == convIdx[j]) {
00756                   found = true;
00757                   break;
00758                 }
00759               }
00760               if (!found) {
00761                 currIdx2[have] = currIdx2[i];
00762                 currRHSIdx[have++] = currRHSIdx[i];
00763               }
00764             }
00765             currRHSIdx.resize(have);
00766             currIdx2.resize(have);
00767 
00768             // Set the remaining indices after deflation.
00769             problem_->setLSIndex( currRHSIdx );
00770 
00771             // Get the current residual vector.
00772             std::vector<MagnitudeType> norms;
00773             R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
00774             for (int i=0; i<have; ++i) { currIdx2[i] = i; }
00775 
00776             // Set the new state and initialize the solver.
00777             CGIterationState<ScalarType,MV> defstate;
00778             defstate.R = R_0;
00779             block_cg_iter->initializeCG(defstate);
00780           }
00781 
00783           //
00784           // check for maximum iterations
00785           //
00787           else if ( maxIterTest_->getStatus() == Passed ) {
00788             // we don't have convergence
00789             isConverged = false;
00790             break;  // break from while(1){block_cg_iter->iterate()}
00791           }
00792 
00794           //
00795           // we returned from iterate(), but none of our status tests Passed.
00796           // something is wrong, and it is probably our fault.
00797           //
00799 
00800           else {
00801             TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00802                                "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
00803           }
00804         }
00805         catch (const std::exception &e) {
00806           printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration "
00807                                    << block_cg_iter->getNumIters() << std::endl
00808                                    << e.what() << std::endl;
00809           throw;
00810         }
00811       }
00812 
00813       // Inform the linear problem that we are finished with this block linear system.
00814       problem_->setCurrLS();
00815 
00816       // Update indices for the linear systems to be solved.
00817       startPtr += numCurrRHS;
00818       numRHS2Solve -= numCurrRHS;
00819 
00820       if ( numRHS2Solve > 0 ) {
00821 
00822         numCurrRHS = numRHS2Solve;
00823         currIdx.resize( numCurrRHS );
00824         currIdx2.resize( numCurrRHS );
00825         for (int i=0; i<numCurrRHS; ++i)
00826           { currIdx[i] = startPtr+i; currIdx2[i] = i; }
00827 
00828         // Set the next indices.
00829         problem_->setLSIndex( currIdx );
00830       }
00831       else {
00832         currIdx.resize( numRHS2Solve );
00833       }
00834 
00835     }// while ( numRHS2Solve > 0 )
00836 
00837   }
00838 
00839   // print final summary
00840   sTest_->print( printer_->stream(FinalSummary) );
00841 
00842   // print timing information
00843 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00844   // Calling summarize() can be expensive, so don't call unless the
00845   // user wants to print out timing details.  summarize() will do all
00846   // the work even if it's passed a "black hole" output stream.
00847   if (verbosity_ & TimingDetails)
00848     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
00849 #endif
00850 
00851   // get iteration information for this solve
00852   numIters_ = maxIterTest_->getNumIters();
00853 
00854   if (!isConverged ) {
00855     return Unconverged; // return from PseudoBlockCGSolMgr::solve()
00856   }
00857   return Converged; // return from PseudoBlockCGSolMgr::solve()
00858 }
00859 
00860 //  This method requires the solver manager to return a std::string that describes itself.
00861 template<class ScalarType, class MV, class OP>
00862 std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::description() const
00863 {
00864   std::ostringstream oss;
00865   oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00866   oss << "{";
00867   oss << "}";
00868   return oss.str();
00869 }
00870 
00871 } // end Belos namespace
00872 
00873 #endif /* BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines