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