Tpetra_MixedScalarMultiplyOp.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef TPETRA_MIXEDSCALARMULTIPLYOP_HPP
#define TPETRA_MIXEDSCALARMULTIPLYOP_HPP


#include "Tpetra_Operator.hpp"
#include "Tpetra_Util.hpp"
#include "Tpetra_Details_Behavior.hpp"
#include "Tpetra_Details_Profiling.hpp"

namespace Tpetra {

template <class Scalar,
          class OpScalar,
          class LocalOrdinal,
          class GlobalOrdinal,
          class Node>
class MixedScalarMultiplyOp : public Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
 public:
  using op_type =
      Operator<OpScalar, LocalOrdinal, GlobalOrdinal, Node>;
  using map_type = Map<LocalOrdinal, GlobalOrdinal, Node>;

 public:


  MixedScalarMultiplyOp(const Teuchos::RCP<const op_type>& A)
    : op_(A) {
    Allocate(1);
  }

  ~MixedScalarMultiplyOp() override = default;



  void
  apply(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& X,
        MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Y,
        Teuchos::ETransp mode = Teuchos::NO_TRANS,
        Scalar alpha          = Teuchos::ScalarTraits<Scalar>::one(),
        Scalar beta           = Teuchos::ScalarTraits<Scalar>::zero()) const override {
    TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
                               Teuchos::typeName(*this) << "::apply(X,Y): X and Y must have the same number of vectors.");

    if (X.getNumVectors() != X_->getNumVectors()) {
      Allocate(X.getNumVectors());
    }

    Tpetra::deep_copy(*X_, X);
    op_->apply(*X_, *Y_, mode, Teuchos::as<OpScalar>(alpha), Teuchos::as<OpScalar>(beta));
    Tpetra::deep_copy(Y, *Y_);
  }

  bool hasTransposeApply() const override {
    return true;
  }

  Teuchos::RCP<const map_type> getDomainMap() const override {
    return op_->getDomainMap();
  }

  Teuchos::RCP<const map_type> getRangeMap() const override {
    return op_->getRangeMap();
  }


 protected:
  typedef MultiVector<OpScalar, LocalOrdinal, GlobalOrdinal, Node> MV;

  const Teuchos::RCP<const op_type> op_;
  mutable Teuchos::RCP<MV> X_, Y_;

  void
  Allocate(int numVecs) const {
    X_ = rcp(new MV(op_->getDomainMap(), numVecs));
    Y_ = rcp(new MV(op_->getRangeMap(), numVecs));
  }
};

template <class Scalar,
          class OpScalar,
          class LocalOrdinal,
          class GlobalOrdinal,
          class Node>
Teuchos::RCP<
    MixedScalarMultiplyOp<Scalar, OpScalar, LocalOrdinal, GlobalOrdinal, Node> >
createMixedScalarMultiplyOp(const Teuchos::RCP<
                            const Operator<OpScalar, LocalOrdinal, GlobalOrdinal, Node> >& A) {
  typedef MixedScalarMultiplyOp<Scalar, OpScalar, LocalOrdinal,
                                GlobalOrdinal, Node>
      op_type;
  return Teuchos::rcp(new op_type(A));
}

}  // end of namespace Tpetra

#endif  // TPETRA_MIXEDSCALARMULTIPLYOP_HPP