42#ifndef TPETRA_CRSMATRIXMULTIPLYOP_HPP
43#define TPETRA_CRSMATRIXMULTIPLYOP_HPP
51#include "Tpetra_CrsMatrix.hpp"
55#include "Tpetra_LocalCrsMatrixOperator.hpp"
95 template <
class Scalar,
101 public Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>
111 using local_matrix_device_type =
125 A->getLocalMatrixDevice ()))
140 Teuchos::ETransp mode = Teuchos::NO_TRANS,
141 Scalar alpha = Teuchos::ScalarTraits<Scalar>::one (),
142 Scalar beta = Teuchos::ScalarTraits<Scalar>::zero ())
const override
144 TEUCHOS_TEST_FOR_EXCEPTION
145 (!
matrix_->isFillComplete (), std::runtime_error,
146 Teuchos::typeName (*
this) <<
"::apply(): underlying matrix is not fill-complete.");
147 TEUCHOS_TEST_FOR_EXCEPTION
149 Teuchos::typeName (*
this) <<
"::apply(X,Y): X and Y must have the same number of vectors.");
150 TEUCHOS_TEST_FOR_EXCEPTION
151 (Teuchos::ScalarTraits<Scalar>::isComplex && mode == Teuchos::TRANS, std::logic_error,
152 Teuchos::typeName (*
this) <<
"::apply() does not currently support transposed multiplications for complex scalar types.");
153 if (mode == Teuchos::NO_TRANS) {
172 return matrix_->getDomainMap ();
177 return matrix_->getRangeMap ();
186 const Teuchos::RCP<const crs_matrix_type>
matrix_;
204 mutable Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
importMV_;
218 mutable Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
exportMV_;
225 Teuchos::ETransp mode,
234 using STS = Teuchos::ScalarTraits<Scalar>;
241 RCP<const import_type> importer =
matrix_->getGraph()->getImporter();
242 RCP<const export_type> exporter =
matrix_->getGraph()->getExporter();
246 const bool Y_is_overwritten = (beta == STS::zero ());
247 const int myRank =
matrix_->getComm ()->getRank ();
248 if (Y_is_replicated && myRank != 0) {
257 X = Teuchos::rcp (
new MV (X_in, Teuchos::Copy));
261 X = Teuchos::rcpFromRef (X_in);
265 if (importer != null) {
273 if (exporter != null) {
283 if (exporter != null) {
288 auto X_lcl = X->getLocalViewDevice(Access::ReadOnly);
292 if (importer != null) {
294 auto Y_lcl =
importMV_->getLocalViewDevice(Access::OverwriteAll);
297 if (Y_is_overwritten) {
310 MV Y (Y_in, Teuchos::Copy);
311 nonconst_view_type Y_lcl;
319 nonconst_view_type Y_lcl;
327 if (Y_is_replicated) {
340 using Teuchos::NO_TRANS;
343 using Teuchos::rcp_const_cast;
344 using Teuchos::rcpFromRef;
347 typedef Teuchos::ScalarTraits<Scalar> STS;
350 if (alpha == STS::zero ()) {
351 if (beta == STS::zero ()) {
353 }
else if (beta != STS::one ()) {
368 RCP<const import_type> importer =
matrix_->getGraph()->getImporter();
369 RCP<const export_type> exporter =
matrix_->getGraph()->getExporter();
375 const bool Y_is_overwritten = (beta == STS::zero());
378 const bool Y_is_replicated =
388 if (Y_is_replicated &&
matrix_->getComm ()->getRank () > 0) {
395 RCP<const MV> X_colMap;
396 if (importer.is_null ()) {
404 RCP<MV> X_colMapNonConst = getColumnMapMultiVector (X_in,
true);
406 X_colMap = rcp_const_cast<const MV> (X_colMapNonConst);
411 X_colMap = rcpFromRef (X_in);
415 ProfilingRegion regionImport (
"Tpetra::CrsMatrixMultiplyOp::apply: Import");
421 RCP<MV> X_colMapNonConst = getColumnMapMultiVector (X_in);
424 X_colMapNonConst->doImport (X_in, *importer,
INSERT);
425 X_colMap = rcp_const_cast<const MV> (X_colMapNonConst);
432 RCP<MV> Y_rowMap = getRowMapMultiVector (Y_in);
434 auto X_lcl = X_colMap->getLocalViewDevice(Access::ReadOnly);
441 if (! exporter.is_null ()) {
442 auto Y_lcl = Y_rowMap->getLocalViewDevice(Access::OverwriteAll);
445 alpha, STS::zero ());
447 ProfilingRegion regionExport (
"Tpetra::CrsMatrixMultiplyOp::apply: Export");
453 if (Y_is_overwritten) {
480 Y_rowMap = getRowMapMultiVector (Y_in,
true);
484 if (beta != STS::zero ()) {
487 nonconst_view_type Y_lcl;
488 if(Y_is_overwritten) Y_lcl = Y_rowMap->getLocalViewDevice(Access::OverwriteAll);
489 else Y_lcl = Y_rowMap->getLocalViewDevice(Access::ReadWrite);
495 nonconst_view_type Y_lcl;
507 if (Y_is_replicated) {
508 ProfilingRegion regionReduce (
"Tpetra::CrsMatrixMultiplyOp::apply: Reduce Y");
535 const bool force =
false)
const
543 RCP<const import_type> importer =
matrix_->getGraph ()->getImporter ();
544 RCP<const map_type> colMap =
matrix_->getColMap ();
557 if (! importer.is_null () || force) {
559 X_colMap = rcp (
new MV (colMap, numVecs));
598 getRowMapMultiVector (
const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y_rangeMap,
599 const bool force =
false)
const
604 typedef Export<LocalOrdinal,GlobalOrdinal,Node> export_type;
606 const size_t numVecs = Y_rangeMap.getNumVectors ();
607 RCP<const export_type> exporter =
matrix_->getGraph ()->getExporter ();
608 RCP<const map_type> rowMap =
matrix_->getRowMap ();
620 if (! exporter.is_null () || force) {
622 Y_rowMap = rcp (
new MV (rowMap, numVecs));
640 template <
class OpScalar,
646 CrsMatrixMultiplyOp<OpScalar, MatScalar, LocalOrdinal, GlobalOrdinal, Node> >
651 GlobalOrdinal, Node> op_type;
652 return Teuchos::rcp (
new op_type (A));
Forward declaration of Tpetra::CrsMatrixMultiplyOp.
Declaration of Tpetra::Details::Behavior, a class that describes Tpetra's behavior.
Declaration of Tpetra::Details::Profiling, a scope guard for Kokkos Profiling.
Stand-alone utility functions and macros.
A class for wrapping a CrsMatrix multiply in a Operator.
Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > importMV_
Column Map MultiVector used in apply().
bool hasTransposeApply() const override
Whether this Operator's apply() method can apply the transpose or conjugate transpose.
void applyNonTranspose(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X_in, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y_in, Scalar alpha, Scalar beta) const
Apply the matrix (not its transpose) to X_in, producing Y_in.
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
Compute Y = beta*Y + alpha*Op(A)*X, where Op(A) is either A, , or .
Teuchos::RCP< const map_type > getRangeMap() const override
The range Map of this Operator.
void applyTranspose(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X_in, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y_in, Teuchos::ETransp mode, Scalar alpha, Scalar beta) const
Apply the transpose or conjugate transpose of the matrix to X_in, producing Y_in.
CrsMatrixMultiplyOp(const Teuchos::RCP< const crs_matrix_type > &A)
Constructor.
~CrsMatrixMultiplyOp() override=default
Destructor (virtual for memory safety of derived classes).
Teuchos::RCP< CrsMatrixMultiplyOp< OpScalar, MatScalar, LocalOrdinal, GlobalOrdinal, Node > > createCrsMatrixMultiplyOp(const Teuchos::RCP< const CrsMatrix< MatScalar, LocalOrdinal, GlobalOrdinal, Node > > &A)
Non-member function to create a CrsMatrixMultiplyOp.
const Teuchos::RCP< const crs_matrix_type > matrix_
The underlying CrsMatrix object.
Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > exportMV_
Row Map MultiVector used in apply().
LocalCrsMatrixOperator< Scalar, MatScalar, typename crs_matrix_type::device_type > localMultiply_
Implementation of local sparse matrix-vector multiply.
Teuchos::RCP< const map_type > getDomainMap() const override
The domain Map of this Operator.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries.
typename Node::device_type device_type
The Kokkos device type.
KokkosSparse::CrsMatrix< impl_scalar_type, local_ordinal_type, device_type, void, typename local_graph_device_type::size_type > local_matrix_device_type
The specialization of Kokkos::CrsMatrix that represents the part of the sparse matrix on each MPI pro...
void doExport(const SrcDistObject &source, const Export< LocalOrdinal, GlobalOrdinal, Node > &exporter, const CombineMode CM, const bool restrictedMode=false)
Export data into this object using an Export object ("forward mode").
bool isDistributed() const
Whether this is a globally distributed object.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
Abstract interface for local operators (e.g., matrices and preconditioners).
A parallel distribution of indices over processes.
One or more distributed dense vectors.
void reduce()
Sum values of a locally replicated multivector across all processes.
void scale(const Scalar &alpha)
Scale in place: this = alpha*this.
size_t getNumVectors() const
Number of columns in the multivector.
dual_view_type::t_dev::const_type getLocalViewDevice(Access::ReadOnlyStruct) const
Return a read-only, up-to-date view of this MultiVector's local data on device. This requires that th...
bool isConstantStride() const
Whether this multivector has constant stride between columns.
void putScalar(const Scalar &value)
Set all values in the multivector with the given value.
Abstract interface for operators (e.g., matrices and preconditioners).
Namespace Tpetra contains the class and methods constituting the Tpetra library.
void deep_copy(MultiVector< DS, DL, DG, DN > &dst, const MultiVector< SS, SL, SG, SN > &src)
Copy the contents of the MultiVector src into dst.
@ ADD_ASSIGN
Accumulate new values into existing values (may not be supported in all classes)
@ INSERT
Insert new values that don't currently exist.