Belos Version of the Day
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     // Method to convert std::string to enumerated type for residual.
00236     Belos::ScaleType convertStringToScaleType( const std::string& scaleType ) {
00237       if (scaleType == "Norm of Initial Residual") {
00238         return Belos::NormOfInitRes;
00239       } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00240         // This algorithm only provides true residuals, not preconditioned residuals, 
00241         // so the residuals should only be scaled by the unpreconditioned initial residual.
00242         return Belos::NormOfInitRes;
00243       } else if (scaleType == "Norm of RHS") {
00244         return Belos::NormOfRHS;
00245       } else if (scaleType == "None") {
00246         return Belos::None;
00247       } else 
00248         TEST_FOR_EXCEPTION( true ,std::logic_error,
00249           "Belos::PseudoBlockCGSolMgr(): Invalid residual scaling type.");
00250     }
00251 
00252     // Linear problem.
00253     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00254     
00255     // Output manager.
00256     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00257     Teuchos::RCP<std::ostream> outputStream_;
00258 
00259     // Status test.
00260     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00261     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00262     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00263     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00264 
00265     // Orthogonalization manager.
00266     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00267 
00268      // Current parameter list.
00269     Teuchos::RCP<Teuchos::ParameterList> params_;
00270    
00271     // Default solver values.
00272     static const MagnitudeType convtol_default_;
00273     static const int maxIters_default_;
00274     static const bool assertPositiveDefiniteness_default_;
00275     static const bool showMaxResNormOnly_default_;
00276     static const int verbosity_default_;
00277     static const int outputStyle_default_;
00278     static const int outputFreq_default_;
00279     static const int defQuorum_default_;
00280     static const std::string resScale_default_; 
00281     static const std::string label_default_;
00282     static const Teuchos::RCP<std::ostream> outputStream_default_;
00283 
00284     // Current solver values.
00285     MagnitudeType convtol_;
00286     int maxIters_, numIters_;
00287     int verbosity_, outputStyle_, outputFreq_, defQuorum_;
00288     bool assertPositiveDefiniteness_, showMaxResNormOnly_;
00289     std::string resScale_;       
00290  
00291     // Timers.
00292     std::string label_;
00293     Teuchos::RCP<Teuchos::Time> timerSolve_;
00294 
00295     // Internal state variables.
00296     bool isSet_;
00297   };
00298   
00299   
00300 // Default solver values.
00301 template<class ScalarType, class MV, class OP>
00302 const typename PseudoBlockCGSolMgr<ScalarType,MV,OP>::MagnitudeType PseudoBlockCGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00303 
00304 template<class ScalarType, class MV, class OP>
00305 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00306 
00307 template<class ScalarType, class MV, class OP>
00308 const bool PseudoBlockCGSolMgr<ScalarType,MV,OP>::assertPositiveDefiniteness_default_ = true;
00309 
00310 template<class ScalarType, class MV, class OP>
00311 const bool PseudoBlockCGSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00312 
00313 template<class ScalarType, class MV, class OP>
00314 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00315 
00316 template<class ScalarType, class MV, class OP>
00317 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00318 
00319 template<class ScalarType, class MV, class OP>
00320 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00321 
00322 template<class ScalarType, class MV, class OP>
00323 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::defQuorum_default_ = 1;
00324 
00325 template<class ScalarType, class MV, class OP>
00326 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::resScale_default_ = "Norm of Initial Residual";
00327 
00328 template<class ScalarType, class MV, class OP>
00329 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00330 
00331 template<class ScalarType, class MV, class OP>
00332 const Teuchos::RCP<std::ostream> PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00333 
00334 
00335 // Empty Constructor
00336 template<class ScalarType, class MV, class OP>
00337 PseudoBlockCGSolMgr<ScalarType,MV,OP>::PseudoBlockCGSolMgr() :
00338   outputStream_(outputStream_default_),
00339   convtol_(convtol_default_),
00340   maxIters_(maxIters_default_),
00341   verbosity_(verbosity_default_),
00342   outputStyle_(outputStyle_default_),
00343   outputFreq_(outputFreq_default_),
00344   defQuorum_(defQuorum_default_),
00345   assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
00346   showMaxResNormOnly_(showMaxResNormOnly_default_),
00347   resScale_(resScale_default_),
00348   label_(label_default_),
00349   isSet_(false)
00350 {}
00351 
00352 // Basic Constructor
00353 template<class ScalarType, class MV, class OP>
00354 PseudoBlockCGSolMgr<ScalarType,MV,OP>::PseudoBlockCGSolMgr( 
00355                const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00356                const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00357   problem_(problem),
00358   outputStream_(outputStream_default_),
00359   convtol_(convtol_default_),
00360   maxIters_(maxIters_default_),
00361   verbosity_(verbosity_default_),
00362   outputStyle_(outputStyle_default_),
00363   outputFreq_(outputFreq_default_),
00364   defQuorum_(defQuorum_default_),
00365   assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
00366   showMaxResNormOnly_(showMaxResNormOnly_default_),
00367   resScale_(resScale_default_),
00368   label_(label_default_),
00369   isSet_(false)
00370 {
00371   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00372 
00373   // If the parameter list pointer is null, then set the current parameters to the default parameter list.
00374   if (!is_null(pl)) {
00375     // Set the parameters using the list that was passed in.
00376     setParameters( pl );  
00377   }
00378 }
00379 
00380 template<class ScalarType, class MV, class OP>
00381 void PseudoBlockCGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00382 {
00383   using Teuchos::ParameterList;
00384   using Teuchos::parameterList;
00385   using Teuchos::RCP;
00386 
00387   RCP<const ParameterList> defaultParams = getValidParameters();
00388 
00389   // Create the internal parameter list if one doesn't already exist.
00390   if (params_.is_null()) {
00391     params_ = parameterList (*defaultParams);
00392   } else {
00393     params->validateParameters (*defaultParams);
00394   }
00395 
00396   // Check for maximum number of iterations
00397   if (params->isParameter("Maximum Iterations")) {
00398     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00399 
00400     // Update parameter in our list and in status test.
00401     params_->set("Maximum Iterations", maxIters_);
00402     if (maxIterTest_!=Teuchos::null)
00403       maxIterTest_->setMaxIters( maxIters_ );
00404   }
00405 
00406   // Check if positive definiteness assertions are to be performed
00407   if (params->isParameter("Assert Positive Definiteness")) {
00408     assertPositiveDefiniteness_ = params->get("Assert Positive Definiteness",assertPositiveDefiniteness_default_);
00409 
00410     // Update parameter in our list.
00411     params_->set("Assert Positive Definiteness", assertPositiveDefiniteness_);
00412   }
00413 
00414   // Check to see if the timer label changed.
00415   if (params->isParameter("Timer Label")) {
00416     std::string tempLabel = params->get("Timer Label", label_default_);
00417 
00418     // Update parameter in our list and solver timer
00419     if (tempLabel != label_) {
00420       label_ = tempLabel;
00421       params_->set("Timer Label", label_);
00422       std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00423 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00424       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00425 #endif
00426     }
00427   }
00428 
00429   // Check for a change in verbosity level
00430   if (params->isParameter("Verbosity")) {
00431     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00432       verbosity_ = params->get("Verbosity", verbosity_default_);
00433     } else {
00434       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00435     }
00436 
00437     // Update parameter in our list.
00438     params_->set("Verbosity", verbosity_);
00439     if (printer_ != Teuchos::null)
00440       printer_->setVerbosity(verbosity_);
00441   }
00442 
00443   // Check for a change in output style
00444   if (params->isParameter("Output Style")) {
00445     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00446       outputStyle_ = params->get("Output Style", outputStyle_default_);
00447     } else {
00448       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00449     }
00450 
00451     // Reconstruct the convergence test if the explicit residual test is not being used.
00452     params_->set("Output Style", outputStyle_);
00453     outputTest_ = Teuchos::null;
00454   }
00455 
00456   // output stream
00457   if (params->isParameter("Output Stream")) {
00458     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00459 
00460     // Update parameter in our list.
00461     params_->set("Output Stream", outputStream_);
00462     if (printer_ != Teuchos::null)
00463       printer_->setOStream( outputStream_ );
00464   }
00465 
00466   // frequency level
00467   if (verbosity_ & Belos::StatusTestDetails) {
00468     if (params->isParameter("Output Frequency")) {
00469       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00470     }
00471 
00472     // Update parameter in out list and output status test.
00473     params_->set("Output Frequency", outputFreq_);
00474     if (outputTest_ != Teuchos::null)
00475       outputTest_->setOutputFrequency( outputFreq_ );
00476   }
00477 
00478   // Create output manager if we need to.
00479   if (printer_ == Teuchos::null) {
00480     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00481   }
00482   
00483   // Convergence
00484   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00485   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00486 
00487   // Check for convergence tolerance
00488   if (params->isParameter("Convergence Tolerance")) {
00489     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00490 
00491     // Update parameter in our list and residual tests.
00492     params_->set("Convergence Tolerance", convtol_);
00493     if (convTest_ != Teuchos::null)
00494       convTest_->setTolerance( convtol_ );
00495   }
00496 
00497   if (params->isParameter("Show Maximum Residual Norm Only")) {
00498     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00499 
00500     // Update parameter in our list and residual tests
00501     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00502     if (convTest_ != Teuchos::null)
00503       convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00504   }
00505 
00506   // Check for a change in scaling, if so we need to build new residual tests.
00507   bool newResTest = false;
00508   {
00509     // "Residual Scaling" is the old parameter name; "Implicit
00510     // Residual Scaling" is the new name.  We support both options for
00511     // backwards compatibility.
00512     std::string tempResScale = resScale_;
00513     bool implicitResidualScalingName = false;
00514     if (params->isParameter ("Residual Scaling")) {
00515       tempResScale = params->get<std::string> ("Residual Scaling");
00516     } 
00517     else if (params->isParameter ("Implicit Residual Scaling")) {
00518       tempResScale = params->get<std::string> ("Implicit Residual Scaling");
00519       implicitResidualScalingName = true;
00520     }
00521 
00522     // Only update the scaling if it's different.
00523     if (resScale_ != tempResScale) {
00524       Belos::ScaleType resScaleType = convertStringToScaleType( tempResScale );
00525       resScale_ = tempResScale;
00526 
00527       // Update parameter in our list and residual tests, using the
00528       // given parameter name.
00529       if (implicitResidualScalingName) {
00530   params_->set ("Implicit Residual Scaling", resScale_);
00531       } 
00532       else {
00533   params_->set ("Residual Scaling", resScale_); 
00534       }
00535 
00536       if (! convTest_.is_null()) {
00537         try {
00538           convTest_->defineScaleForm( resScaleType, Belos::TwoNorm );
00539         }
00540         catch (std::exception& e) {
00541     // Make sure the convergence test gets constructed again.
00542     newResTest = true;
00543         }
00544       }
00545     }      
00546   }
00547 
00548   // Get the deflation quorum, or number of converged systems before deflation is allowed
00549   if (params->isParameter("Deflation Quorum")) {
00550     defQuorum_ = params->get("Deflation Quorum", defQuorum_);
00551     params_->set("Deflation Quorum", defQuorum_);
00552     if (convTest_ != Teuchos::null)
00553       convTest_->setQuorum( defQuorum_ );
00554   }
00555 
00556   // Create status tests if we need to.
00557 
00558   // Basic test checks maximum iterations and native residual.
00559   if (maxIterTest_ == Teuchos::null)
00560     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00561   
00562   // Implicit residual test, using the native residual to determine if convergence was achieved.
00563   if (convTest_ == Teuchos::null || newResTest) {
00564     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, defQuorum_, showMaxResNormOnly_ ) );
00565     convTest_->defineScaleForm( convertStringToScaleType( resScale_ ), Belos::TwoNorm );
00566   } 
00567 
00568   if (sTest_ == Teuchos::null || newResTest)
00569     sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00570   
00571   if (outputTest_ == Teuchos::null || newResTest) {
00572 
00573     // Create the status test output class.
00574     // This class manages and formats the output from the status test.
00575     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00576     outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00577 
00578     // Set the solver string for the output test
00579     std::string solverDesc = " Pseudo Block CG ";
00580     outputTest_->setSolverDesc( solverDesc );
00581 
00582   }
00583 
00584   // Create the timer if we need to.
00585   if (timerSolve_ == Teuchos::null) {
00586     std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00587 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00588     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00589 #endif
00590   }
00591 
00592   // Inform the solver manager that the current parameters were set.
00593   isSet_ = true;
00594 }
00595 
00596     
00597 template<class ScalarType, class MV, class OP>
00598 Teuchos::RCP<const Teuchos::ParameterList>
00599 PseudoBlockCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00600 {
00601   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00602   if (is_null(validPL)) {
00603     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00604     // Set all the valid parameters and their default values.
00605     pl= Teuchos::rcp( new Teuchos::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     validPL = pl;
00647   }
00648   return validPL;
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   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       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