Belos Package Browser (Single Doxygen Collection) Development
BelosGmresSolMgr.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_GmresSolMgr_hpp
00043 #define __Belos_GmresSolMgr_hpp
00044 
00048 
00049 #include <BelosGmresBaseIteration.hpp>
00050 
00051 #include <BelosOutputManager.hpp>
00052 #include <BelosSolverManager.hpp>
00053 #include <BelosStatusTestGenResNorm.hpp>
00054 #include <BelosStatusTestOutputFactory.hpp>
00055 #include <BelosOrthoManagerFactory.hpp>
00056 
00057 #include <functional> // std::logical_and
00058 #include <numeric> // std::accumulate
00059 
00060 namespace Belos {
00061 
00072   template<class Scalar, class MV, class OP>
00073   class GmresSolMgr : public SolverManager<Scalar, MV, OP> {
00074   public:
00075     typedef Scalar scalar_type;
00076     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00077     typedef MV multivector_type;
00078     typedef OP operator_type;
00079 
00080   private:
00081     typedef MultiVecTraits<Scalar,MV> MVT;
00082     typedef OperatorTraits<Scalar,MV,OP> OPT;
00083     typedef Teuchos::ScalarTraits<Scalar> STS;
00084     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00085     typedef GmresBaseIteration<Scalar, MV, OP> iteration_type;
00086 
00087   public:
00089 
00090 
00096     GmresSolMgr (const Teuchos::RCP<LinearProblem<Scalar,MV,OP> >& problem,
00097      const Teuchos::RCP<Teuchos::ParameterList>& params,
00098      const bool debug = false) : 
00099       problem_ (validatedProblem (problem)), 
00100       debug_ (debug)
00101     { 
00102       setParametersImpl (params);
00103     }
00104 
00106     GmresSolMgr() 
00107     {
00108       throw std::logic_error("TODO: implement deferring setting parameters until solve()");
00109     }
00110 
00112     virtual ~GmresSolMgr() {}
00113 
00115 
00117 
00118 
00120     const LinearProblem<Scalar,MV,OP>& getProblem() const {
00121       return *problem_;
00122     }
00123 
00132     static Teuchos::RCP<const Teuchos::ParameterList> getDefaultParameters();
00133 
00139     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
00140       return getDefaultParameters();
00141     }
00142 
00144     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const {
00145       return params_;
00146     }
00147 
00154     int getNumIters() const { 
00155       if (iter_.is_null())
00156   throw std::logic_error("The number of iterations is undefined, since "
00157              "the iteration has not yet been initialized.");
00158       return iter_->getNumIters();
00159     }
00160 
00168     bool isLOADetected() const {
00169       return false;
00170     }
00171 
00179     std::vector<std::pair<int, int> >
00180     totalNumIters () const {
00181       return totalNumIters_;
00182     }
00183 
00185 
00187 
00188 
00192     void 
00193     setProblem (const Teuchos::RCP<LinearProblem<Scalar,MV,OP> >& problem)
00194     {
00195       problem_ = validatedProblem (problem);
00196       // Changing the problem requires (re)initializing possibly a lot
00197       // of different things, including the status tests,
00198       // orthogonalization method, and the GmresBaseIteration object.
00199       reset (Belos::Problem);
00200     }
00201 
00213     void 
00214     setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params);
00215 
00226     void 
00227     setUserConvStatusTest (const Teuchos::RCP<StatusTest<Scalar,MV,OP> >& userConvTest);
00229 
00231 
00232 
00248     void reset (const ResetType type);
00250 
00252 
00253 
00264     ReturnType solve();
00265 
00284     void 
00285     setRHS (const Teuchos::RCP<const MV>& B) 
00286     {
00287       TEUCHOS_TEST_FOR_EXCEPTION(B.is_null(), std::invalid_argument,
00288        "Belos::GmresSolMgr::setRHS(): the given new right-"
00289        "hand side B is null.");
00290       // This "unsets" the problem: it makes the initial residual (and
00291       // the left-preconditioned initial residual, if applicable)
00292       // invalid.  We don't have to call setProblem() to reset the
00293       // problem, since setProblem() will be called in solve().
00294       problem_->setRHS (B);
00295     }
00296 
00312     void 
00313     setLHS (const Teuchos::RCP<MV>& X) 
00314     {
00315       TEUCHOS_TEST_FOR_EXCEPTION(X.is_null(), std::invalid_argument,
00316        "Belos::GmresSolMgr::setLHS(): the given new initial "
00317        "guess X is null.");
00318       // This "unsets" the problem; it makes the initial residual (and
00319       // the left-preconditioned initial residual, if applicable)
00320       // invalid.  We don't have to call setProblem() to reset the
00321       // problem, since setProblem() will be called in solve().
00322       problem_->setLHS (X);
00323     }
00324 
00333     void
00334     changeStoppingCriteria (const magnitude_type convTol,
00335           const int maxItersPerRestart,
00336           const int maxNumRestarts)
00337     {
00338       // Testing maxNumRestarts now, rather than below, ensures the
00339       // strong exception guarantee: either no state will be changed,
00340       // or no exceptions will be thrown.
00341       TEUCHOS_TEST_FOR_EXCEPTION(maxNumRestarts < 0, std::invalid_argument,
00342        "The maximum number of restart cycle(s) must be "
00343        "nonnegative, but a value of " << maxNumRestarts 
00344        << " was specified.");
00345       rebuildStatusTests (convTol, maxItersPerRestart);
00346       params_->set ("Maximum Restarts", maxNumRestarts,
00347         "Maximum number of restart cycle(s) allowed for each "
00348         "right-hand side solved.");
00349     }
00351     
00352   private:
00354     Teuchos::RCP<LinearProblem<Scalar, MV, OP> > problem_;
00355 
00357     Teuchos::RCP<OutputManager<Scalar> > outMan_;
00358 
00360     Teuchos::RCP<Teuchos::ParameterList> params_;
00361 
00363     Teuchos::RCP<StatusTest<Scalar,MV,OP> > convTest_;
00364 
00366     Teuchos::RCP<StatusTest<Scalar,MV,OP> > statusTest_;
00367 
00374     Teuchos::RCP<StatusTest<Scalar,MV,OP> > userConvTest_;
00375 
00377     Teuchos::RCP<StatusTestOutput<Scalar,MV,OP> > outTest_;
00378 
00380     Teuchos::RCP<const OrthoManager<Scalar, MV> > orthoMan_;
00381 
00383     Teuchos::RCP<iteration_type> iter_;
00384 
00390     std::vector<std::pair<int, int> > totalNumIters_;
00391 
00393     bool debug_;
00394 
00404     void 
00405     setParametersImpl (const Teuchos::RCP<Teuchos::ParameterList>& params);
00406 
00413     void 
00414     rebuildStatusTests (Teuchos::RCP<Teuchos::ParameterList> plist = 
00415       Teuchos::null);
00416 
00429     void
00430     rebuildStatusTests (const magnitude_type convTol,
00431       const int maxItersPerRestart);
00432 
00439     void
00440     rebuildOrthoManager (Teuchos::RCP<Teuchos::ParameterList> plist = 
00441        Teuchos::null);
00442 
00449     void 
00450     rebuildIteration (Teuchos::RCP<Teuchos::ParameterList> plist = 
00451           Teuchos::null);
00452 
00464     static Teuchos::RCP<LinearProblem<Scalar, MV, OP> >
00465     validatedProblem (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem)
00466     {
00467       const char prefix[] = "Belos::GmresSolMgr: ";
00468       TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
00469        prefix << "The linear problem (Belos::LinearProblem "
00470        "instance) that you want me to solve is null.");
00471       if (! problem->isProblemSet())
00472   {
00473     problem->setProblem();
00474     const bool notInitialized = problem->getOperator().is_null() || 
00475       problem->getRHS().is_null() || problem->getLHS().is_null();
00476     if (notInitialized)
00477       {
00478         std::ostringstream os;
00479         os << prefix << "The given linear problem (Belos::Linear"
00480     "Problem) instance is not fully initialized: "
00481      << (problem->getOperator().is_null() ? "the operator A is null; " : "")
00482      << (problem->getRHS().is_null() ? "the right-hand side B is null; " : "")
00483      << (problem->getLHS().is_null() ? "the initial guess X is null" : "")
00484      << ".";
00485         TEUCHOS_TEST_FOR_EXCEPTION(notInitialized, std::invalid_argument, os.str());
00486       }
00487     TEUCHOS_TEST_FOR_EXCEPTION(! problem->isProblemSet(), std::logic_error,
00488            prefix << "Although the given LinearProblem instance "
00489            "has non-null operator (matrix A), right-hand side B,"
00490            " and initial guess X, and although its setProblem() "
00491            "method has been called, its isProblemSet() method "
00492            "returns false.  This should not happen.");
00493   }
00494       return problem;
00495     }
00496 
00498     static ScaleType
00499     stringToScaleType (const std::string& scaleType) 
00500     {
00501       if (scaleType == "Norm of Initial Residual") 
00502         return Belos::NormOfInitRes;
00503       else if (scaleType == "Norm of Preconditioned Initial Residual") 
00504         return Belos::NormOfPrecInitRes;
00505       else if (scaleType == "Norm of RHS" || scaleType == "Norm of Right-Hand Side")
00506         return Belos::NormOfRHS;
00507       else if (scaleType == "None")
00508         return Belos::None;
00509       else {
00510         TEUCHOS_TEST_FOR_EXCEPTION (true, std::logic_error,
00511           "Invalid residual scaling type \"" << scaleType 
00512           << "\".");
00513       }
00514     }
00515 
00527     static Teuchos::RCP<StatusTestOutput<Scalar,MV,OP> > 
00528     initOutputTest (const Teuchos::RCP<OutputManager<Scalar> >& outMan,
00529         const Teuchos::RCP<StatusTest<Scalar,MV,OP> >& statusTest,
00530         const OutputType outStyle,
00531         const int outFreq,
00532         const std::string& solverDesc = "",  // e.g., " CA-GMRES "
00533         const std::string& precondDesc = ""); // Preconditioner description
00534 
00540     //
00555     static Teuchos::RCP<StatusTest<Scalar,MV,OP> >
00556     initConvTest (const magnitude_type convTol,
00557       const bool haveLeftPreconditioner,
00558       const std::string& implicitScaleType,
00559       const std::string& explicitScaleType,
00560       Teuchos::RCP<StatusTest<Scalar,MV,OP> > userConvTest);
00561 
00574     static Teuchos::RCP<StatusTest<Scalar,MV,OP> >
00575     initStatusTest (Teuchos::RCP<StatusTest<Scalar,MV,OP> > convTest,
00576         const int maxIters);
00577 
00596     static Teuchos::RCP<OutputManager<Scalar> >
00597     initOutputManager (const Teuchos::RCP<OutputManager<Scalar> >& outMan,
00598            const MsgType verbosity,
00599            const Teuchos::RCP<std::ostream>& outStream);
00600 
00601   };
00602 
00603 
00604   template<class Scalar, class MV, class OP>
00605   Teuchos::RCP<OutputManager<Scalar> >
00606   GmresSolMgr<Scalar,MV,OP>::
00607   initOutputManager (const Teuchos::RCP<OutputManager<Scalar> >& outMan,
00608          const MsgType verbosity,
00609          const Teuchos::RCP<std::ostream>& outStream)
00610   {
00611     // If outStream is null, use std::cout as the default.
00612     Teuchos::RCP<std::ostream> out = 
00613       outStream.is_null() ? Teuchos::rcpFromRef(std::cout) : outStream;
00614     if (outMan.is_null())
00615       return Teuchos::rcp (new OutputManager<Scalar> (verbosity, out));
00616     else
00617       {
00618   outMan->setVerbosity (verbosity);
00619   outMan->setOStream (out);
00620   return outMan;
00621       }
00622   }
00623 
00624   template<class Scalar, class MV, class OP>
00625   Teuchos::RCP<StatusTestOutput<Scalar,MV,OP> > 
00626   GmresSolMgr<Scalar,MV,OP>::
00627   initOutputTest (const Teuchos::RCP<OutputManager<Scalar> >& outMan,
00628       const Teuchos::RCP<StatusTest<Scalar,MV,OP> >& statusTest,
00629       const OutputType outStyle,
00630       const int outFreq,
00631       const std::string& solverDesc,
00632       const std::string& precondDesc)
00633   {
00634     using Teuchos::RCP;
00635 
00636     TEUCHOS_TEST_FOR_EXCEPTION(outMan.is_null(), std::logic_error,
00637       "Construction / reinitialization of the output status "
00638       "test depends on the OutputManager being initialized, "
00639       "but it has not yet been initialized.");
00640     TEUCHOS_TEST_FOR_EXCEPTION(statusTest.is_null(), std::logic_error,
00641            "Construction / reinitialization of the output status "
00642            "test depends on the status test being initialized, "
00643            "but it has not yet been initialized.");
00644     StatusTestOutputFactory<Scalar,MV,OP> factory (outStyle);
00645 
00646     RCP<StatusTestOutput<Scalar, MV, OP> > newOutTest = 
00647       factory.create (outMan, statusTest, outFreq, 
00648           Passed+Failed+Undefined);
00649 
00650     // newOutTest->setOutputManager (outMan);
00651     // newOutTest->setOutputFrequency (outFreq);
00652     // newOutTest->setChild (statusTest);
00653 
00654     // The factory doesn't know how to set the description strings,
00655     // so we do so here.
00656     if (solverDesc != "")
00657       newOutTest->setSolverDesc (solverDesc);
00658     if (precondDesc != "")
00659       newOutTest->setPrecondDesc (precondDesc);
00660 
00661     return newOutTest;
00662   }
00663 
00664 
00665   template<class Scalar, class MV, class OP>
00666   Teuchos::RCP<StatusTest<Scalar,MV,OP> >
00667   GmresSolMgr<Scalar,MV,OP>::
00668   initConvTest (const magnitude_type convTol,
00669     const bool haveLeftPreconditioner,
00670     const std::string& implicitScaleType,
00671     const std::string& explicitScaleType,
00672     Teuchos::RCP<StatusTest<Scalar,MV,OP> > userConvTest)
00673   {
00674     using Teuchos::RCP;
00675     using Teuchos::rcp;
00676     typedef StatusTest<Scalar,MV,OP> base_test;
00677     typedef StatusTestGenResNorm<Scalar,MV,OP> res_norm_test;
00678     typedef StatusTestCombo<Scalar,MV,OP> combo_test;
00679 
00680     // "Deflation Quorum": number of converged systems before
00681     // deflation is allowed.  Cannot be larger than "Block Size",
00682     // which in our case is 1.  -1 is the default in
00683     // StatusTestGenResNorm.
00684     const int defQuorum = -1;
00685 
00686     // "Show Maximum Residual Norm Only" is only meaningful when the
00687     // "Block Size" parameter is > 1.  Our solver only supports a
00688     // Block Size of 1.
00689     const bool showMaxResNormOnly = false;
00690 
00691     // The "implicit" residual test checks the "native" residual
00692     // norm to determine if convergence was achieved.  It is less
00693     // expensive than the "explicit" residual test.
00694     RCP<res_norm_test> implicitTest (new res_norm_test (convTol, defQuorum));
00695     implicitTest->defineScaleForm (stringToScaleType (implicitScaleType), 
00696            Belos::TwoNorm);
00697     implicitTest->setShowMaxResNormOnly (showMaxResNormOnly);
00698 
00699     // If there's a left preconditioner, create a combined status
00700     // test that check first the "explicit," then the "implicit"
00701     // residual norm, requiring that both have converged to within
00702     // the specified tolerance.  Otherwise, we only perform the
00703     // "implicit" test.
00704     RCP<res_norm_test> explicitTest;
00705     if (haveLeftPreconditioner) // ! problem_->getLeftPrec().is_null()
00706       {
00707   explicitTest = rcp (new res_norm_test (convTol, defQuorum));
00708   explicitTest->defineResForm (res_norm_test::Explicit, Belos::TwoNorm);
00709   explicitTest->defineScaleForm (stringToScaleType (explicitScaleType),
00710                Belos::TwoNorm);
00711   explicitTest->setShowMaxResNormOnly (showMaxResNormOnly);
00712       }
00713     // The "final" convergence test: 
00714     //
00715     // First, the implicit residual norm test,
00716     // Followed by the explicit residual norm test if applicable,
00717     // Followed by the user-defined convergence test if supplied.
00718     RCP<base_test> convTest;
00719     if (explicitTest.is_null())
00720       convTest = implicitTest;
00721     else
00722       // The "explicit" residual test is only performed once the
00723       // native ("implicit") residual is below the convergence
00724       // tolerance.
00725       convTest = rcp (new combo_test (combo_test::SEQ, 
00726               implicitTest, 
00727               explicitTest));
00728     if (! userConvTest.is_null())
00729       convTest = rcp (new combo_test (combo_test::SEQ, 
00730               convTest, 
00731               userConvTest));
00732     return convTest;
00733   }
00734 
00735 
00736   template<class Scalar, class MV, class OP>
00737   Teuchos::RCP<StatusTest<Scalar,MV,OP> >
00738   GmresSolMgr<Scalar,MV,OP>::
00739   initStatusTest (Teuchos::RCP<StatusTest<Scalar,MV,OP> > convTest,
00740       const int maxIters)
00741   {
00742     using Teuchos::RCP;
00743     using Teuchos::rcp;
00744     typedef StatusTest<Scalar,MV,OP> base_test;
00745     typedef StatusTestMaxIters<Scalar,MV,OP> max_iter_test;
00746     typedef StatusTestCombo<Scalar,MV,OP> combo_test;
00747 
00748     // Stopping criterion for maximum number of iterations.
00749     //
00750     // We could just keep this around and call setMaxIters() on it
00751     // when the maximum number of iterations changes, but why
00752     // bother?  It's not a heavyweight object.  Just make a new one.
00753     RCP<max_iter_test> maxIterTest (new max_iter_test (maxIters));
00754 
00755     // The "final" stopping criterion:
00756     //
00757     // Either we've run out of iterations, OR we've converged.
00758     return rcp (new combo_test (combo_test::OR, maxIterTest, convTest));
00759   }
00760 
00761   template<class Scalar, class MV, class OP>
00762   Teuchos::RCP<const Teuchos::ParameterList>
00763   GmresSolMgr<Scalar,MV,OP>::getDefaultParameters()
00764   {
00765     using Teuchos::RCP;
00766     using Teuchos::ParameterList;
00767     using Teuchos::parameterList;
00768 
00769     // FIXME (mfh 15 Feb 2011) This makes the routine nonreentrant.
00770     static RCP<const ParameterList> defaultParams;
00771     if (defaultParams.is_null())
00772       {
00773   RCP<ParameterList> params = parameterList();
00774   params->setName ("GmresSolMgr");
00775   //
00776   // The OutputManager only depends on the verbosity and the
00777   // output stream.
00778   //
00779   // Only display errors by default.
00780   //
00781   // FIXME (mfh 16 Feb 2011) Annoyingly, I have to write this as
00782   // an int and not a MsgType, otherwise it's saved as a string.
00783   MsgType defaultVerbosity = Belos::Errors;
00784   params->set ("Verbosity", static_cast<int>(defaultVerbosity),
00785          "What type(s) of solver information should be written "
00786          "to the output stream.");
00787   // Output goes to std::cout by default.
00788   RCP<std::ostream> defaultOutStream = Teuchos::rcpFromRef (std::cout);
00789   params->set ("Output Stream", defaultOutStream, 
00790          "A reference-counted pointer to the output stream where "
00791          "all solver output is sent.");
00792   //
00793   // The status test output class (created via
00794   // StatusTestOutputFactory) depends on the output style and
00795   // frequency.
00796   //
00797   // FIXME (mfh 16 Feb 2011) Annoyingly, even when I insist (via
00798   // a cast to OutputType) that the value has type OutputType,
00799   // it's stored in the XML file as a Teuchos::any, and read
00800   // back in as a string.  Thus, I cast to int instead.
00801   //
00802   // Default output style.
00803   OutputType defaultOutputStyle = Belos::General;
00804   params->set ("Output Style", (int) defaultOutputStyle,
00805          "What style is used for the solver information written "
00806          "to the output stream.");
00807   // Output frequency level refers to the number of iterations
00808   // between status outputs.  The default is -1, which means
00809   // "never."
00810   const int defaultOutFreq = -1;
00811   params->set ("Output Frequency", defaultOutFreq,
00812          "How often (in terms of number of iterations) "
00813          "convergence information should be written to "
00814          "the output stream.  (-1 means \"never.\")");
00815 
00816   // Default orthogonalization type is "Simple," which does
00817   // block Gram-Schmidt for project() and MGS for normalize().
00818   const std::string defaultOrthoType ("Simple");
00819   params->set ("Orthogonalization", defaultOrthoType);
00820 
00821   // Parameters for the orthogonalization.  These depend on
00822   // the orthogonalization type, so if you change that, you
00823   // should also change its parameters.
00824   OrthoManagerFactory<Scalar, MV, OP> orthoFactory;
00825   RCP<const ParameterList> defaultOrthoParams = 
00826     orthoFactory.getDefaultParameters (defaultOrthoType);
00827   // Storing the ParameterList as a sublist, rather than as an
00828   // RCP, ensures correct input and output.
00829   params->set ("Orthogonalization Parameters", *defaultOrthoParams);
00830 
00831   // Maximum number of iterations per restart cycle.
00832   const int defaultMaxIters = 100;
00833   params->set ("Maximum Iterations", defaultMaxIters,
00834          "Maximum number of iterations allowed per restart cycle, "
00835          "for each right-hand side solved.");
00836   // Maximum number of restart cycle(s) per right-hand side.
00837   const int defaultMaxRestarts = 10;
00838   params->set ("Maximum Restarts", defaultMaxRestarts,
00839          "Maximum number of restart cycle(s) allowed for each "
00840          "right-hand side solved.");
00841 
00842   // Default convergence tolerance is the square root of
00843   // machine precision.
00844   const magnitude_type defaultConvTol = STM::squareroot (STM::eps());
00845   params->set ("Convergence Tolerance", defaultConvTol,
00846          "Relative residual tolerance that must be attained by "
00847          "the iterative solver in order for the linear system to "
00848          "be declared Converged.");
00849   //
00850   // The implicit (and explicit, if applicable) relative
00851   // residual tests needs an initial scaling, which makes the
00852   // norms "relative."  The default for the implicit test is
00853   // to use the (left-)preconditioned initial residual (which
00854   // is just the initial residual in the case of no left
00855   // preconditioning).  The default for the explicit test is
00856   // to use the non-preconditioned initial residual.
00857   //
00858   const std::string defaultImplicitScaleType = 
00859     "Norm of Preconditioned Initial Residual";
00860   params->set ("Implicit Residual Scaling", 
00861          defaultImplicitScaleType,      
00862          "The type of scaling used in the implicit residual "
00863          "convergence test.");
00864   const std::string defaultExplicitScaleType =
00865     "Norm of Initial Residual";
00866   params->set ("Explicit Residual Scaling", 
00867          defaultExplicitScaleType,
00868          "The type of scaling used in the explicit residual "
00869          "convergence test.");
00870 
00871   // TODO (mfh 21 Feb 2011) Fill in default values for the GMRES
00872   // iteration parameters.  Right now, the parameter list exists
00873   // but is empty.
00874   ParameterList iterParams;
00875   params->set ("Iteration Parameters", iterParams);
00876 
00877   // TODO (mfh 15 Feb 2011) Set any the other default parameters.
00878 
00879   defaultParams = params;
00880       }
00881     return defaultParams;
00882   }
00883 
00884   template<class Scalar, class MV, class OP>
00885   void
00886   GmresSolMgr<Scalar,MV,OP>::
00887   setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
00888   {
00889     using Teuchos::ParameterList;
00890     using Teuchos::rcp_const_cast;
00891     setParametersImpl (params);
00892   }
00893 
00894   template<class Scalar, class MV, class OP>
00895   void
00896   GmresSolMgr<Scalar,MV,OP>::
00897   setParametersImpl (const Teuchos::RCP<Teuchos::ParameterList>& params)
00898   {
00899     using Teuchos::Exceptions::InvalidParameter;
00900     using Teuchos::Exceptions::InvalidParameterType;
00901     using Teuchos::null;
00902     using Teuchos::ParameterList;
00903     using Teuchos::parameterList;
00904     using Teuchos::RCP;
00905     using Teuchos::rcp;
00906     using std::endl;
00907 
00908     RCP<const ParameterList> defaultParams = getDefaultParameters();
00909     TEUCHOS_TEST_FOR_EXCEPTION(defaultParams.is_null(), std::logic_error, 
00910            "Belos::GmresSolMgr::setParametersImpl: The default "
00911            "parameter list is null; this should never happen.  "
00912            "Please report this bug to the Belos developers.");
00913     RCP<ParameterList> actualParams = params.is_null() ? 
00914       parameterList (*defaultParams) : params;
00915 
00916     // Fill in default values for parameters that aren't provided,
00917     // and make sure that all the provided parameters' values are
00918     // correct.
00919     //
00920     // FIXME (mfh 16 Feb 2011) Reading the output stream (which has
00921     // type RCP<std::ostream>) from a ParameterList may be impossible
00922     // if the ParameterList was read in from a file.  We hackishly
00923     // test for this by catching InvalidParameterType, setting the
00924     // output stream in actualParams to its default value, and redoing
00925     // the validation.  This is a hack because we don't know whether
00926     // the "Output Stream" parameter really caused that exception to
00927     // be thrown.
00928     bool success = false;
00929     try {
00930       actualParams->validateParametersAndSetDefaults (*defaultParams);
00931       success = true;
00932     } catch (InvalidParameterType&) {
00933       success = false;
00934     }
00935     if (! success) {
00936       RCP<std::ostream> outStream = 
00937   defaultParams->get<RCP<std::ostream> > ("Output Stream");
00938       actualParams->set ("Output Stream", outStream, 
00939        "A reference-counted pointer to the output "
00940        "stream where all solver output is sent.");
00941       // Retry the validation.
00942       actualParams->validateParametersAndSetDefaults (*defaultParams);
00943       success = true;
00944     }
00945 
00946     // If we don't have a problem to solve yet, we haven't previously
00947     // gotten parameters, or the maximum number of iterations per
00948     // restart cycle has changed, we will need to reset the solver.
00949     //
00950     // Changing the maximum number of iterations requires resetting
00951     // the solver, since GMRES usually preallocates this number of
00952     // vectors.
00953     //
00954     // FIXME (mfh 21 Feb 2011) Perhaps we should ask the specific
00955     // GMRES implementation whether it needs to be reset?
00956     const bool needToResetSolver = problem_.is_null() || params_.is_null() || 
00957       params_->get<int> ("Maximum Iterations") != 
00958       actualParams->get<int> ("Maximum Iterations");
00959 
00960     // Initialize the OutputManager.    
00961     {
00962       // FIXME (mfh 16 Feb 2011) Annoyingly, I have to read this back
00963       // in as an int rather than a Belos::MsgType.
00964       //
00965       // TODO (mfh 15 Feb 2011) Validate verbosity; MsgType is really
00966       // just an int, so it might have an invalid value.  (C++
00967       // typedefs don't have a nice type theory.)
00968       const MsgType verbosity = (MsgType) actualParams->get<int> ("Verbosity");
00969 
00970       // Reading the output stream from a ParameterList can be tricky
00971       // if the ParameterList was read in from a file.  Serialization
00972       // of pointers to arbitrary data is very, very difficult, and
00973       // ParameterList doesn't try.  If the XML file shows the type as
00974       // "any", the get() call below may throw InvalidParameterType.
00975       // In that case, we set the outStream to point to std::cout, as
00976       // a reasonable default.
00977       RCP<std::ostream> outStream;
00978       try {
00979   outStream = actualParams->get<RCP<std::ostream> > ("Output Stream");
00980       } catch (InvalidParameterType&) {
00981   outStream = Teuchos::rcpFromRef (std::cout);
00982       }
00983   
00984       // Sanity check
00985       //
00986       // FIXME (mfh 16 Feb 2011) Should this be a "black hole" stream
00987       // on procs other than Proc 0?
00988       if (outStream.is_null())
00989   outStream = Teuchos::rcpFromRef (std::cout);
00990 
00991       // This will do the right thing whether or not outMan_ has
00992       // previously been initialized (== is not null).
00993       outMan_ = initOutputManager (outMan_, verbosity, outStream);
00994     }
00995     TEUCHOS_TEST_FOR_EXCEPTION(outMan_.is_null(), std::logic_error,
00996            "Belos::GmresSolMgr::setParametersImpl: OutputManager "
00997            "instance is null after its initialization; this should"
00998            " never happen.  Please report this bug to the Belos "
00999            "developers.");
01000 
01001     // Now that we've initialized the output manager, we can print
01002     // debug output.
01003     std::ostream& dbg = outMan_->stream(Debug);
01004     dbg << "Belos::GmresSolMgr::setParametersImpl:" << endl
01005   << "-- Initialized output manager; now we can print debug output." 
01006   << endl;
01007 
01008     // (Re)initialize the stopping criteria.
01009     //
01010     // FIXME (mfh 21 Feb 2011) Currently, userConvTest_ remains unset
01011     // until setUserConvStatusTest() is called.  Should we look for a
01012     // custom test in the parameter list?
01013     rebuildStatusTests (actualParams);
01014     dbg << "-- Initialized status tests." << endl;
01015 
01016     // Initialize the "output status test."  This must be done after
01017     // initializing statusTest_.
01018     {
01019       // FIXME (mfh 16 Feb 2011) Annoyingly, I have to read this back
01020       // in as an int rather than a Belos::OutputType.
01021       //
01022       // TODO (mfh 15 Feb 2011) Validate.
01023       const OutputType outStyle = 
01024   (OutputType) actualParams->get<int> ("Output Style");
01025 
01026       // TODO (mfh 15 Feb 2011) Validate.
01027       const int outFreq = actualParams->get<int> ("Output Frequency");
01028 
01029       // TODO (mfh 15 Feb 2011) Set the solver description according
01030       // to the specific kind of GMRES (CA-GMRES, standard GMRES,
01031       // Flexible GMRES, Flexible CA-GMRES) being used.
01032       const std::string solverDesc (" GMRES ");
01033       // TODO (mfh 15 Feb 2011) Get a preconditioner description.
01034       const std::string precondDesc;
01035       outTest_ = initOutputTest (outMan_, statusTest_, outStyle, outFreq, 
01036          solverDesc, precondDesc);
01037     }
01038     dbg << "-- Initialized output status test." << endl;
01039 
01040     // (Re)initialize the orthogonalization.
01041     rebuildOrthoManager (actualParams);
01042     dbg << "-- Initialized OrthoManager subclass instance." << endl;
01043 
01044     // We don't initialize the Iteration subclass (GmresBaseIteration)
01045     // here; we do that on demand.  This is because GmresBaseIteration
01046     // instantiates GmresBase, and GmresBase wants setLSIndex() to
01047     // have been called on the LinearProblem.  setLSIndex() won't be
01048     // called until the solve() method.
01049 
01050     // Note that the parameter list we store contains the actual
01051     // parameter values used by this solver manager, not necessarily
01052     // those supplied by the user.  We reserve the right to fill in
01053     // default values and silently correct bad values.
01054     params_ = actualParams;
01055 
01056     // If necessary, restart the current solve that might be in
01057     // progress.
01058     if (needToResetSolver)
01059       reset (Belos::Problem);
01060 
01061     dbg << "-- Done with setParametersImpl()." << endl;
01062   }
01063 
01064   template<class Scalar, class MV, class OP>
01065   void
01066   GmresSolMgr<Scalar,MV,OP>::
01067   setUserConvStatusTest (const Teuchos::RCP<StatusTest<Scalar,MV,OP> >& userConvTest)
01068   {
01069     userConvTest_ = userConvTest;
01070 
01071     // For now, we just rebuild the entire stopping criterion from
01072     // scratch.  We could instead keep the parts of the stopping
01073     // criterion that don't change, and just change the user-defined
01074     // convergence criterion.
01075     rebuildStatusTests ();
01076   }
01077 
01078   template<class Scalar, class MV, class OP>
01079   void
01080   GmresSolMgr<Scalar,MV,OP>::
01081   rebuildStatusTests (Teuchos::RCP<Teuchos::ParameterList> plist)
01082   {
01083     using Teuchos::rcp_const_cast;
01084     using Teuchos::ParameterList;
01085     using Teuchos::RCP;
01086 
01087     // Default value for plist is null, in which case we use the
01088     // stored parameter list.  One of those two should be non-null.
01089     RCP<const ParameterList> theParams = plist.is_null() ? params_ : plist;
01090 
01091     TEUCHOS_TEST_FOR_EXCEPTION(theParams.is_null(),
01092            std::logic_error,
01093            "Belos::GmresSolMgr::rebuildStatusTests: We can't (re)"
01094            "build the status tests without any parameters.");
01095     const magnitude_type convTol = 
01096       theParams->get<magnitude_type> ("Convergence Tolerance");
01097     TEUCHOS_TEST_FOR_EXCEPTION(convTol < STM::zero(), std::invalid_argument,
01098            "Convergence tolerance " << convTol << " is negative.");
01099     const int maxIters = 
01100       theParams->get<int> ("Maximum Iterations");
01101     TEUCHOS_TEST_FOR_EXCEPTION(maxIters < 0, std::invalid_argument,
01102            "Maximum number of iterations " << maxIters 
01103            << " is negative.");
01104     // TODO (mfh 15 Feb 2011) Validate.
01105     const std::string implicitScaleType = 
01106       theParams->get<std::string> ("Implicit Residual Scaling");
01107     // TODO (mfh 15 Feb 2011) Validate.
01108     const std::string explicitScaleType = 
01109       theParams->get<std::string> ("Explicit Residual Scaling");
01110 
01111     // If we don't have a problem to solve yet (problem_ is null), do
01112     // both an implicit and an explicit convergence test by default.
01113     // Then, when the user later calls setProblem() with a problem to
01114     // solve, that method will call reset(Belos::Problem), which in
01115     // turn will call rebuildStatusTests() again.
01116     const bool haveLeftPrecond = problem_.is_null() ||
01117       ! problem_->getLeftPrec().is_null();
01118 
01119     convTest_ = initConvTest (convTol, haveLeftPrecond,
01120             implicitScaleType, explicitScaleType,
01121             userConvTest_);
01122     statusTest_ = initStatusTest (convTest_, maxIters);
01123   }
01124 
01125 
01126   template<class Scalar, class MV, class OP>
01127   void
01128   GmresSolMgr<Scalar,MV,OP>::
01129   rebuildStatusTests (const typename GmresSolMgr<Scalar,MV,OP>::magnitude_type convTol,
01130           const int maxItersPerRestart)
01131   {
01132     using Teuchos::ParameterList;
01133     using Teuchos::RCP;
01134 
01135     TEUCHOS_TEST_FOR_EXCEPTION(convTol < STM::zero(), std::invalid_argument,
01136            "Convergence tolerance " << convTol << " is negative.");
01137     params_->set ("Convergence Tolerance", convTol);
01138     TEUCHOS_TEST_FOR_EXCEPTION(maxItersPerRestart < 0, std::invalid_argument,
01139            "Maximum number of iterations " << maxItersPerRestart
01140            << " per restart cycle is negative.");
01141     params_->set ("Maximum Iterations", maxItersPerRestart,
01142       "Maximum number of iterations allowed per restart cycle, "
01143       "for each right-hand side solved.");
01144     // If we don't have a problem to solve yet (problem_ is null), do
01145     // both an implicit and an explicit convergence test by default.
01146     // Then, when the user later calls setProblem() with a problem to
01147     // solve, that method will call reset(Belos::Problem), which in
01148     // turn will call rebuildStatusTests() again.
01149     const bool haveLeftPrecond = problem_.is_null() ||
01150       ! problem_->getLeftPrec().is_null();
01151     const std::string implicitScaleType = 
01152       params_->get<std::string> ("Implicit Residual Scaling");
01153     const std::string explicitScaleType = 
01154       params_->get<std::string> ("Explicit Residual Scaling");
01155     convTest_ = initConvTest (convTol, haveLeftPrecond,
01156             implicitScaleType, explicitScaleType,
01157             userConvTest_);
01158     statusTest_ = initStatusTest (convTest_, maxItersPerRestart);
01159   }
01160 
01161 
01162   template<class Scalar, class MV, class OP>
01163   void
01164   GmresSolMgr<Scalar,MV,OP>::
01165   rebuildOrthoManager (Teuchos::RCP<Teuchos::ParameterList> plist)
01166   {
01167     using Teuchos::null;
01168     using Teuchos::ParameterList;
01169     using Teuchos::parameterList;
01170     using Teuchos::RCP;
01171     using Teuchos::rcp;
01172     using Teuchos::sublist;
01173 
01174     const char prefix[] = "Belos::GmresSolMgr::rebuildOrthoManager: ";
01175     TEUCHOS_TEST_FOR_EXCEPTION(outMan_.is_null(), std::logic_error,
01176            prefix << "Cannot (re)build the orthogonalization "
01177            "method unless the OutputManager has been "
01178            "instantiated.");
01179 
01180     // Default value for plist is null, in which case we use the
01181     // stored parameter list.  One of those two should be non-null.
01182     RCP<ParameterList> actualParams = plist.is_null() ? params_ : plist;
01183 
01184     TEUCHOS_TEST_FOR_EXCEPTION(actualParams.is_null(), std::logic_error,
01185            prefix << "We can't (re)build the orthogonalization "
01186            "method without any parameters.");
01187 
01188     // Get the orthogonalization method name.
01189     const std::string orthoType = 
01190       actualParams->get<std::string> ("Orthogonalization");
01191     
01192     // Validate the orthogonalization method name.
01193     //
01194     // TODO (mfh 16 Feb 2011) Encode the validator in the default
01195     // parameter list.
01196     OrthoManagerFactory<Scalar, MV, OP> factory;
01197     if (! factory.isValidName (orthoType))
01198       {
01199   std::ostringstream os;
01200   os << prefix << "Invalid orthogonalization type \"" << orthoType 
01201      << "\".  Valid orthogonalization types: "
01202      << factory.printValidNames(os) << ".";
01203   throw std::invalid_argument (os.str());
01204       }
01205 
01206     // Get the parameters for that orthogonalization method.
01207     RCP<ParameterList> orthoParams = 
01208       sublist (actualParams, "Orthogonalization Parameters");
01209       
01210     // (Re)instantiate the orthogonalization manager.  Don't bother
01211     // caching this, since it's too much of a pain to check whether
01212     // any of the parameters have changed.
01213 
01214     // Set the timer label for orthogonalization
01215     std::ostringstream os; 
01216     os << "Orthogonalization (method \"" << orthoType << "\")";
01217     orthoMan_ = factory.makeOrthoManager (orthoType, null, outMan_, 
01218             os.str(), orthoParams);
01219   }
01220 
01221   template<class Scalar, class MV, class OP>
01222   void
01223   GmresSolMgr<Scalar,MV,OP>::
01224   rebuildIteration (Teuchos::RCP<Teuchos::ParameterList> plist)
01225   {
01226     using Teuchos::rcp_const_cast;
01227     using Teuchos::ParameterList;
01228     using Teuchos::parameterList;
01229     using Teuchos::RCP;
01230     using Teuchos::rcp;
01231     using Teuchos::sublist;
01232     using std::endl;
01233     const char prefix[] = "Belos::GmresSolMgr::rebuildIteration: ";
01234 
01235     std::ostream& dbg = outMan_->stream(Debug);
01236     dbg << prefix << endl;
01237 
01238     RCP<ParameterList> theParams = plist.is_null() ? params_ : plist;
01239     TEUCHOS_TEST_FOR_EXCEPTION(theParams.is_null(), std::logic_error,
01240            prefix << "We can't (re)build the Iteration subclass "
01241            "instance without any parameters.");
01242     TEUCHOS_TEST_FOR_EXCEPTION(! theParams->isSublist ("Iteration Parameters"),
01243            std::logic_error,
01244            prefix << "The parameter list needs an \"Iteration "
01245            "Parameters\" sublist.");
01246 
01247     RCP<ParameterList> iterParams = sublist (theParams, "Iteration Parameters");
01248     dbg << "-- Instantiating GmresBaseIteration instance:" << endl;
01249     TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::logic_error, 
01250            prefix << "LinearProblem instance is null.");
01251     TEUCHOS_TEST_FOR_EXCEPTION(orthoMan_.is_null(), std::logic_error, 
01252            prefix << "OrthoManager instance is null.");
01253     TEUCHOS_TEST_FOR_EXCEPTION(outMan_.is_null(), std::logic_error, 
01254            prefix << "OutputManager instance is null.");
01255     TEUCHOS_TEST_FOR_EXCEPTION(statusTest_.is_null(), std::logic_error, 
01256            prefix << "StatusTest instance is null.");
01257     TEUCHOS_TEST_FOR_EXCEPTION(iterParams.is_null(), std::logic_error, 
01258            prefix << "Iteration parameters list is null.");
01259     iter_ = rcp (new iteration_type (problem_, orthoMan_, outMan_, 
01260              statusTest_, iterParams));
01261     dbg << "-- Done instantiating GmresBaseIteration instance." << endl;
01262   }
01263 
01264   template<class Scalar, class MV, class OP>
01265   void
01266   GmresSolMgr<Scalar,MV,OP>::reset (const ResetType type) 
01267   {
01268     using Teuchos::rcp;
01269     using std::endl;
01270 
01271     const char prefix[] = "Belos::GmresSolMgr::reset: ";
01272     std::ostream& dbg = outMan_->stream(Debug);
01273     dbg << prefix << endl;
01274 
01275     if ((type & Belos::RecycleSubspace) != 0)
01276       // TODO (mfh 15 Feb 2011) Do we want to support recycling GMRES here?
01277       throw std::invalid_argument ("This version of the GMRES solution "
01278            "manager does not currently support "
01279            "Krylov subspace recycling.");
01280     else if ((type & Belos::Problem) != 0)
01281       { // Recompute the initial residual(s) (both left-
01282   // preconditioned (if applicable) and unpreconditioned).
01283   //
01284   // FIXME (mfh 21 Feb 2011) This recomputes _all_ the initial
01285   // residuals, since it calls setProblem() on the
01286   // LinearProblem instance.  Broken!  We should only compute
01287   // initial residuals for the unsolved right-hand sides.
01288   dbg << "-- Recomputing initial residual(s)" << endl;
01289   problem_ = validatedProblem (problem_);
01290 
01291   // Rebuild the status tests, using the solver manager's
01292   // current list of parameters.
01293   dbg << "-- Rebuilding status tests" << endl;
01294   rebuildStatusTests ();
01295 
01296   // Rebuild the orthogonalization method.  This is a bit of a
01297   // hack; TsqrOrthoManager specifically needs to be
01298   // reinitialized if the data layout of the vectors have
01299   // changed, but the MVT and OPT interfaces don't let me
01300   // check for that.
01301   dbg << "-- Rebuilding OrthoManager" << endl;
01302   rebuildOrthoManager ();
01303       }
01304     else
01305       {
01306   std::ostringstream os;
01307   os << prefix << "Invalid ResetType argument type = " 
01308      << type << ".  Currently, this method only accepts accept type = "
01309     "Belos::Problem (== " << Belos::Problem << ").";
01310   throw std::invalid_argument (os.str());
01311       }
01312   }
01313 
01314   template<class Scalar, class MV, class OP>
01315   ReturnType
01316   GmresSolMgr<Scalar,MV,OP>::solve ()
01317   {
01318     using Teuchos::RCP;
01319     using std::endl;
01320     using std::make_pair;
01321     using std::pair;
01322     using std::vector;
01323 
01324     const char prefix[] = "Belos::GmresSolMgr::solve(): ";
01325     std::ostream& dbg = outMan_->stream(Debug);
01326     dbg << prefix << endl;
01327 
01328     // Reset the status test output.
01329     dbg << "-- Resetting status test output" << endl;
01330     outTest_->reset();
01331     // Reset the number of calls about which the status test output knows.
01332     dbg << "-- Resetting status test number of calls" << endl;
01333     outTest_->resetNumCalls();
01334 
01335     TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::logic_error,
01336            prefix << "The LinearProblem instance is null.");
01337     TEUCHOS_TEST_FOR_EXCEPTION(problem_->getOperator().is_null(), std::logic_error,
01338            prefix << "The LinearProblem's operator (the matrix A "
01339            "in the equation AX=B) is null.");
01340     TEUCHOS_TEST_FOR_EXCEPTION(problem_->getRHS().is_null(), std::logic_error,
01341            prefix << "The LinearProblem's initial guess X is "
01342            "null.");
01343     TEUCHOS_TEST_FOR_EXCEPTION(problem_->getRHS().is_null(), std::logic_error,
01344            prefix << "The LinearProblem's right-hand side B is "
01345            "null.");
01346 
01347     // Set up the linear problem by computing initial residual(s)
01348     // for all right-hand side(s).  setProblem() takes optional
01349     // argument(s) if we want to change the initial guess(es) or the
01350     // right-hand side(s).  At some point we might want to change
01351     // the initial guess(es) for future solve(s), if we know that
01352     // the right-hand side(s) are related.
01353     if (! problem_->isProblemSet())
01354       {
01355   dbg << "-- (Re)computing initial residual(s) for all right-hand "
01356     "side(s)" << endl;
01357   problem_->setProblem ();
01358       }
01359     // Total number of linear systems (right-hand sides) to solve.
01360     const int numRHS = MVT::GetNumberVecs (*(problem_->getRHS()));
01361     dbg << "-- There are " << numRHS << " total right-hand side" 
01362   << (numRHS != 1 ? "s" : "") << " to solve." << endl;
01363 
01364     // Keep track of which linear system(s) converged.  Initially,
01365     // none of them have yet converged, since we haven't tested for
01366     // their convergence yet.
01367     //
01368     // FIXME (mfh 21 Feb 2011) Oh dear, std::vector<bool>.... full
01369     // of hackish badness with respect to iterators.
01370     vector<bool> converged (numRHS, false);
01371     
01372     // For each right-hand side, keep track of the total number of
01373     // restart cycle(s) (first element in the pair), and the total
01374     // number of iterations over all the restart cycle(s) (second
01375     // element in the pair).
01376     vector<pair<int, int> > totalNumIters (numRHS, make_pair (0, 0));
01377 
01378     // Maximum number of restart cycles.
01379     const int maxNumRestarts = params_->get<int> ("Maximum Restarts");
01380     dbg << "-- Max number of restart cycles: " << maxNumRestarts << endl;
01381       
01382     // For each right-hand side, restart until GMRES converges or we
01383     // run out of restart cycles.
01384     for (int curRHS = 0; curRHS < numRHS; ++curRHS)
01385       {
01386   dbg << "-- Solving for right-hand side " << (curRHS+1) 
01387       << " of " << numRHS << ":" << endl;
01388 
01389   // Tell the linear problem for which right-hand side we are
01390   // currently solving.  A block solver could solve for more
01391   // than one right-hand side at a time; the GMRES solvers that
01392   // this solution manager knows how to manage can only solve
01393   // for one right-hand side at a time.
01394   vector<int> curProbIndex (1);
01395   curProbIndex[0] = curRHS;
01396   // This method ensures that problem_->getCurr{LHS,RHS}Vec()
01397   // return the right vector, corresponding to the current
01398   // linear system to solve.
01399   problem_->setLSIndex (curProbIndex);
01400 
01401   // Sanity checks.
01402   {
01403     // Make sure that setLSIndex() has been called on the linear
01404     // problem.  Otherwise, getCurr{L,R}HSVec() (corresponding to
01405     // the current linear system to solve) will return null.
01406     RCP<const MV> X = problem_->getCurrLHSVec();
01407     RCP<const MV> B = problem_->getCurrRHSVec();
01408     TEUCHOS_TEST_FOR_EXCEPTION(X.is_null() || B.is_null(), std::invalid_argument,
01409            prefix << "setLSIndex() must not yet have been "
01410            "called on the given Belos::LinearProblem "
01411            "instance, since getCurrLHSVec() and/or "
01412            "getCurrRHSVec() return null.");
01413     // Our GMRES implementations only know how to solve for one
01414     // right-hand side at a time.  Make sure that X and B each have
01415     // exactly one column.
01416     const int nLHS = MVT::GetNumberVecs(*X);
01417     const int nRHS = MVT::GetNumberVecs(*B);
01418     TEUCHOS_TEST_FOR_EXCEPTION(nLHS != 1 || nRHS != 1, std::invalid_argument,
01419            prefix << "The current linear system to solve has "
01420            << nLHS << " initial guess" 
01421            << (nLHS != 1 ? "es" : "")
01422            << " and " << nRHS << " right-hand side"
01423            << (nLHS != 1 ? "s" : "") 
01424            << ", but our GMRES solvers only know how to "
01425            "solve for one right-hand side at a time.");
01426     if (debug_)
01427       { // In debug mode, compute and print the initial
01428         // residual independently of the solver framework, as
01429         // a sanity check.
01430         RCP<const OP> A = problem_->getOperator ();
01431         RCP<MV> R = MVT::Clone (*B, MVT::GetNumberVecs (*B));
01432         // R := A * X
01433         OPT::Apply (*A, *X, *R);
01434         // R := B - R
01435         MVT::MvAddMv (STS::one(), *B, -STS::one(), *R, *R);
01436         std::vector<magnitude_type> theNorm (MVT::GetNumberVecs (*R));
01437         MVT::MvNorm (*R, theNorm);
01438         dbg << "-- Initial residual norm: ||B - A*X||_2 = " 
01439       << theNorm[0] << endl;
01440       }
01441   }
01442 
01443   // (Re)initialize the GmresBaseIteration, and therefore the
01444   // GmresBase subclass instance, with the current linear
01445   // system to solve.  (We need to do this since we've moved
01446   // on to the next right-hand side, which means we're solving
01447   // a different linear system, even though the
01448   // "LinearProblem" object expresses the same set of linear
01449   // system(s) to solve.
01450   rebuildIteration ();
01451 
01452   // Restart until GMRES converges or we run out of restart
01453   // cycles.
01454   for (int restartCycle = 0; 
01455        convTest_->getStatus() != Passed && restartCycle < maxNumRestarts;
01456        ++restartCycle)
01457     {
01458       dbg << "-- Restart cycle " << (restartCycle+1) << " of " 
01459     << maxNumRestarts << ":" << endl;
01460       // reset() restarts the iteration, so we don't need to
01461       // restart the first time.
01462       if (restartCycle > 0)
01463         iter_->restart ();
01464       // Iterate to convergence or maximum number of
01465       // iterations.
01466       iter_->iterate ();
01467       // For the current right-hand side, keep track of the
01468       // number of restart cycles and the total number of
01469       // iterations over all restart cycles.
01470       totalNumIters[curRHS].first++;
01471       totalNumIters[curRHS].second += iter_->getNumIters();
01472 
01473       // Update the current approximate solution in the linear problem.
01474       iter_->updateSolution ();
01475     }
01476   // Remember whether the current linear system converged.
01477   converged[curRHS] = (convTest_->getStatus() == Passed);
01478 
01479   // Tell the linear problem that we have finished attempting
01480   // to solve the current linear system, and are ready to move
01481   // on to the next system (if there is one).
01482   //
01483   // FIXME (mfh 21 Feb 2011) Yeah, I know, this LinearProblem
01484   // method name doesn't make any sense.  Do _you_ have time to
01485   // rewrite most of Belos?
01486   problem_->setCurrLS();
01487 
01488   // Remember the total number of restarts and iterations.
01489   totalNumIters_ = totalNumIters;
01490       }
01491 
01492     // We've Converged if all of the right-hand sides converged.
01493     if (std::accumulate (converged.begin(), converged.end(), 
01494        true, std::logical_and<bool>()))
01495       return Converged;
01496     else
01497       return Unconverged;
01498   }
01499 
01500 
01501 } // namespace Belos
01502 
01503 #endif // __Belos_GmresSolMgr_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines