BelosGmresInnerSolver.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_GmresInnerSolver_hpp
#define __Belos_GmresInnerSolver_hpp

#include <BelosInnerSolver.hpp>
#include <BelosGmresSolMgr.hpp>
#include <Teuchos_TypeNameTraits.hpp>

namespace Belos {

  template<class Scalar, class MV, class OP>
  class GmresInnerSolver : public InnerSolver<Scalar, MV, OP> {
  public:
    typedef InnerSolver<Scalar, MV, OP> base_type;
    typedef typename base_type::scalar_type scalar_type;
    typedef typename base_type::magnitude_type magnitude_type;
    typedef typename base_type::multivector_type multivector_type;
    typedef typename base_type::operator_type operator_type;

  private:
    static Teuchos::RCP<LinearProblem<Scalar, MV, OP> >
    validProblem (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem)
    {
      if (problem.is_null())
  throw std::invalid_argument("GmresInnerSolver constructor: The linear "
            "problem to solve is null.");
      else if (problem->getOperator().is_null())
  throw std::invalid_argument("GmresInnerSolver constructor: the matrix A"
            " in the linear system AX=B to solve is "
            "null.");
      else
  return problem;
    }

  public:
    GmresInnerSolver (const Teuchos::RCP<LinearProblem<Scalar,MV,OP> >& problem,
          const Teuchos::RCP<const Teuchos::ParameterList>& params,
          const bool debug = false)
      : solMgr_ (validProblem (problem), params, debug)
    {}

    Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const {
      return solMgr_.getCurrentParameters();
    }

    InnerSolveResult
    solve (const Teuchos::RCP<MV>& X,
     const Teuchos::RCP<const MV>& B,
     const magnitude_type convTol,
     const int maxItersPerRestart,
     const int maxNumRestarts)
    {
      solMgr_.setRHS (B);
      solMgr_.setLHS (X);
      // FIXME (mfh 07 Mar 2011) 
      //
      // I have a vague feeling that the stopping criteria should be
      // changed _after_ setting the RHS, but I'm not sure if this is
      // true.  (Does it have something to do with one of the stopping
      // criteria computing or storing its own residual vector(s)?)
      solMgr_.changeStoppingCriteria (convTol, maxItersPerRestart, 
              maxNumRestarts);
      return invokeSolver ();
    }

    InnerSolveResult
    solve (const Teuchos::RCP<MV>& X,
     const Teuchos::RCP<const MV>& B)
    {
      // The solver manager retains the last configuration of
      // convergence tolerance, max number of iterations per restart
      // cycle, and max number of restart cycles. Just set the
      // left-hand side X and right-hand side B, and call solve().
      solMgr_.setLHS (X);
      solMgr_.setRHS (B);
      return invokeSolver ();
    }

  private:
    GmresInnerSolver();

    GmresSolMgr<Scalar, MV, OP> solMgr_;

    InnerSolveResult
    invokeSolver ()
    {
      using std::pair;
      using std::vector;

      ReturnType result = solMgr_.solve ();
      TEUCHOS_TEST_FOR_EXCEPTION(result != Converged && result != Unconverged,
       std::logic_error,
       "The solver manager returned an invalid ResultType "
       "enum value " << result << ".  Valid values are "
       "Converged=" << Converged << " and Unconverged=" 
       << Unconverged << ".");
      //
      // Compute max (total number of iterations, number of restart cycles)
      // for all right-hand side(s).
      //
      const vector<pair<int, int> > iterInfo = solMgr_.totalNumIters();
      // If we could use lambdas, or if C++ had better iteration
      // primitives, we wouldn't need a for loop here.
      typedef vector<pair<int, int>> >::const_iterator iter_type;
      int totalNumIters = 0, numRestartCycles = 0;
      for (iter_type it = iterInfo.begin(); it != iterInfo.end(); ++it)
  {
    numRestartCycles = std::max (it->first, numRestartCycles);
    totalNumIters = std::max (it->second, totalNumIters);
  }
      return InnerSolveResult (result, numRestartCycles, totalNumIters);
    }

  };

} // namespace Belos

#endif // __Belos_GmresInnerSolver_hpp