MueLu_TekoSmoother_decl.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_TEKOSMOOTHER_DECL_HPP_
#define MUELU_TEKOSMOOTHER_DECL_HPP_
 
#ifdef HAVE_MUELU_TEKO
 
#include "Teko_Utilities.hpp"
 
#include "Teko_InverseLibrary.hpp"
#include "Teko_InverseFactory.hpp"
 
#include "MueLu_ConfigDefs.hpp"
 
#include <Teuchos_ParameterList.hpp>
 
#include <Xpetra_MapExtractor_fwd.hpp>
 
#include "MueLu_TekoSmoother_fwd.hpp"
#include "MueLu_SmootherPrototype.hpp"
#include "MueLu_FactoryBase_fwd.hpp"
#include "MueLu_Monitor.hpp"
 
namespace MueLu {
 
template <class Scalar        = SmootherPrototype<>::scalar_type,
          class LocalOrdinal  = typename SmootherPrototype<Scalar>::local_ordinal_type,
          class GlobalOrdinal = typename SmootherPrototype<Scalar, LocalOrdinal>::global_ordinal_type,
          class Node          = typename SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal>::node_type>
class TekoSmoother : public SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
  typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> MapExtractorClass;
 
#undef MUELU_TEKOSMOOTHER_SHORT
#include "MueLu_UseShortNames.hpp"
 
 public:
 
 
  TekoSmoother()
    : type_("Teko smoother") {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::TekoSmoother: Teko can only be used with SC=double. For more information refer to the doxygen documentation of TekoSmoother.");
  };
 
  TekoSmoother(const Teuchos::ParameterList &paramList)
    : type_("Teko smoother") {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::TekoSmoother: Teko can only be used with SC=double. For more information refer to the doxygen documentation of TekoSmoother.");
  }
 
  virtual ~TekoSmoother() {}
 
 
  RCP<const ParameterList> GetValidParameterList() const {
    RCP<ParameterList> validParamList = rcp(new ParameterList());
    return validParamList;
  }
 
  void DeclareInput(Level &currentLevel) const {}
 
  void SetTekoParameters(RCP<ParameterList> tekoParams){};
 
 
 
  void Setup(Level &currentLevel) {}
 
  void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero = false) const {}
 
  RCP<SmootherPrototype> Copy() const { return Teuchos::null; }
 
 
 
  std::string description() const {
    std::ostringstream out;
    out << SmootherPrototype::description();
    out << "{type = " << type_ << "}";
    return out.str();
  }
 
  // using MueLu::Describable::describe; // overloading, not hiding
  void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel = Default) const {
    MUELU_DESCRIBE;
 
    if (verbLevel & Parameters0)
      out0 << "Prec. type: " << type_ << std::endl;
 
    if (verbLevel & Debug)
      out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
  }
 
  size_t getNodeSmootherComplexity() const;
 
 
 private:
  std::string type_;
};  // class TekoSmoother
 
template <class GlobalOrdinal,
          class Node>
class TekoSmoother<double, int, GlobalOrdinal, Node> : public SmootherPrototype<double, int, GlobalOrdinal, Node> {
  typedef int LocalOrdinal;
  typedef double Scalar;
  typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> MapExtractorClass;
 
#undef MUELU_TEKOSMOOTHER_SHORT
#include "MueLu_UseShortNames.hpp"
 
 public:
 
 
  TekoSmoother()
    : type_("Teko smoother")
    , A_(Teuchos::null)
    , bA_(Teuchos::null)
    , bThyOp_(Teuchos::null)
    , tekoParams_(Teuchos::null)
    , inverseOp_(Teuchos::null){};
 
  TekoSmoother(const Teuchos::ParameterList &paramList)
    : type_("Teko smoother")
    , A_(Teuchos::null)
    , bA_(Teuchos::null)
    , bThyOp_(Teuchos::null)
    , tekoParams_(Teuchos::null)
    , inverseOp_(Teuchos::null) {
    this->SetParameter("Inverse Type", paramList.getEntry("Inverse Type"));
    const auto &tekoSettings = paramList.sublist("Inverse Factory Library");
    tekoParams_              = Teuchos::make_rcp<Teuchos::ParameterList>(tekoSettings);
  };
 
  virtual ~TekoSmoother() {}
 
 
  RCP<const ParameterList> GetValidParameterList() const {
    RCP<ParameterList> validParamList = rcp(new ParameterList());
 
    validParamList->set<RCP<const FactoryBase>>("A", null, "Generating factory of the matrix A");
    validParamList->set<std::string>("Inverse Type", "", "Name of parameter list within 'Teko parameters' containing the Teko smoother parameters.");
 
    return validParamList;
  }
 
  void DeclareInput(Level &currentLevel) const {
    this->Input(currentLevel, "A");
  }
 
  void SetTekoParameters(RCP<ParameterList> tekoParams) { tekoParams_ = tekoParams; };
 
 
 
  void Setup(Level &currentLevel) {
    RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
 
    FactoryMonitor m(*this, "Setup TekoSmoother", currentLevel);
    if (this->IsSetup() == true)
      this->GetOStream(Warnings0) << "MueLu::TekoSmoother::Setup(): Setup() has already been called";
 
    // extract blocked operator A from current level
    A_  = Factory::Get<RCP<Matrix>>(currentLevel, "A");  // A needed for extracting map extractors
    bA_ = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
    TEUCHOS_TEST_FOR_EXCEPTION(bA_.is_null(), Exceptions::BadCast,
                               "MueLu::TekoSmoother::Build: input matrix A is not of type BlockedCrsMatrix.");
 
    bThyOp_ = bA_->getThyraOperator();
    TEUCHOS_TEST_FOR_EXCEPTION(bThyOp_.is_null(), Exceptions::BadCast,
                               "MueLu::TekoSmoother::Build: Could not extract thyra operator from BlockedCrsMatrix.");
 
    Teuchos::RCP<const Thyra::LinearOpBase<Scalar>> thyOp = Teuchos::rcp_dynamic_cast<const Thyra::LinearOpBase<Scalar>>(bThyOp_);
    TEUCHOS_TEST_FOR_EXCEPTION(thyOp.is_null(), Exceptions::BadCast,
                               "MueLu::TekoSmoother::Build: Downcast of Thyra::BlockedLinearOpBase to Teko::LinearOp failed.");
 
    // parameter list contains TekoSmoother parameters but does not handle the Teko parameters itself!
    const ParameterList &pL  = Factory::GetParameterList();
    std::string smootherType = pL.get<std::string>("Inverse Type");
    TEUCHOS_TEST_FOR_EXCEPTION(smootherType.empty(), Exceptions::RuntimeError,
                               "MueLu::TekoSmoother::Build: You must provide a 'Smoother Type' name that is defined in the 'Teko parameters' sublist.");
    type_ = smootherType;
 
    TEUCHOS_TEST_FOR_EXCEPTION(tekoParams_.is_null(), Exceptions::BadCast,
                               "MueLu::TekoSmoother::Build: No Teko parameters have been set.");
 
    Teuchos::RCP<Teko::InverseLibrary> invLib  = Teko::InverseLibrary::buildFromParameterList(*tekoParams_);
    Teuchos::RCP<Teko::InverseFactory> inverse = invLib->getInverseFactory(smootherType);
 
    inverseOp_ = Teko::buildInverse(*inverse, thyOp);
    TEUCHOS_TEST_FOR_EXCEPTION(inverseOp_.is_null(), Exceptions::BadCast,
                               "MueLu::TekoSmoother::Build: Failed to build Teko inverse operator. Probably a problem with the Teko parameters.");
 
    this->IsSetup(true);
  }
 
  void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero = false) const {
    TEUCHOS_TEST_FOR_EXCEPTION(!this->IsSetup(), Exceptions::RuntimeError,
                               "MueLu::TekoSmoother::Apply(): Setup() has not been called");
 
    auto comm           = X.getMap()->getComm();
    auto rgMapExtractor = bA_->getRangeMapExtractor();
    TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgMapExtractor));
 
    using STS = Teuchos::ScalarTraits<Scalar>;
 
    auto createProductMultiVector = [](const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> &space,
                                       int numVecs,
                                       const char *castErrorMsg) {
      auto mv     = Thyra::createMembers(space, numVecs);
      auto prodMv = Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<Scalar>>(mv);
      TEUCHOS_TEST_FOR_EXCEPTION(prodMv.is_null(), Exceptions::BadCast, castErrorMsg);
      return std::make_pair(mv, prodMv);
    };
 
    auto solveWithThyra = [&](const Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &rhs,
                              bool initializeFromX,
                              const typename STS::magnitudeType alphaX) {
      auto [thyB, thyProdB] = createProductMultiVector(
          Teuchos::rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar>>(bThyOp_->productRange()),
          Teuchos::as<int>(rhs->getNumVectors()),
          "MueLu::TekoSmoother::Apply: Failed to cast range space to product range space.");
 
      Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::updateThyra(rhs, rgMapExtractor, thyProdB);
 
      auto [thyX, thyProdX] = createProductMultiVector(
          Teuchos::rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar>>(bThyOp_->productDomain()),
          Teuchos::as<int>(X.getNumVectors()),
          "MueLu::TekoSmoother::Apply: Failed to cast domain space to product domain space.");
 
      if (initializeFromX) {
        Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::updateThyra(Teuchos::rcpFromRef(X), rgMapExtractor, thyProdX);
      }
 
      inverseOp_->apply(
          Thyra::NOTRANS,
          *thyB,
          thyX.ptr(),
          STS::one(),
          STS::zero());
 
      auto XX = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(thyX, comm);
      X.update(STS::one(), *XX, alphaX);
    };
 
    if (InitialGuessIsZero) {
      solveWithThyra(Teuchos::rcpFromRef(B), /*initializeFromX=*/true, STS::zero());
      return;
    }
 
    auto residual = Utilities::Residual(*A_, X, B);
    solveWithThyra(residual, /*initializeFromX=*/false, STS::one());
  }
 
  RCP<SmootherPrototype> Copy() const { return Teuchos::rcp(new MueLu::TekoSmoother<double, int, GlobalOrdinal, Node>(*this)); }
 
 
 
  std::string description() const {
    std::ostringstream out;
    out << SmootherPrototype::description();
    out << "{type = " << type_ << "}";
    return out.str();
  }
 
  // using MueLu::Describable::describe; // overloading, not hiding
  void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel = Default) const {
    MUELU_DESCRIBE;
 
    if (verbLevel & Parameters0)
      out0 << "Prec. type: " << type_ << std::endl;
 
    if (verbLevel & Debug)
      out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
  }
 
  size_t getNodeSmootherComplexity() const {
    size_t cplx = 0;
    return cplx;
  }
 
 
 private:
  std::string type_;
 
  RCP<Matrix> A_;  // < ! internal blocked operator "A" generated by AFact_
  RCP<BlockedCrsMatrix> bA_;
  RCP<const Thyra::BlockedLinearOpBase<Scalar>> bThyOp_;
 
  RCP<ParameterList> tekoParams_;  // < ! parameter list containing Teko parameters. These parameters are not administrated by the factory and not validated.
 
  Teko::LinearOp inverseOp_;  // < ! Teko inverse operator
};                            // class TekoSmoother (specialization on SC=double)
}  // namespace MueLu
 
#define MUELU_TEKOSMOOTHER_SHORT
 
#endif  // HAVE_MUELU_TEKO
 
#endif /* MUELU_TEKOSMOOTHER_DECL_HPP_ */