42 #ifndef TPETRA_CRSMATRIXSOLVEOP_HPP 43 #define TPETRA_CRSMATRIXSOLVEOP_HPP 50 #include <Tpetra_CrsMatrix.hpp> 88 template <
class Scalar,
89 class MatScalar = Scalar,
91 class GlobalOrdinal = Details::DefaultTypes::global_ordinal_type,
94 public Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
121 Teuchos::ETransp mode = Teuchos::NO_TRANS,
122 Scalar alpha = Teuchos::ScalarTraits<Scalar>::one (),
123 Scalar beta = Teuchos::ScalarTraits<Scalar>::zero ())
const 125 typedef Teuchos::ScalarTraits<Scalar> STOS;
126 const char prefix[] =
"Tpetra::CrsMatrixSolveOp::apply: ";
128 TEUCHOS_TEST_FOR_EXCEPTION
129 (!
matrix_->isFillComplete (), std::runtime_error,
130 prefix <<
"Underlying matrix is not fill complete.");
131 TEUCHOS_TEST_FOR_EXCEPTION
132 (!
matrix_->isLowerTriangular () && !
matrix_->isUpperTriangular (),
133 std::runtime_error, prefix <<
"The matrix is neither lower nor upper " 134 "triangular. Remember that in Tpetra, triangular-ness is a local " 135 "(per MPI process) property.");
136 TEUCHOS_TEST_FOR_EXCEPTION
138 prefix <<
"X and Y must have the same number of columns (vectors). " 141 TEUCHOS_TEST_FOR_EXCEPTION
142 (alpha != STOS::one () || beta != STOS::zero (), std::logic_error,
143 prefix <<
"The case alpha != 1 or beta != 0 has not yet been " 144 "implemented. Please speak with the Tpetra developers.");
145 if (mode == Teuchos::NO_TRANS) {
147 }
else if (mode == Teuchos::TRANS || mode == Teuchos::CONJ_TRANS) {
150 TEUCHOS_TEST_FOR_EXCEPTION
151 (
true, std::invalid_argument, prefix <<
"The 'mode' argument has an " 152 "invalid value " << mode <<
". Valid values are Teuchos::NO_TRANS=" 153 << Teuchos::NO_TRANS <<
", Teuchos::TRANS=" << Teuchos::TRANS <<
", " 154 "and Teuchos::CONJ_TRANS=" << Teuchos::CONJ_TRANS <<
".");
167 return matrix_->getRangeMap ();
173 return matrix_->getDomainMap ();
181 const Teuchos::RCP<const crs_matrix_type>
matrix_;
191 using Teuchos::NO_TRANS;
193 typedef Teuchos::ScalarTraits<Scalar> ST;
199 Teuchos::RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > importer =
200 matrix_->getGraph ()->getImporter ();
201 Teuchos::RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > exporter =
202 matrix_->getGraph ()->getExporter ();
203 Teuchos::RCP<const MV> X;
210 if (importer != null) {
211 if (importMV_ != null && importMV_->getNumVectors () != numVectors) {
214 if (importMV_ == null) {
215 importMV_ = Teuchos::rcp (
new MV (matrix_->getColMap (), numVectors));
218 if (exporter != null) {
219 if (exportMV_ != null && exportMV_->getNumVectors () != numVectors) {
222 if (exportMV_ == null) {
223 exportMV_ = Teuchos::rcp (
new MV (matrix_->getRowMap (), numVectors));
237 if (exporter != null) {
238 exportMV_->doImport (X_in, *exporter,
INSERT);
244 X = Teuchos::rcp (
new MV (X_in));
248 X = Teuchos::rcpFromRef (X_in);
254 if (importer != null) {
255 matrix_->template localSolve<Scalar, Scalar> (*X, *
importMV_, NO_TRANS);
266 matrix_->template localSolve<Scalar, Scalar> (*X, Y, NO_TRANS);
270 matrix_->template localSolve<Scalar, Scalar> (*X, Y_in, NO_TRANS);
276 void applyTranspose (
const MV& X_in, MV& Y_in,
const Teuchos::ETransp mode)
const 278 typedef Teuchos::ScalarTraits<Scalar> ST;
281 TEUCHOS_TEST_FOR_EXCEPTION
282 (mode != Teuchos::TRANS && mode != Teuchos::CONJ_TRANS, std::logic_error,
283 "Tpetra::CrsMatrixSolveOp::applyTranspose: mode is neither TRANS nor " 284 "CONJ_TRANS. Should never get here! Please report this bug to the " 285 "Tpetra developers.");
288 Teuchos::RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > importer =
289 matrix_->getGraph ()->getImporter ();
290 Teuchos::RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > exporter =
291 matrix_->getGraph ()->getExporter ();
292 Teuchos::RCP<const MV> X;
299 if (importer != null) {
300 if (importMV_ != null && importMV_->getNumVectors() != numVectors) {
303 if (importMV_ == null) {
304 importMV_ = Teuchos::rcp(
new MV(matrix_->getColMap(),numVectors) );
307 if (exporter != null) {
308 if (exportMV_ != null && exportMV_->getNumVectors() != numVectors) {
311 if (exportMV_ == null) {
312 exportMV_ = Teuchos::rcp(
new MV(matrix_->getRowMap(),numVectors) );
326 if (importer != null) {
327 importMV_->doImport(X_in,*importer,
INSERT);
333 X = Teuchos::rcp(
new MV(X_in));
337 X = Teuchos::rcpFromRef (X_in);
344 if (exporter != null) {
345 matrix_->template localSolve<Scalar, Scalar> (*X, *
exportMV_,
346 Teuchos::CONJ_TRANS);
356 matrix_->template localSolve<Scalar, Scalar> (*X, Y, Teuchos::CONJ_TRANS);
360 matrix_->template localSolve<Scalar, Scalar> (*X, Y_in, Teuchos::CONJ_TRANS);
373 template<
class OpScalar,
378 Teuchos::RCP<CrsMatrixSolveOp<OpScalar, MatScalar, LocalOrdinal, GlobalOrdinal, Node> >
386 #endif // TPETRA_CRSMATRIXSOLVEOP_HPP Namespace Tpetra contains the class and methods constituting the Tpetra library.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
KokkosClassic::DefaultNode::DefaultNodeType node_type
Default value of Node template parameter.
const Teuchos::RCP< const crs_matrix_type > matrix_
The underlying CrsMatrix.
One or more distributed dense vectors.
void deep_copy(MultiVector< DS, DL, DG, DN, dstClassic > &dst, const MultiVector< SS, SL, SG, SN, srcClassic > &src)
Copy the contents of the MultiVector src into dst.
CrsMatrix< MatScalar, LocalOrdinal, GlobalOrdinal, Node > crs_matrix_type
The specialization of CrsMatrix which this class wraps.
bool isConstantStride() const
Whether this multivector has constant stride between columns.
int local_ordinal_type
Default value of LocalOrdinal template parameter.
Teuchos::RCP< CrsMatrixSolveOp< OpScalar, MatScalar, LocalOrdinal, GlobalOrdinal, Node > > createCrsMatrixSolveOp(const Teuchos::RCP< const CrsMatrix< MatScalar, LocalOrdinal, GlobalOrdinal, Node > > &A)
Nonmember function that wraps a CrsMatrix in a CrsMatrixSolveOp.
Map< LocalOrdinal, GlobalOrdinal, Node > map_type
The specialization of Map which this class uses.
Teuchos::RCP< const map_type > getRangeMap() const
The range Map of this operator. This is the domain Map of the underlying CrsMatrix.
Insert new values that don't currently exist.
void applyNonTranspose(const MV &X_in, MV &Y_in) const
Do the non-transpose solve.
Abstract interface for operators (e.g., matrices and preconditioners).
Sum new values into existing values.
Teuchos::RCP< const map_type > getDomainMap() const
The domain Map of this operator. This is the range map of the underlying CrsMatrix.
Teuchos::RCP< MV > importMV_
Cached temporary destination of Import operation in apply().
void putScalar(const Scalar &value)
Set all values in the multivector with the given value.
bool hasTransposeApply() const
Whether apply() can solve with the (conjugate) transpose of the matrix.
Describes a parallel distribution of objects over processes.
Wrap a CrsMatrix instance's triangular solve in an Operator.
void applyTranspose(const MV &X_in, MV &Y_in, const Teuchos::ETransp mode) const
Do the transpose or conjugate transpose solve.
size_t getNumVectors() const
Number of columns in the multivector.
virtual ~CrsMatrixSolveOp()
Destructor (virtual for memory safety of derived classes).
void doExport(const SrcDistObject &source, const Export< LocalOrdinal, GlobalOrdinal, Node > &exporter, CombineMode CM)
Export data into this object using an Export object ("forward mode").
CrsMatrixSolveOp(const Teuchos::RCP< const crs_matrix_type > &A)
Constructor; takes a CrsMatrix to use for local triangular solves.
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
Compute , where represents the result of the local triangular solve.
Teuchos::RCP< MV > exportMV_
Cached temporary source of Export operation in apply().