MueLu_Details_LinearSolverFactory_def.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//        MueLu: A package for multigrid based preconditioning
//
// Copyright 2012 NTESS and the MueLu contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER


#ifndef MUELU_DETAILS_LINEARSOLVERFACTORY_DEF_HPP
#define MUELU_DETAILS_LINEARSOLVERFACTORY_DEF_HPP

#include "MueLu_config.hpp"
#include "Trilinos_Details_LinearSolver.hpp"
#include "Trilinos_Details_LinearSolverFactory.hpp"
#include <type_traits>

#ifdef HAVE_MUELU_EPETRA
#include "Epetra_CrsMatrix.h"
#include "MueLu_CreateEpetraPreconditioner.hpp"
#endif  // HAVE_MUELU_EPETRA

#include "Tpetra_Operator.hpp"
#include "MueLu_CreateTpetraPreconditioner.hpp"

namespace MueLu {
namespace Details {

template <class MV, class OP, class NormType>
class LinearSolver : public Trilinos::Details::LinearSolver<MV, OP, NormType>,
                     virtual public Teuchos::Describable {
 public:
  LinearSolver() {}

  virtual ~LinearSolver() {}

  void setMatrix(const Teuchos::RCP<const OP>& A);

  Teuchos::RCP<const OP> getMatrix() const {
    return A_;
  }

  void solve(MV& X, const MV& B);

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

  void symbolic() {}

  void numeric();

  std::string description() const;

  void
  describe(Teuchos::FancyOStream& out,
           const Teuchos::EVerbosityLevel verbLevel =
               Teuchos::Describable::verbLevel_default) const;

 private:
  Teuchos::RCP<const OP> A_;
  Teuchos::RCP<Teuchos::ParameterList> params_;
};

// Why does MueLu_EpetraOperator insist on HAVE_MUELU_SERIAL?
#if defined(HAVE_MUELU_SERIAL) and defined(HAVE_MUELU_EPETRA)
template <>
class LinearSolver<Epetra_MultiVector, Epetra_Operator, double> : public Trilinos::Details::LinearSolver<Epetra_MultiVector, Epetra_Operator, double>,
                                                                  virtual public Teuchos::Describable {
 public:
  LinearSolver()
    : changedA_(false)
    , changedParams_(false) {}

  virtual ~LinearSolver() {}

  void setMatrix(const Teuchos::RCP<const Epetra_Operator>& A) {
    const char prefix[] = "MueLu::Details::LinearSolver::setMatrix: ";

    if (A != A_) {
      if (solver_ != Teuchos::null)
        changedA_ = true;

      A_ = rcp_dynamic_cast<const Epetra_CrsMatrix>(A);
      TEUCHOS_TEST_FOR_EXCEPTION(A_.is_null(), std::runtime_error, prefix << "MueLu requires "
                                                                             "an Epetra_CrsMatrix, but the matrix you provided is of a "
                                                                             "different type.  Please provide an Epetra_CrsMatrix instead.");
    }
  }

  Teuchos::RCP<const Epetra_Operator> getMatrix() const {
    return A_;
  }

  void solve(Epetra_MultiVector& X, const Epetra_MultiVector& B) {
    // TODO amk: Do we assume the user has called numeric before solve, or should we call it for them?
    const char prefix[] = "MueLu::Details::LinearSolver::solve: ";
    TEUCHOS_TEST_FOR_EXCEPTION(solver_.is_null(), std::runtime_error, prefix << "The solver does not "
                                                                                "exist yet.  You must call numeric() before you may call this method.");
    TEUCHOS_TEST_FOR_EXCEPTION(changedA_, std::runtime_error, prefix << "The matrix A has been reset "
                                                                        "since the last call to numeric().  Please call numeric() again.");
    TEUCHOS_TEST_FOR_EXCEPTION(changedParams_, std::runtime_error, prefix << "The parameters have been reset "
                                                                             "since the last call to numeric().  Please call numeric() again.");

    int err = solver_->ApplyInverse(B, X);

    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, prefix << "EpetraOperator::ApplyInverse returned "
                                                                       "nonzero error code "
                                                                    << err);
  }

  void setParameters(const Teuchos::RCP<Teuchos::ParameterList>& params) {
    if (solver_ != Teuchos::null && params != params_)
      changedParams_ = true;

    params_ = params;
  }

  void symbolic() {}

  void numeric() {
    const char prefix[] = "MueLu::Details::LinearSolver::numeric: ";

    // If the solver is up-to-date, leave it alone
    if (solver_ == Teuchos::null || changedA_ || changedParams_) {
      changedA_      = false;
      changedParams_ = false;

      TEUCHOS_TEST_FOR_EXCEPTION(A_ == Teuchos::null, std::runtime_error, prefix << "The matrix has not been "
                                                                                    "set yet.  You must call setMatrix() with a nonnull matrix before you may "
                                                                                    "call this method.");

      // TODO: We should not have to cast away the constness here
      // TODO: See bug 6462
      if (params_ != Teuchos::null)
        solver_ = CreateEpetraPreconditioner(rcp_const_cast<Epetra_CrsMatrix>(A_), *params_);
      else
        solver_ = CreateEpetraPreconditioner(rcp_const_cast<Epetra_CrsMatrix>(A_));
    }
  }

  std::string description() const {
    if (solver_.is_null()) {
      return "\"MueLu::Details::LinearSolver\": {MV: Epetra_MultiVector, OP: Epetra_Operator, NormType: double}";
    } else {
      return solver_->GetHierarchy()->description();
    }
  }

  void
  describe(Teuchos::FancyOStream& out,
           const Teuchos::EVerbosityLevel verbLevel =
               Teuchos::Describable::verbLevel_default) const {
    using std::endl;
    if (solver_.is_null()) {
      if (verbLevel > Teuchos::VERB_NONE) {
        Teuchos::OSTab tab0(out);
        out << "\"MueLu::Details::LinearSolver\":" << endl;
        Teuchos::OSTab tab1(out);
        out << "MV: Epetra_MultiVector" << endl
            << "OP: Epetra_Operator" << endl
            << "NormType: double" << endl;
      }
    } else {
      solver_->GetHierarchy()->describe(out, verbLevel);
    }
  }

 private:
  Teuchos::RCP<const Epetra_CrsMatrix> A_;
  Teuchos::RCP<Teuchos::ParameterList> params_;
  Teuchos::RCP<EpetraOperator> solver_;
  bool changedA_;
  bool changedParams_;
};
#endif  // HAVE_MUELU_EPETRA

template <class Scalar, class LO, class GO, class Node>
class LinearSolver<Tpetra::MultiVector<Scalar, LO, GO, Node>,
                   Tpetra::Operator<Scalar, LO, GO, Node>,
                   typename Teuchos::ScalarTraits<Scalar>::magnitudeType> : public Trilinos::Details::LinearSolver<Tpetra::MultiVector<Scalar, LO, GO, Node>,
                                                                                                                   Tpetra::Operator<Scalar, LO, GO, Node>,
                                                                                                                   typename Teuchos::ScalarTraits<Scalar>::magnitudeType>,
                                                                            virtual public Teuchos::Describable {
 public:
  LinearSolver()
    : changedA_(false)
    , changedParams_(false) {}

  virtual ~LinearSolver() {}

  void setMatrix(const Teuchos::RCP<const Tpetra::Operator<Scalar, LO, GO, Node> >& A) {
    if (A != A_) {
      if (solver_ != Teuchos::null)
        changedA_ = true;

      A_ = A;
    }
  }

  Teuchos::RCP<const Tpetra::Operator<Scalar, LO, GO, Node> > getMatrix() const {
    return A_;
  }

  void solve(Tpetra::MultiVector<Scalar, LO, GO, Node>& X, const Tpetra::MultiVector<Scalar, LO, GO, Node>& B) {
    // TODO amk: Do we assume the user has called numeric before solve, or should we call it for them?
    const char prefix[] = "MueLu::Details::LinearSolver::solve: ";
    TEUCHOS_TEST_FOR_EXCEPTION(solver_.is_null(), std::runtime_error, prefix << "The solver does not "
                                                                                "exist yet.  You must call numeric() before you may call this method.");
    TEUCHOS_TEST_FOR_EXCEPTION(changedA_, std::runtime_error, prefix << "The matrix A has been reset "
                                                                        "since the last call to numeric().  Please call numeric() again.");
    TEUCHOS_TEST_FOR_EXCEPTION(changedParams_, std::runtime_error, prefix << "The parameters have been reset "
                                                                             "since the last call to numeric().  Please call numeric() again.");

    solver_->apply(B, X);
  }

  void setParameters(const Teuchos::RCP<Teuchos::ParameterList>& params) {
    if (solver_ != Teuchos::null && params != params_)
      changedParams_ = true;

    params_ = params;
  }

  void symbolic() {}

  void numeric() {
    const char prefix[] = "MueLu::Details::LinearSolver::numeric: ";

    // If the solver is up-to-date, leave it alone
    if (solver_ == Teuchos::null || changedParams_) {
      TEUCHOS_TEST_FOR_EXCEPTION(A_ == Teuchos::null, std::runtime_error, prefix << "The matrix has not been "
                                                                                    "set yet.  You must call setMatrix() with a nonnull matrix before you may "
                                                                                    "call this method.");

      // TODO: We should not have to cast away the constness here
      // TODO: See bug 6462
      if (params_ != Teuchos::null)
        solver_ = CreateTpetraPreconditioner(rcp_const_cast<Tpetra::Operator<Scalar, LO, GO, Node> >(A_), *params_);
      else
        solver_ = CreateTpetraPreconditioner(rcp_const_cast<Tpetra::Operator<Scalar, LO, GO, Node> >(A_));
    } else if (changedA_) {
      TEUCHOS_TEST_FOR_EXCEPTION(A_ == Teuchos::null, std::runtime_error, prefix << "The matrix has not been "
                                                                                    "set yet.  You must call setMatrix() with a nonnull matrix before you may "
                                                                                    "call this method.");

      // TODO: We should not have to cast away the constness here
      // TODO: See bug 6462
      RCP<const Tpetra::CrsMatrix<Scalar, LO, GO, Node> > helperMat;
      helperMat = rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LO, GO, Node> >(A_);
      TEUCHOS_TEST_FOR_EXCEPTION(helperMat.is_null(), std::runtime_error, prefix << "MueLu requires "
                                                                                    "a Tpetra::CrsMatrix, but the matrix you provided is of a "
                                                                                    "different type.  Please provide a Tpetra::CrsMatrix instead.");
      ReuseTpetraPreconditioner(rcp_const_cast<Tpetra::CrsMatrix<Scalar, LO, GO, Node> >(helperMat), *solver_);
    }

    changedA_      = false;
    changedParams_ = false;
  }

  std::string description() const {
    using Teuchos::TypeNameTraits;
    if (solver_.is_null()) {
      std::ostringstream os;
      os << "\"MueLu::Details::LinearSolver\": {"
         << "MV: " << TypeNameTraits<Tpetra::MultiVector<Scalar, LO, GO, Node> >::name()
         << "OP: " << TypeNameTraits<Tpetra::Operator<Scalar, LO, GO, Node> >::name()
         << "NormType: " << TypeNameTraits<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>::name()
         << "}";
      return os.str();
    } else {
      return solver_->GetHierarchy()->description();
    }
  }

  void
  describe(Teuchos::FancyOStream& out,
           const Teuchos::EVerbosityLevel verbLevel =
               Teuchos::Describable::verbLevel_default) const {
    using std::endl;
    using Teuchos::TypeNameTraits;
    if (solver_.is_null()) {
      if (verbLevel > Teuchos::VERB_NONE) {
        Teuchos::OSTab tab0(out);
        out << "\"MueLu::Details::LinearSolver\":" << endl;
        Teuchos::OSTab tab1(out);
        out << "MV: " << TypeNameTraits<Tpetra::MultiVector<Scalar, LO, GO, Node> >::name() << endl
            << "OP: " << TypeNameTraits<Tpetra::Operator<Scalar, LO, GO, Node> >::name() << endl
            << "NormType: " << TypeNameTraits<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>::name() << endl;
      }
    } else {
      solver_->GetHierarchy()->describe(out, verbLevel);
    }
  }

 private:
  Teuchos::RCP<const Tpetra::Operator<Scalar, LO, GO, Node> > A_;
  Teuchos::RCP<Teuchos::ParameterList> params_;
  Teuchos::RCP<TpetraOperator<Scalar, LO, GO, Node> > solver_;
  bool changedA_;
  bool changedParams_;
};

template <class MV, class OP, class NormType>
Teuchos::RCP<Trilinos::Details::LinearSolver<MV, OP, NormType> >
LinearSolverFactory<MV, OP, NormType>::
    getLinearSolver(const std::string& solverName) {
  using Teuchos::rcp;
  return rcp(new MueLu::Details::LinearSolver<MV, OP, NormType>());
}

template <class MV, class OP, class NormType>
void LinearSolverFactory<MV, OP, NormType>::
    registerLinearSolverFactory() {
#ifdef HAVE_TEUCHOSCORE_CXX11
  typedef std::shared_ptr<MueLu::Details::LinearSolverFactory<MV, OP, NormType> > ptr_type;
  // typedef std::shared_ptr<Trilinos::Details::LinearSolverFactory<MV, OP> > base_ptr_type;
#else
  typedef Teuchos::RCP<MueLu::Details::LinearSolverFactory<MV, OP, NormType> > ptr_type;
  // typedef Teuchos::RCP<Trilinos::Details::LinearSolverFactory<MV, OP> > base_ptr_type;
#endif  // HAVE_TEUCHOSCORE_CXX11

  ptr_type factory(new MueLu::Details::LinearSolverFactory<MV, OP, NormType>());
  Trilinos::Details::registerLinearSolverFactory<MV, OP, NormType>("MueLu", factory);
}

}  // namespace Details
}  // namespace MueLu

// Macro for doing explicit instantiation of
// MueLu::Details::LinearSolverFactory, for Tpetra objects, with
// given Tpetra template parameters (SC = Scalar, LO = LocalOrdinal,
// GO = GlobalOrdinal, NT = Node).
//
// We don't have to protect use of Tpetra objects here, or include
// any header files for them, because this is a macro definition.
#define MUELU_DETAILS_LINEARSOLVERFACTORY_INSTANT(SC, LO, GO, NT)                         \
  template class MueLu::Details::LinearSolverFactory<Tpetra::MultiVector<SC, LO, GO, NT>, \
                                                     Tpetra::Operator<SC, LO, GO, NT>,    \
                                                     typename Tpetra::MultiVector<SC, LO, GO, NT>::mag_type>;

#endif  // MUELU_DETAILS_LINEARSOLVERFACTORY_DEF_HPP