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),
68 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A);
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>
83 isInitialized_ (false),
88 initializeTime_ (0.0),
92 if (! out_.is_null ()) {
93 *out_ <<
">>> DEBUG Ifpack2::LocalSparseTriangularSolver constructor" 100 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A);
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;
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 ()) {
138 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A_);
140 (A_crs.
is_null (), std::invalid_argument, prefix <<
141 "The input matrix A is not a Tpetra::CrsMatrix.");
144 auto G = A_->getGraph ();
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.");
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;
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().");
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.");
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_) {
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>
206 using Teuchos::rcpFromRef;
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;
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.");
235 (A_.is_null (), std::logic_error, prefix <<
"A_ is null. " 236 "Please report this bug to the Ifpack2 developers.");
238 (A_crs_.is_null (), std::logic_error, prefix <<
"A_crs_ is null. " 239 "Please report this bug to the Ifpack2 developers.");
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 ();
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>
303 const scalar_type& alpha,
304 const scalar_type& beta)
const 307 typedef scalar_type ST;
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>
423 auto comm = A_.is_null () ?
424 Tpetra::DefaultPlatform::getDefaultPlatform ().getComm () :
429 if (! comm.is_null () && comm->getRank () == 0) {
434 out <<
"\"Ifpack2::LocalSparseTriangularSolver\":" << endl;
444 template <
class MatrixType>
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>
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>
474 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::setMatrix: ";
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;
504 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A);
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
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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
static std::string name()