BelosMinresSolMgr.hpp

Go to the documentation of this file.

//@HEADER
// ************************************************************************
//
//                 Belos: Block Linear Solvers Package
//                  Copyright 2004 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER

#ifndef BELOS_MINRES_SOLMGR_HPP
#define BELOS_MINRES_SOLMGR_HPP


#include "BelosConfigDefs.hpp"
#include "BelosTypes.hpp"

#include "BelosLinearProblem.hpp"
#include "BelosSolverManager.hpp"

#include "BelosMinresIter.hpp"
#include "BelosStatusTestMaxIters.hpp"
#include "BelosStatusTestGenResNorm.hpp"
#include "BelosStatusTestCombo.hpp"
#include "BelosStatusTestOutputFactory.hpp"
#include "BelosOutputManager.hpp"
#include "Teuchos_BLAS.hpp"
#include "Teuchos_LAPACK.hpp"
#ifdef BELOS_TEUCHOS_TIME_MONITOR
#include "Teuchos_TimeMonitor.hpp"
#endif

#include "Teuchos_StandardParameterEntryValidators.hpp"
// Teuchos::ScalarTraits<int> doesn't define rmax(), alas, so we get
// INT_MAX from here.
#include <climits>

namespace Belos {



  //
  class MinresSolMgrLinearProblemFailure : public BelosError {
  public:
    MinresSolMgrLinearProblemFailure (const std::string& what_arg) :
      BelosError(what_arg)
    {}
  };

  template<class ScalarType, class MV, class OP>
  class MinresSolMgr : public SolverManager<ScalarType,MV,OP> {

  private:
    typedef MultiVecTraits<ScalarType,MV> MVT;
    typedef OperatorTraits<ScalarType,MV,OP> OPT;
    typedef Teuchos::ScalarTraits<ScalarType> SCT;
    typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
    typedef Teuchos::ScalarTraits< MagnitudeType > MT;

  public:

    static Teuchos::RCP<const Teuchos::ParameterList> defaultParameters();



    MinresSolMgr();

    MinresSolMgr (const Teuchos::RCP<LinearProblem< ScalarType, MV, OP> > &problem,
                  const Teuchos::RCP<Teuchos::ParameterList> &params);

    virtual ~MinresSolMgr() {};

    Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
      return Teuchos::rcp(new MinresSolMgr<ScalarType,MV,OP>);
    }



    const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
      return *problem_;
    }

    Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override {
      if (defaultParams_.is_null()) {
        defaultParams_ = defaultParameters ();
      }
      return defaultParams_;
    }

    Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override {
      return params_;
    }

    Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
      return Teuchos::tuple (timerSolve_);
    }

    MagnitudeType achievedTol() const override {
      return achievedTol_;
    }

    int getNumIters() const override {
      return numIters_;
    }

    bool isLOADetected() const override { return false; }




    void
    setProblem (const Teuchos::RCP<LinearProblem<ScalarType, MV, OP> > &problem) override
    {
      problem_ = problem;
    }

    void
    setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params) override;




    void
    reset (const ResetType type) override
    {
      if ((type & Belos::Problem) && ! problem_.is_null()) {
        problem_->setProblem ();
      }
    }



    ReturnType solve() override;



    std::string description() const override;


  private:
    Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;

    Teuchos::RCP<OutputManager<ScalarType> > printer_;
    Teuchos::RCP<std::ostream> outputStream_;

    Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;

    Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;

    Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;

    Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > impConvTest_;

    Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_;

    Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;

    mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;

    Teuchos::RCP<Teuchos::ParameterList> params_;

    MagnitudeType convtol_;

    MagnitudeType achievedTol_;

    int maxIters_;

    int numIters_;

    int blockSize_;

    int verbosity_;

    int outputStyle_;

    int outputFreq_;

    std::string label_;

    Teuchos::RCP<Teuchos::Time> timerSolve_;

    bool parametersSet_;

    static void
    validateProblem (const Teuchos::RCP<LinearProblem<ScalarType, MV, OP> >& problem);
  };


  template<class ScalarType, class MV, class OP>
  Teuchos::RCP<const Teuchos::ParameterList>
  MinresSolMgr<ScalarType, MV, OP>::defaultParameters()
  {
    using Teuchos::ParameterList;
    using Teuchos::parameterList;
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::rcpFromRef;
    using Teuchos::EnhancedNumberValidator;
    typedef MagnitudeType MT;
    typedef Teuchos::ScalarTraits<MT> MST;

    // List of parameters accepted by MINRES, and their default values.
    RCP<ParameterList> pl = parameterList ("MINRES");

    pl->set ("Convergence Tolerance", MST::squareroot (MST::eps()),
             "Relative residual tolerance that needs to be achieved by "
             "the iterative solver, in order for the linear system to be "
             "declared converged.",
             rcp (new EnhancedNumberValidator<MT> (MST::zero(), MST::rmax())));
    pl->set ("Maximum Iterations", static_cast<int>(1000),
             "Maximum number of iterations allowed for each right-hand "
             "side solved.",
             rcp (new EnhancedNumberValidator<int> (0, INT_MAX)));
    pl->set ("Num Blocks", static_cast<int> (-1),
             "Ignored, but permitted, for compatibility with other Belos "
             "solvers.");
    pl->set ("Block Size", static_cast<int> (1),
             "Number of vectors in each block.  WARNING: The current "
             "implementation of MINRES only accepts a block size of 1, "
             "since it can only solve for 1 right-hand side at a time.",
             rcp (new EnhancedNumberValidator<int> (1, 1)));
    pl->set ("Verbosity", (int) Belos::Errors,
             "The type(s) of solver information that should "
             "be written to the output stream.");
    pl->set ("Output Style", (int) Belos::General,
             "What style is used for the solver information written "
             "to the output stream.");
    pl->set ("Output Frequency", static_cast<int>(-1),
             "How often (in terms of number of iterations) intermediate "
             "convergence information should be written to the output stream."
             "  -1 means never.");
    pl->set ("Output Stream", rcpFromRef(std::cout),
             "A reference-counted pointer to the output stream where all "
             "solver output is sent.  The output stream defaults to stdout.");
    pl->set ("Timer Label", std::string("Belos"),
             "The string to use as a prefix for the timer labels.");
    return pl;
  }

  //
  // Empty Constructor
  //
  template<class ScalarType, class MV, class OP>
  MinresSolMgr<ScalarType,MV,OP>::MinresSolMgr () :
    convtol_(0.0),
    achievedTol_(0.0),
    maxIters_(0),
    numIters_ (0),
    blockSize_(0),
    verbosity_(0),
    outputStyle_(0),
    outputFreq_(0),
    parametersSet_ (false)
  {}

  //
  // Primary constructor (use this one)
  //
  template<class ScalarType, class MV, class OP>
  MinresSolMgr<ScalarType, MV, OP>::
  MinresSolMgr (const Teuchos::RCP<LinearProblem<ScalarType, MV, OP> > &problem,
                const Teuchos::RCP<Teuchos::ParameterList>& params) :
    problem_ (problem),
    numIters_ (0),
    parametersSet_ (false)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument,
                               "MinresSolMgr: The version of the constructor "
                               "that takes a LinearProblem to solve was given a "
                               "null LinearProblem.");
    setParameters (params);
  }

  template<class ScalarType, class MV, class OP>
  void
  MinresSolMgr<ScalarType, MV, OP>::
  validateProblem (const Teuchos::RCP<LinearProblem<ScalarType, MV, OP> >& problem)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(),
      MinresSolMgrLinearProblemFailure,
      "MINRES requires that you have provided a nonnull LinearProblem to the "
      "solver manager, before you call the solve() method.");
    TEUCHOS_TEST_FOR_EXCEPTION(problem->getOperator().is_null(),
      MinresSolMgrLinearProblemFailure,
      "MINRES requires a LinearProblem object with a non-null operator (the "
      "matrix A).");
    TEUCHOS_TEST_FOR_EXCEPTION(problem->getRHS().is_null(),
      MinresSolMgrLinearProblemFailure,
      "MINRES requires a LinearProblem object with a non-null right-hand side.");
    TEUCHOS_TEST_FOR_EXCEPTION( ! problem->isProblemSet(),
      MinresSolMgrLinearProblemFailure,
      "MINRES requires that before you give it a LinearProblem to solve, you "
      "must first call the linear problem's setProblem() method.");
  }

  template<class ScalarType, class MV, class OP>
  void
  MinresSolMgr<ScalarType, MV, OP>::
  setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
  {
    using Teuchos::ParameterList;
    using Teuchos::parameterList;
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::rcpFromRef;
    using Teuchos::null;
    using Teuchos::is_null;
    using std::string;
    using std::ostream;
    using std::endl;

    if (params_.is_null()) {
      params_ = parameterList (*getValidParameters());
    }
    RCP<ParameterList> pl = params;
    pl->validateParametersAndSetDefaults (*params_);

    //
    // Read parameters from the parameter list.  We have already
    // populated it with defaults.
    //
    blockSize_ = pl->get<int> ("Block Size");
    verbosity_ = pl->get<int> ("Verbosity");
    outputStyle_ = pl->get<int> ("Output Style");
    outputFreq_ = pl->get<int>("Output Frequency");
    outputStream_ = pl->get<RCP<std::ostream> > ("Output Stream");
    convtol_ = pl->get<MagnitudeType> ("Convergence Tolerance");
    maxIters_ = pl->get<int> ("Maximum Iterations");
    //
    // All done reading parameters from the parameter list.
    // Now we know it's valid and we can store it.
    //
    params_ = pl;

    // Change the timer label, and create the timer if necessary.
    const string newLabel = pl->get<string> ("Timer Label");
    {
      if (newLabel != label_ || timerSolve_.is_null()) {
        label_ = newLabel;
#ifdef BELOS_TEUCHOS_TIME_MONITOR
        const string solveLabel = label_ + ": MinresSolMgr total solve time";
        // Unregister the old timer before creating a new one.
        if (! timerSolve_.is_null()) {
          Teuchos::TimeMonitor::clearCounter (label_);
          timerSolve_ = Teuchos::null;
        }
        timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
#endif // BELOS_TEUCHOS_TIME_MONITOR
      }
    }

    // Create output manager, if necessary; otherwise, set its parameters.
    bool recreatedPrinter = false;
    if (printer_.is_null()) {
      printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
      recreatedPrinter = true;
    } else {
      // Set the output stream's verbosity level.
      printer_->setVerbosity (verbosity_);
      // Tell the output manager about the new output stream.
      printer_->setOStream (outputStream_);
    }

    //
    // Set up the convergence tests
    //
    typedef StatusTestGenResNorm<ScalarType, MV, OP> res_norm_type;
    typedef StatusTestCombo<ScalarType, MV, OP> combo_type;

    // Do we need to allocate at least one of the implicit or explicit
    // residual norm convergence tests?
    const bool allocatedConvergenceTests =
      impConvTest_.is_null() || expConvTest_.is_null();

    // Allocate or set the tolerance of the implicit residual norm
    // convergence test.
    if (impConvTest_.is_null()) {
      impConvTest_ = rcp (new res_norm_type (convtol_));
      impConvTest_->defineResForm (res_norm_type::Implicit, TwoNorm);
      // TODO (mfh 03 Nov 2011) Allow users to define the type of
      // scaling (or a custom scaling factor).
      impConvTest_->defineScaleForm (NormOfInitRes, TwoNorm);
    } else {
      impConvTest_->setTolerance (convtol_);
    }

    // Allocate or set the tolerance of the explicit residual norm
    // convergence test.
    if (expConvTest_.is_null()) {
      expConvTest_ = rcp (new res_norm_type (convtol_));
      expConvTest_->defineResForm (res_norm_type::Explicit, TwoNorm);
      // TODO (mfh 03 Nov 2011) Allow users to define the type of
      // scaling (or a custom scaling factor).
      expConvTest_->defineScaleForm (NormOfInitRes, TwoNorm);
    } else {
      expConvTest_->setTolerance (convtol_);
    }

    // Whether we need to recreate the full status test.  We only need
    // to do that if at least one of convTest_ or maxIterTest_ had to
    // be reallocated.
    bool needToRecreateFullStatusTest = sTest_.is_null();

    // Residual status test is a combo of the implicit and explicit
    // convergence tests.
    if (convTest_.is_null() || allocatedConvergenceTests) {
      convTest_ = rcp (new combo_type (combo_type::SEQ, impConvTest_, expConvTest_));
      needToRecreateFullStatusTest = true;
    }

    // Maximum number of iterations status test.  It tells the solver to
    // stop iteration, if the maximum number of iterations has been
    // exceeded.  Initialize it if we haven't yet done so, otherwise
    // tell it the new maximum number of iterations.
    if (maxIterTest_.is_null()) {
      maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
      needToRecreateFullStatusTest = true;
    } else {
      maxIterTest_->setMaxIters (maxIters_);
    }

    // Create the full status test if we need to.
    //
    // The full status test: the maximum number of iterations have
    // been reached, OR the residual has converged.
    //
    // "If we need to" means either that the status test was never
    // created before, or that its two component tests had to be
    // reallocated.
    if (needToRecreateFullStatusTest) {
      sTest_ = rcp (new combo_type (combo_type::OR, maxIterTest_, convTest_));
    }

    // If necessary, create the status test output class.  This class
    // manages and formats the output from the status test.  We have
    // to recreate the output test if we had to (re)allocate either
    // printer_ or sTest_.
    if (outputTest_.is_null() || needToRecreateFullStatusTest || recreatedPrinter) {
      StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
      outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
                                       Passed+Failed+Undefined);
    } else {
      outputTest_->setOutputFrequency (outputFreq_);
    }
    // Set the solver string for the output test.
    // StatusTestOutputFactory has no constructor argument for this.
    outputTest_->setSolverDesc (std::string (" MINRES "));

    // Inform the solver manager that the current parameters were set.
    parametersSet_ = true;

    if (verbosity_ & Debug) {
      using std::endl;

      std::ostream& dbg = printer_->stream (Debug);
      dbg << "MINRES parameters:" << endl << params_ << endl;
    }
  }


  template<class ScalarType, class MV, class OP>
  ReturnType MinresSolMgr<ScalarType,MV,OP>::solve()
  {
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::rcp_const_cast;
    using std::endl;

    if (! parametersSet_) {
      setParameters (params_);
    }
    std::ostream& dbg = printer_->stream (Debug);

#ifdef BELOS_TEUCHOS_TIME_MONITOR
    Teuchos::TimeMonitor solveTimerMonitor (*timerSolve_);
#endif // BELOS_TEUCHOS_TIME_MONITOR

    // We need a problem to solve, else we can't solve it.
    validateProblem (problem_);

    // Reset the status test for this solve.
    outputTest_->reset();

    // The linear problem has this many right-hand sides to solve.
    // MINRES can solve only one at a time, so we solve for each
    // right-hand side in succession.
    const int numRHS2Solve = MVT::GetNumberVecs (*(problem_->getRHS()));

    // Create MINRES iteration object.  Pass along the solver
    // manager's parameters, which have already been validated.
    typedef MinresIter<ScalarType, MV, OP> iter_type;
    RCP<iter_type> minres_iter =
      rcp (new iter_type (problem_, printer_, outputTest_, *params_));

    // The index/indices of the right-hand sides for which MINRES did
    // _not_ converge.  Hopefully this is empty after the for loop
    // below!  If it is not empty, at least one right-hand side did
    // not converge.
    std::vector<int> notConverged;
    std::vector<int> currentIndices(1);

    numIters_ = 0;

    // Solve for each right-hand side in turn.
    for (int currentRHS = 0; currentRHS < numRHS2Solve; ++currentRHS) {
      // Inform the linear problem of the right-hand side(s) currently
      // being solved.  MINRES only knows how to solve linear problems
      // with one right-hand side, so we only include one index, which
      // is the index of the current right-hand side.
      currentIndices[0] = currentRHS;
      problem_->setLSIndex (currentIndices);

      dbg << "-- Current right-hand side index being solved: "
          << currentRHS << endl;

      // Reset the number of iterations.
      minres_iter->resetNumIters();
      // Reset the number of calls that the status test output knows about.
      outputTest_->resetNumCalls();
      // Set the new state and initialize the solver.
      MinresIterationState<ScalarType, MV> newstate;

      // Get the residual vector for the current linear system
      // (that is, for the current right-hand side).
      newstate.Y = MVT::CloneViewNonConst (*(rcp_const_cast<MV> (problem_->getInitResVec())), currentIndices);
      minres_iter->initializeMinres (newstate);

      // Attempt to solve for the solution corresponding to the
      // current right-hand side.
      while (true) {
        try {
          minres_iter->iterate();

          // First check for convergence
          if (convTest_->getStatus() == Passed) {
            dbg << "---- Converged after " << maxIterTest_->getNumIters()
                << " iterations" << endl;
            break;
          }
          // Now check for max # of iterations
          else if (maxIterTest_->getStatus() == Passed) {
            dbg << "---- Did not converge after " << maxIterTest_->getNumIters()
                << " iterations" << endl;
            // This right-hand side didn't converge!
            notConverged.push_back (currentRHS);
            break;
          } else {
            // If we get here, we returned from iterate(), but none of
            // our status tests Passed.  Something is wrong, and it is
            // probably our fault.
            TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
              "Belos::MinresSolMgr::solve(): iterations neither converged, "
              "nor reached the maximum number of iterations " << maxIters_
              << ".  That means something went wrong.");
          }
        } catch (const std::exception &e) {
          printer_->stream (Errors)
            << "Error! Caught std::exception in MinresIter::iterate() at "
            << "iteration " << minres_iter->getNumIters() << endl
            << e.what() << endl;
          throw e;
        }
      }

      // Inform the linear problem that we are finished with the
      // current right-hand side.  It may or may not have converged,
      // but we don't try again if the first time didn't work.
      problem_->setCurrLS();

      // Get iteration information for this solve: total number of
      // iterations for all right-hand sides.
      numIters_ += maxIterTest_->getNumIters();
    }

    // Print final summary of the solution process
    sTest_->print (printer_->stream (FinalSummary));

    // Print timing information, if the corresponding compile-time and
    // run-time options are enabled.
#ifdef BELOS_TEUCHOS_TIME_MONITOR
    // Calling summarize() can be expensive, so don't call unless the
    // user wants to print out timing details.  summarize() will do all
    // the work even if it's passed a "black hole" output stream.
    if (verbosity_ & TimingDetails) {
      Teuchos::TimeMonitor::summarize (printer_->stream (TimingDetails));
    }
#endif // BELOS_TEUCHOS_TIME_MONITOR

    // Save the convergence test value ("achieved tolerance") for this
    // solve.  This solver always has two residual norm status tests:
    // an explicit and an implicit test.  The master convergence test
    // convTest_ is a SEQ combo of the implicit resp. explicit tests.
    // If the implicit test never passes, then the explicit test won't
    // ever be executed.  This manifests as
    // expConvTest_->getTestValue()->size() < 1.  We deal with this
    // case by using the values returned by
    // impConvTest_->getTestValue().
    {
      const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
      if (pTestValues == NULL || pTestValues->size() < 1) {
        pTestValues = impConvTest_->getTestValue();
      }
      TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
        "Belos::MinresSolMgr::solve(): The implicit convergence test's getTestValue() "
        "method returned NULL.  Please report this bug to the Belos developers.");
      TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
        "Belos::MinresSolMgr::solve(): The implicit convergence test's getTestValue() "
        "method returned a vector of length zero.  Please report this bug to the "
        "Belos developers.");

      // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
      // achieved tolerances for all vectors in the current solve(), or
      // just for the vectors from the last deflation?
      achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
    }

    if (notConverged.size() > 0) {
      return Unconverged;
    } else {
      return Converged;
    }
  }

  //  This method requires the solver manager to return a std::string that describes itself.
  template<class ScalarType, class MV, class OP>
  std::string MinresSolMgr<ScalarType,MV,OP>::description() const
  {
    std::ostringstream oss;
    oss << "Belos::MinresSolMgr< "
        << Teuchos::ScalarTraits<ScalarType>::name()
        <<", MV, OP >";
    // oss << "{";
    // oss << "Block Size=" << blockSize_;
    // oss << "}";
    return oss.str();
  }

} // end Belos namespace

#endif /* BELOS_MINRES_SOLMGR_HPP */