43 #ifndef IFPACK2_LOCALSPARSETRIANGULARSOLVER_DEF_HPP 44 #define IFPACK2_LOCALSPARSETRIANGULARSOLVER_DEF_HPP 46 #include "Tpetra_CrsMatrix.hpp" 47 #include "Tpetra_DefaultPlatform.hpp" 51 template<
class MatrixType>
55 isInitialized_ (false),
60 initializeTime_ (0.0),
67 Teuchos::RCP<const crs_matrix_type> A_crs =
68 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A);
69 TEUCHOS_TEST_FOR_EXCEPTION
70 (A_crs.is_null (), std::invalid_argument,
71 "Ifpack2::LocalSparseTriangularSolver constructor: " 72 "The input matrix A is not a Tpetra::CrsMatrix.");
77 template<
class MatrixType>
80 const Teuchos::RCP<Teuchos::FancyOStream>& out) :
83 isInitialized_ (false),
88 initializeTime_ (0.0),
92 if (! out_.is_null ()) {
93 *out_ <<
">>> DEBUG Ifpack2::LocalSparseTriangularSolver constructor" 99 Teuchos::RCP<const crs_matrix_type> A_crs =
100 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A);
101 TEUCHOS_TEST_FOR_EXCEPTION
102 (A_crs.is_null (), std::invalid_argument,
103 "Ifpack2::LocalSparseTriangularSolver constructor: " 104 "The input matrix A is not a Tpetra::CrsMatrix.");
109 template<
class MatrixType>
114 template<
class MatrixType>
120 template<
class MatrixType>
125 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::initialize: ";
126 if (! out_.is_null ()) {
127 *out_ <<
">>> DEBUG " << prefix << std::endl;
130 TEUCHOS_TEST_FOR_EXCEPTION
131 (A_.is_null (), std::runtime_error, prefix <<
"You must call " 132 "setMatrix() with a nonnull input matrix before you may call " 133 "initialize() or compute().");
134 if (A_crs_.is_null ()) {
137 Teuchos::RCP<const crs_matrix_type> A_crs =
138 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A_);
139 TEUCHOS_TEST_FOR_EXCEPTION
140 (A_crs.is_null (), std::invalid_argument, prefix <<
141 "The input matrix A is not a Tpetra::CrsMatrix.");
144 auto G = A_->getGraph ();
145 TEUCHOS_TEST_FOR_EXCEPTION
146 (G.is_null (), std::logic_error, prefix <<
"A_ and A_crs_ are nonnull, " 147 "but A_'s RowGraph G is null. " 148 "Please report this bug to the Ifpack2 developers.");
152 TEUCHOS_TEST_FOR_EXCEPTION
153 (! G->isFillComplete (), std::runtime_error,
"If you call this method, " 154 "the matrix's graph must be fill complete. It is not.");
156 isInitialized_ =
true;
160 template<
class MatrixType>
165 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::compute: ";
166 if (! out_.is_null ()) {
167 *out_ <<
">>> DEBUG " << prefix << std::endl;
170 TEUCHOS_TEST_FOR_EXCEPTION
171 (A_.is_null (), std::runtime_error, prefix <<
"You must call " 172 "setMatrix() with a nonnull input matrix before you may call " 173 "initialize() or compute().");
174 TEUCHOS_TEST_FOR_EXCEPTION
175 (A_crs_.is_null (), std::logic_error, prefix <<
"A_ is nonnull, but " 176 "A_crs_ is null. Please report this bug to the Ifpack2 developers.");
178 TEUCHOS_TEST_FOR_EXCEPTION
179 (! A_crs_->isFillComplete (), std::runtime_error,
"If you call this " 180 "method, the matrix must be fill complete. It is not.");
182 if (! isInitialized_) {
185 TEUCHOS_TEST_FOR_EXCEPTION
186 (! isInitialized_, std::logic_error, prefix <<
"initialize() should have " 187 "been called by this point, but isInitialized_ is false. " 188 "Please report this bug to the Ifpack2 developers.");
194 template<
class MatrixType>
200 Teuchos::ETransp mode,
206 using Teuchos::rcpFromRef;
208 typedef Teuchos::ScalarTraits<ST> STS;
209 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::apply: ";
210 if (! out_.is_null ()) {
211 *out_ <<
">>> DEBUG " << prefix;
212 if (A_crs_.is_null ()) {
213 *out_ <<
"A_crs_ is null!" << std::endl;
216 const std::string uplo = A_crs_->isUpperTriangular () ?
"U" :
217 (A_crs_->isLowerTriangular () ?
"L" :
"N");
218 const std::string trans = (mode == Teuchos::CONJ_TRANS) ?
"C" :
219 (mode == Teuchos::TRANS ?
"T" :
"N");
220 const std::string diag =
221 (A_crs_->getNodeNumDiags () < A_crs_->getNodeNumRows ()) ?
"U" :
"N";
222 *out_ <<
"uplo=\"" << uplo
223 <<
"\", trans=\"" << trans
224 <<
"\", diag=\"" << diag <<
"\"" << std::endl;
228 TEUCHOS_TEST_FOR_EXCEPTION
229 (! isComputed (), std::runtime_error, prefix <<
"If compute() has not yet " 230 "been called, or if you have changed the matrix via setMatrix(), you must " 231 "call compute() before you may call this method.");
234 TEUCHOS_TEST_FOR_EXCEPTION
235 (A_.is_null (), std::logic_error, prefix <<
"A_ is null. " 236 "Please report this bug to the Ifpack2 developers.");
237 TEUCHOS_TEST_FOR_EXCEPTION
238 (A_crs_.is_null (), std::logic_error, prefix <<
"A_crs_ is null. " 239 "Please report this bug to the Ifpack2 developers.");
242 TEUCHOS_TEST_FOR_EXCEPTION
243 (! A_crs_->isFillComplete (), std::runtime_error,
"If you call this " 244 "method, the matrix must be fill complete. It is not. This means that " 245 " you must have called resumeFill() on the matrix before calling apply(). " 246 "This is NOT allowed. Note that this class may use the matrix's data in " 247 "place without copying it. Thus, you cannot change the matrix and expect " 248 "the solver to stay the same. If you have changed the matrix, first call " 249 "fillComplete() on it, then call compute() on this object, before you call" 250 " apply(). You do NOT need to call setMatrix, as long as the matrix " 251 "itself (that is, its address in memory) is the same.");
253 auto G = A_->getGraph ();
254 TEUCHOS_TEST_FOR_EXCEPTION
255 (G.is_null (), std::logic_error, prefix <<
"A_ and A_crs_ are nonnull, " 256 "but A_'s RowGraph G is null. " 257 "Please report this bug to the Ifpack2 developers.");
258 auto importer = G->getImporter ();
259 auto exporter = G->getExporter ();
261 if (! importer.is_null ()) {
262 if (X_colMap_.is_null () || X_colMap_->getNumVectors () != X.getNumVectors ()) {
263 X_colMap_ = rcp (
new MV (importer->getTargetMap (), X.getNumVectors ()));
266 X_colMap_->putScalar (STS::zero ());
271 X_colMap_->doImport (X, *importer, Tpetra::ZERO);
273 RCP<const MV> X_cur = importer.is_null () ? rcpFromRef (X) :
274 Teuchos::rcp_const_cast<
const MV> (X_colMap_);
276 if (! exporter.is_null ()) {
277 if (Y_rowMap_.is_null () || Y_rowMap_->getNumVectors () != Y.getNumVectors ()) {
278 Y_rowMap_ = rcp (
new MV (exporter->getSourceMap (), Y.getNumVectors ()));
281 Y_rowMap_->putScalar (STS::zero ());
283 Y_rowMap_->doExport (Y, *importer, Tpetra::ADD);
285 RCP<MV> Y_cur = exporter.is_null () ? rcpFromRef (Y) : Y_rowMap_;
287 localApply (*X_cur, *Y_cur, mode, alpha, beta);
289 if (! exporter.is_null ()) {
290 Y.putScalar (STS::zero ());
291 Y.doExport (*Y_cur, *exporter, Tpetra::ADD);
297 template<
class MatrixType>
302 const Teuchos::ETransp mode,
303 const scalar_type& alpha,
304 const scalar_type& beta)
const 307 typedef scalar_type ST;
308 typedef Teuchos::ScalarTraits<ST> STS;
310 if (beta == STS::zero ()) {
311 if (alpha == STS::zero ()) {
312 Y.putScalar (STS::zero ());
315 A_crs_->template localSolve<ST, ST> (X, Y, mode);
316 if (alpha != STS::one ()) {
322 if (alpha == STS::zero ()) {
326 MV Y_tmp (Y, Teuchos::Copy);
327 A_crs_->template localSolve<ST, ST> (X, Y_tmp, mode);
328 Y.update (alpha, Y_tmp, beta);
334 template <
class MatrixType>
338 return numInitialize_;
341 template <
class MatrixType>
348 template <
class MatrixType>
355 template <
class MatrixType>
359 return initializeTime_;
362 template<
class MatrixType>
369 template<
class MatrixType>
376 template <
class MatrixType>
381 std::ostringstream os;
386 os <<
"\"Ifpack2::LocalSparseTriangularSolver\": {";
387 if (this->getObjectLabel () !=
"") {
388 os <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
390 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", " 391 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
394 os <<
"Matrix: null";
397 os <<
"Matrix: not null" 398 <<
", Global matrix dimensions: [" 399 << A_->getGlobalNumRows () <<
", " 400 << A_->getGlobalNumCols () <<
"]";
407 template <
class MatrixType>
410 const Teuchos::EVerbosityLevel verbLevel)
const 415 const Teuchos::EVerbosityLevel vl
416 = (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
418 if (vl != Teuchos::VERB_NONE) {
423 auto comm = A_.is_null () ?
424 Tpetra::DefaultPlatform::getDefaultPlatform ().getComm () :
429 if (! comm.is_null () && comm->getRank () == 0) {
431 Teuchos::OSTab tab0 (out);
434 out <<
"\"Ifpack2::LocalSparseTriangularSolver\":" << endl;
435 Teuchos::OSTab tab1 (out);
436 out <<
"Scalar: " << Teuchos::TypeNameTraits<scalar_type>::name () << endl
437 <<
"LocalOrdinal: " << Teuchos::TypeNameTraits<local_ordinal_type>::name () << endl
438 <<
"GlobalOrdinal: " << Teuchos::TypeNameTraits<global_ordinal_type>::name () << endl
439 <<
"Node: " << Teuchos::TypeNameTraits<node_type>::name () << endl;
444 template <
class MatrixType>
445 Teuchos::RCP<const typename LocalSparseTriangularSolver<MatrixType>::map_type>
449 TEUCHOS_TEST_FOR_EXCEPTION
450 (A_.is_null (), std::runtime_error,
451 "Ifpack2::LocalSparseTriangularSolver::getDomainMap: " 452 "The matrix is null. Please call setMatrix() with a nonnull input " 453 "before calling this method.");
454 return A_->getDomainMap ();
457 template <
class MatrixType>
458 Teuchos::RCP<const typename LocalSparseTriangularSolver<MatrixType>::map_type>
462 TEUCHOS_TEST_FOR_EXCEPTION
463 (A_.is_null (), std::runtime_error,
464 "Ifpack2::LocalSparseTriangularSolver::getRangeMap: " 465 "The matrix is null. Please call setMatrix() with a nonnull input " 466 "before calling this method.");
467 return A_->getRangeMap ();
470 template<
class MatrixType>
472 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
474 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::setMatrix: ";
480 if (A.getRawPtr () != A_.getRawPtr ()) {
482 TEUCHOS_TEST_FOR_EXCEPTION
483 (! A.is_null () && A->getComm ()->getSize () == 1 &&
484 A->getNodeNumRows () != A->getNodeNumCols (),
485 std::runtime_error, prefix <<
"If A's communicator only contains one " 486 "process, then A must be square. Instead, you provided a matrix A with " 487 << A->getNodeNumRows () <<
" rows and " << A->getNodeNumCols ()
493 isInitialized_ =
false;
499 A_crs_ = Teuchos::null;
503 Teuchos::RCP<const crs_matrix_type> A_crs =
504 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A);
505 TEUCHOS_TEST_FOR_EXCEPTION
506 (A_crs.is_null (), std::invalid_argument, prefix <<
507 "The input matrix A is not a Tpetra::CrsMatrix.");
516 #define IFPACK2_LOCALSPARSETRIANGULARSOLVER_INSTANT(S,LO,GO,N) \ 517 template class Ifpack2::LocalSparseTriangularSolver< Tpetra::RowMatrix<S, LO, GO, N> >; 519 #endif // IFPACK2_LOCALSPARSETRIANGULARSOLVER_DEF_HPP void apply(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Apply the preconditioner to X, and put the result in Y.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:196
int getNumInitialize() const
Return the number of calls to initialize().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:337
Teuchos::RCP< const map_type > getDomainMap() const
The domain of this operator.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:447
Teuchos::RCP< const map_type > getRangeMap() const
The range of this operator.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:460
int getNumCompute() const
Return the number of calls to compute().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:344
double getApplyTime() const
Return the time spent in apply().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:372
void compute()
"Numeric" phase of setup
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:163
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print this object with given verbosity to the given output stream.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:409
MatrixType::global_ordinal_type global_ordinal_type
Type of the global indices of the input matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:99
MatrixType::node_type node_type
Node type of the input matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:101
double getComputeTime() const
Return the time spent in compute().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:365
int getNumApply() const
Return the number of calls to apply().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:351
"Preconditioner" that solves local sparse triangular systems.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:83
void initialize()
"Symbolic" phase of setup
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:123
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set this preconditioner's matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:472
LocalSparseTriangularSolver(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:53
double getInitializeTime() const
Return the time spent in initialize().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:358
MatrixType::local_ordinal_type local_ordinal_type
Type of the local indices of the input matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:97
void setParameters(const Teuchos::ParameterList ¶ms)
Set this object's parameters.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:117
MatrixType::scalar_type scalar_type
Type of the entries of the input matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:95
std::string description() const
A one-line description of this object.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:379
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
virtual ~LocalSparseTriangularSolver()
Destructor (virtual for memory safety).
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:111