Belos_Details_LinearSolver.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_DETAILS_LINEARSOLVER_HPP
#define BELOS_DETAILS_LINEARSOLVER_HPP


#include "BelosSolverFactory.hpp"
#include "Trilinos_Details_LinearSolver.hpp"

namespace Belos {
namespace Details {

template<class MV, class OP, class ScalarType, class NormType>
class LinearSolver :
    public Trilinos::Details::LinearSolver<MV, OP, NormType>
{
private:
  typedef Belos::LinearProblem<ScalarType, MV, OP> problem_type;
  typedef Belos::SolverManager<ScalarType, MV, OP> solver_type;

public:

  LinearSolver (const std::string& solverName) :
    solverName_ (solverName)
  {
    // In alignment with Belos philosophy, we delay initialization
    // (which in this case means creation of the solver) until needed.
  }


  virtual ~LinearSolver () {}

  void setMatrix (const Teuchos::RCP<const OP>& A) {
    if (problem_.is_null ()) {
      problem_ = Teuchos::rcp (new problem_type (A, Teuchos::null, Teuchos::null));
    } else if (A != problem_->getOperator ()) {
      problem_->setOperator (A);
    }
  }

  Teuchos::RCP<const OP> getMatrix () const {
    if (problem_.is_null ()) {
      return Teuchos::null;
    } else {
      return problem_->getOperator ();
    }
  }

  void solve (MV& X, const MV& B) {
    TEUCHOS_TEST_FOR_EXCEPTION
      (problem_.is_null () || problem_->getOperator ().is_null (),
       std::runtime_error, "Belos::Details::LinearSolver::solve: "
       "The matrix A in the linear system to solve has not yet been set.  "
       "Please call setMatrix() with a nonnull input before calling solve().");
    Teuchos::RCP<MV> X_ptr = Teuchos::rcpFromRef (X);
    Teuchos::RCP<const MV> B_ptr = Teuchos::rcpFromRef (B);

    problem_->setLHS (X_ptr);
    problem_->setRHS (B_ptr);
    problem_->setProblem ();

    // We can delay creating the Belos solver until the moment when we
    // actually need it.  This aligns with Belos' preference for lazy
    // initialization.
    if (solver_.is_null ()) {
      Belos::SolverFactory<ScalarType, MV, OP> factory;
      solver_ = factory.create (solverName_, params_);
      solver_->setProblem (problem_);
    }

    //Belos::ReturnType ret = solver_->solve ();
    (void) solver_->solve ();
  }

  void setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params) {
    if (! solver_.is_null () && ! params.is_null ()) {
      solver_->setParameters (params);
    }
    params_ = params;
  }

  void symbolic () {
    TEUCHOS_TEST_FOR_EXCEPTION
      (problem_.is_null () || problem_->getOperator ().is_null (),
       std::runtime_error, "Belos::Details::LinearSolver::symbolic: "
       "The matrix A in the linear system to solve has not yet been set.  "
       "Please call setMatrix() with a nonnull input before calling this method.");

    // Belos solvers can't handle changes to the matrix's domain or
    // range Maps.  It's best in this case to destroy the solver and
    // start over.
    solver_ = Teuchos::null;
  }

  void numeric () {
    TEUCHOS_TEST_FOR_EXCEPTION
      (problem_.is_null () || problem_->getOperator ().is_null (),
       std::runtime_error, "Belos::Details::LinearSolver::numeric: "
       "The matrix A in the linear system to solve has not yet been set.  "
       "Please call setMatrix() with a nonnull input before calling this method.");
    // NOTE (mfh 23 Aug 2015) For the seed or recycling solvers, it
    // would make sense to do something special here.  However, the
    // line below is always correct.  It recomputes the initial
    // residual vector, which is what Belos expects before a solve if
    // the matrix or right-hand side may have changed.
    problem_->setProblem ();
  }

private:
  std::string solverName_;
  Teuchos::RCP<problem_type> problem_;
  Teuchos::RCP<solver_type> solver_;
  Teuchos::RCP<Teuchos::ParameterList> params_;
};

} // namespace Details
} // namespace Belos

#endif /* BELOS_DETAILS_LINEARSOLVER_HPP */