10 #ifndef TPETRA_CRSMATRIXMULTIPLYOP_HPP
11 #define TPETRA_CRSMATRIXMULTIPLYOP_HPP
19 #include "Tpetra_CrsMatrix.hpp"
23 #include "Tpetra_LocalCrsMatrixOperator.hpp"
63 template <
class Scalar,
69 public Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>
79 using local_matrix_device_type =
93 A->getLocalMatrixDevice ()))
108 Teuchos::ETransp mode = Teuchos::NO_TRANS,
109 Scalar alpha = Teuchos::ScalarTraits<Scalar>::one (),
110 Scalar beta = Teuchos::ScalarTraits<Scalar>::zero ())
const override
112 TEUCHOS_TEST_FOR_EXCEPTION
113 (!
matrix_->isFillComplete (), std::runtime_error,
114 Teuchos::typeName (*
this) <<
"::apply(): underlying matrix is not fill-complete.");
115 TEUCHOS_TEST_FOR_EXCEPTION
117 Teuchos::typeName (*
this) <<
"::apply(X,Y): X and Y must have the same number of vectors.");
118 TEUCHOS_TEST_FOR_EXCEPTION
119 (Teuchos::ScalarTraits<Scalar>::isComplex && mode == Teuchos::TRANS, std::logic_error,
120 Teuchos::typeName (*
this) <<
"::apply() does not currently support transposed multiplications for complex scalar types.");
121 if (mode == Teuchos::NO_TRANS) {
140 return matrix_->getDomainMap ();
145 return matrix_->getRangeMap ();
154 const Teuchos::RCP<const crs_matrix_type>
matrix_;
172 mutable Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
importMV_;
186 mutable Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
exportMV_;
193 Teuchos::ETransp mode,
202 using STS = Teuchos::ScalarTraits<Scalar>;
209 RCP<const import_type> importer =
matrix_->getGraph()->getImporter();
210 RCP<const export_type> exporter =
matrix_->getGraph()->getExporter();
214 const bool Y_is_overwritten = (beta == STS::zero ());
215 const int myRank =
matrix_->getComm ()->getRank ();
216 if (Y_is_replicated && myRank != 0) {
225 X = Teuchos::rcp (
new MV (X_in, Teuchos::Copy));
229 X = Teuchos::rcpFromRef (X_in);
233 if (importer != null) {
241 if (exporter != null) {
251 if (exporter != null) {
256 auto X_lcl = X->getLocalViewDevice(Access::ReadOnly);
260 if (importer != null) {
262 auto Y_lcl =
importMV_->getLocalViewDevice(Access::OverwriteAll);
265 if (Y_is_overwritten) {
278 MV Y (Y_in, Teuchos::Copy);
279 nonconst_view_type Y_lcl;
287 nonconst_view_type Y_lcl;
295 if (Y_is_replicated) {
308 using Teuchos::NO_TRANS;
311 using Teuchos::rcp_const_cast;
312 using Teuchos::rcpFromRef;
315 typedef Teuchos::ScalarTraits<Scalar> STS;
318 if (alpha == STS::zero ()) {
319 if (beta == STS::zero ()) {
321 }
else if (beta != STS::one ()) {
336 RCP<const import_type> importer =
matrix_->getGraph()->getImporter();
337 RCP<const export_type> exporter =
matrix_->getGraph()->getExporter();
343 const bool Y_is_overwritten = (beta == STS::zero());
346 const bool Y_is_replicated =
356 if (Y_is_replicated &&
matrix_->getComm ()->getRank () > 0) {
363 RCP<const MV> X_colMap;
364 if (importer.is_null ()) {
372 RCP<MV> X_colMapNonConst = getColumnMapMultiVector (X_in,
true);
374 X_colMap = rcp_const_cast<
const MV> (X_colMapNonConst);
379 X_colMap = rcpFromRef (X_in);
383 ProfilingRegion regionImport (
"Tpetra::CrsMatrixMultiplyOp::apply: Import");
389 RCP<MV> X_colMapNonConst = getColumnMapMultiVector (X_in);
392 X_colMapNonConst->doImport (X_in, *importer,
INSERT);
393 X_colMap = rcp_const_cast<
const MV> (X_colMapNonConst);
400 RCP<MV> Y_rowMap = getRowMapMultiVector (Y_in);
402 auto X_lcl = X_colMap->getLocalViewDevice(Access::ReadOnly);
409 if (! exporter.is_null ()) {
410 auto Y_lcl = Y_rowMap->getLocalViewDevice(Access::OverwriteAll);
413 alpha, STS::zero ());
415 ProfilingRegion regionExport (
"Tpetra::CrsMatrixMultiplyOp::apply: Export");
421 if (Y_is_overwritten) {
448 Y_rowMap = getRowMapMultiVector (Y_in,
true);
452 if (beta != STS::zero ()) {
455 nonconst_view_type Y_lcl;
456 if(Y_is_overwritten) Y_lcl = Y_rowMap->getLocalViewDevice(Access::OverwriteAll);
457 else Y_lcl = Y_rowMap->getLocalViewDevice(Access::ReadWrite);
463 nonconst_view_type Y_lcl;
475 if (Y_is_replicated) {
476 ProfilingRegion regionReduce (
"Tpetra::CrsMatrixMultiplyOp::apply: Reduce Y");
503 const bool force =
false)
const
511 RCP<const import_type> importer =
matrix_->getGraph ()->getImporter ();
512 RCP<const map_type> colMap =
matrix_->getColMap ();
525 if (! importer.is_null () || force) {
527 X_colMap = rcp (
new MV (colMap, numVecs));
566 getRowMapMultiVector (
const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y_rangeMap,
567 const bool force =
false)
const
572 typedef Export<LocalOrdinal,GlobalOrdinal,Node> export_type;
574 const size_t numVecs = Y_rangeMap.getNumVectors ();
575 RCP<const export_type> exporter =
matrix_->getGraph ()->getExporter ();
576 RCP<const map_type> rowMap =
matrix_->getRowMap ();
588 if (! exporter.is_null () || force) {
590 Y_rowMap = rcp (
new MV (rowMap, numVecs));
608 template <
class OpScalar,
614 CrsMatrixMultiplyOp<OpScalar, MatScalar, LocalOrdinal, GlobalOrdinal, Node> >
619 GlobalOrdinal, Node> op_type;
620 return Teuchos::rcp (
new op_type (A));
625 #endif // TPETRA_CRSMATRIXMULTIPLYOP_HPP
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).
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
Teuchos::RCP< const map_type > getRangeMap() const override
The range Map of this Operator.
Declaration of Tpetra::Details::Profiling, a scope guard for Kokkos Profiling.
void putScalar(const Scalar &value)
Set all values in the multivector with the given value.
size_t getNumVectors() const
Number of columns in the multivector.
Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > exportMV_
Row Map MultiVector used in apply().
bool isConstantStride() const
Whether this multivector has constant stride between columns.
bool isDistributed() const
Whether this is a globally distributed object.
Teuchos::RCP< const map_type > getDomainMap() const override
The domain Map of this Operator.
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 .
bool hasTransposeApply() const override
Whether this Operator's apply() method can apply the transpose or conjugate transpose.
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...
Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > importMV_
Column Map MultiVector used in apply().
~CrsMatrixMultiplyOp() override=default
Destructor (virtual for memory safety of derived classes).
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.
CrsMatrixMultiplyOp(const Teuchos::RCP< const crs_matrix_type > &A)
Constructor.
LocalCrsMatrixOperator< Scalar, MatScalar, typename crs_matrix_type::device_type > localMultiply_
Implementation of local sparse matrix-vector multiply.
Insert new values that don't currently exist.
void scale(const Scalar &alpha)
Scale in place: this = alpha*this.
Abstract interface for operators (e.g., matrices and preconditioners).
const Teuchos::RCP< const crs_matrix_type > matrix_
The underlying CrsMatrix object.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
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. ...
typename Node::device_type device_type
The Kokkos device type.
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...
Forward declaration of Tpetra::CrsMatrixMultiplyOp.
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.
A parallel distribution of indices over processes.
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").
A class for wrapping a CrsMatrix multiply in a Operator.
Stand-alone utility functions and macros.
void reduce()
Sum values of a locally replicated multivector across all processes.
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.
Accumulate new values into existing values (may not be supported in all classes)
Declaration of Tpetra::Details::Behavior, a class that describes Tpetra's behavior.