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