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() {};
    


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

    Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
      if (defaultParams_.is_null()) {
  defaultParams_ = defaultParameters ();
      }
      return defaultParams_;
    }
    
    Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { 
      return params_; 
    }
    
    Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
      return Teuchos::tuple (timerSolve_);
    }

    MagnitudeType achievedTol() const {
      return achievedTol_;
    }

    int getNumIters() const {
      return numIters_;
    }

    bool isLOADetected() const { return false; }
 
    

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


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

    
    ReturnType solve();
    
    
    
    std::string description() const;
    
    
  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;
    typedef Teuchos::ScalarTraits<int> IST;

    // 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 ("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 () :
    numIters_ (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;

    RCP<const ParameterList> defaults = getValidParameters ();

    RCP<ParameterList> pl;
    if (params.is_null()) {
      // We don't need to validate the default parameter values.
      pl = parameterList (*defaults);
    } else {
      pl = params;
      pl->validateParametersAndSetDefaults (*defaults);
    }

    //
    // 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 */