43 #ifndef IFPACK2_SPARSECONTAINER_DEF_HPP 44 #define IFPACK2_SPARSECONTAINER_DEF_HPP 49 #include "Teuchos_DefaultMpiComm.hpp" 51 #include "Teuchos_DefaultSerialComm.hpp" 57 template<
class MatrixType,
class InverseType>
60 const Teuchos::ArrayView<const local_ordinal_type>& localRows) :
61 Container<MatrixType> (matrix, localRows),
62 numRows_ (localRows.size ()),
63 IsInitialized_ (false),
66 localComm_ (
Teuchos::rcp (new
Teuchos::MpiComm<int> (MPI_COMM_SELF)))
73 template<
class MatrixType,
class InverseType>
78 template<
class MatrixType,
class InverseType>
89 template<
class MatrixType,
class InverseType>
92 return IsInitialized_;
96 template<
class MatrixType,
class InverseType>
103 template<
class MatrixType,
class InverseType>
110 template<
class MatrixType,
class InverseType>
115 typedef Tpetra::Map<InverseLocalOrdinal, InverseGlobalOrdinal,
116 InverseNode> map_type;
117 typedef Tpetra::CrsMatrix<InverseScalar, InverseLocalOrdinal,
118 InverseGlobalOrdinal, InverseNode> crs_matrix_type;
124 IsInitialized_ =
false;
131 localMap_ = rcp (
new map_type (numRows_, 0, localComm_));
132 diagBlock_ = rcp (
new crs_matrix_type (localMap_, 0));
137 Inverse_ = rcp (
new InverseType (diagBlock_));
138 Inverse_->setParameters (List_);
140 IsInitialized_ =
true;
144 template<
class MatrixType,
class InverseType>
158 Inverse_->initialize ();
159 Inverse_->compute ();
165 template<
class MatrixType,
class InverseType>
167 applyImpl (
const Tpetra::MultiVector<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode>& X,
168 Tpetra::MultiVector<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode>& Y,
169 Teuchos::ETransp mode,
171 InverseScalar beta)
const 173 TEUCHOS_TEST_FOR_EXCEPTION(
174 Inverse_->getDomainMap ()->getNodeNumElements () != X.getLocalLength (),
175 std::logic_error,
"Ifpack2::SparseContainer::apply: Inverse_ " 176 "operator and X have incompatible dimensions (" <<
177 Inverse_->getDomainMap ()->getNodeNumElements () <<
" resp. " 178 << X.getLocalLength () <<
"). Please report this bug to " 179 "the Ifpack2 developers.");
180 TEUCHOS_TEST_FOR_EXCEPTION(
181 Inverse_->getRangeMap ()->getNodeNumElements () != Y.getLocalLength (),
182 std::logic_error,
"Ifpack2::SparseContainer::apply: Inverse_ " 183 "operator and Y have incompatible dimensions (" <<
184 Inverse_->getRangeMap ()->getNodeNumElements () <<
" resp. " 185 << Y.getLocalLength () <<
"). Please report this bug to " 186 "the Ifpack2 developers.");
188 Inverse_->apply (X, Y, mode, alpha, beta);
192 template<
class MatrixType,
class InverseType>
194 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
195 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
196 Teuchos::ETransp mode,
198 scalar_type beta)
const 200 using Teuchos::ArrayView;
215 typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV_mat;
217 typedef Tpetra::MultiVector<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode> MV_inv;
219 const size_t numVecs = X.getNumVectors ();
221 TEUCHOS_TEST_FOR_EXCEPTION(
222 ! IsComputed_, std::runtime_error,
"Ifpack2::SparseContainer::apply: " 223 "You must have called the compute() method before you may call apply(). " 224 "You may call the apply() method as many times as you want after calling " 225 "compute() once, but you must have called compute() at least once.");
226 TEUCHOS_TEST_FOR_EXCEPTION(
227 numVecs != Y.getNumVectors (), std::runtime_error,
228 "Ifpack2::SparseContainer::apply: X and Y have different numbers of " 229 "vectors. X has " << X.getNumVectors ()
230 <<
", but Y has " << X.getNumVectors () <<
".");
262 X_ = rcp (
new MV_inv (Inverse_->getDomainMap (), numVecs));
264 RCP<MV_inv> X_local = X_;
265 TEUCHOS_TEST_FOR_EXCEPTION(
266 X_local->getLocalLength () != numRows_, std::logic_error,
267 "Ifpack2::SparseContainer::apply: " 268 "X_local has length " << X_local->getLocalLength () <<
", which does " 269 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 270 "the Ifpack2 developers.");
271 ArrayView<const local_ordinal_type> localRows = this->
getLocalRows ();
272 mvgs.gather (*X_local, X, localRows);
280 Y_ = rcp (
new MV_inv (Inverse_->getRangeMap (), numVecs));
282 RCP<MV_inv> Y_local = Y_;
283 TEUCHOS_TEST_FOR_EXCEPTION(
284 Y_local->getLocalLength () != numRows_, std::logic_error,
285 "Ifpack2::SparseContainer::apply: " 286 "Y_local has length " << X_local->getLocalLength () <<
", which does " 287 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 288 "the Ifpack2 developers.");
289 mvgs.gather (*Y_local, Y, localRows);
293 this->applyImpl (*X_local, *Y_local, mode, as<InverseScalar> (alpha),
294 as<InverseScalar> (beta));
298 mvgs.scatter (Y, *Y_local, localRows);
302 template<
class MatrixType,
class InverseType>
304 weightedApply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
305 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
306 const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& D,
307 Teuchos::ETransp mode,
309 scalar_type beta)
const 311 using Teuchos::ArrayRCP;
312 using Teuchos::ArrayView;
313 using Teuchos::Range1D;
316 using Teuchos::rcp_const_cast;
319 typedef Teuchos::ScalarTraits<scalar_type> STS;
334 typedef Tpetra::MultiVector<InverseScalar, InverseLocalOrdinal,
335 InverseGlobalOrdinal, InverseNode> MV_inv;
337 typedef Tpetra::Vector<InverseScalar, InverseLocalOrdinal,
338 InverseGlobalOrdinal, InverseNode> V_inv;
340 const size_t numVecs = X.getNumVectors ();
342 TEUCHOS_TEST_FOR_EXCEPTION(
343 ! IsComputed_, std::runtime_error,
"Ifpack2::SparseContainer::" 344 "weightedApply: You must have called the compute() method before you may " 345 "call apply(). You may call the apply() method as many times as you want " 346 "after calling compute() once, but you must have called compute() at least " 348 TEUCHOS_TEST_FOR_EXCEPTION(
349 Inverse_.is_null (), std::logic_error,
"Ifpack2::SparseContainer::" 350 "weightedApply: Inverse_ is null. Please report this bug to the Ifpack2 " 352 TEUCHOS_TEST_FOR_EXCEPTION(
353 numVecs != Y.getNumVectors (), std::runtime_error,
354 "Ifpack2::SparseContainer::weightedApply: X and Y have different numbers " 355 "of vectors. X has " << X.getNumVectors () <<
", but Y has " 356 << X.getNumVectors () <<
".");
388 X_ = rcp (
new MV_inv (Inverse_->getDomainMap (), numVecs));
390 RCP<MV_inv> X_local = X_;
391 TEUCHOS_TEST_FOR_EXCEPTION(
392 X_local->getLocalLength () != numRows_, std::logic_error,
393 "Ifpack2::SparseContainer::weightedApply: " 394 "X_local has length " << X_local->getLocalLength () <<
", which does " 395 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 396 "the Ifpack2 developers.");
397 ArrayView<const local_ordinal_type> localRows = this->
getLocalRows ();
398 mvgs.gather (*X_local, X, localRows);
406 Y_ = rcp (
new MV_inv (Inverse_->getRangeMap (), numVecs));
408 RCP<MV_inv> Y_local = Y_;
409 TEUCHOS_TEST_FOR_EXCEPTION(
410 Y_local->getLocalLength () != numRows_, std::logic_error,
411 "Ifpack2::SparseContainer::weightedApply: " 412 "Y_local has length " << X_local->getLocalLength () <<
", which does " 413 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 414 "the Ifpack2 developers.");
415 mvgs.gather (*Y_local, Y, localRows);
427 RCP<V_inv> D_local = rcp (
new V_inv (Inverse_->getDomainMap ()));
428 TEUCHOS_TEST_FOR_EXCEPTION(
429 D_local->getLocalLength () != numRows_, std::logic_error,
430 "Ifpack2::SparseContainer::weightedApply: " 431 "D_local has length " << X_local->getLocalLength () <<
", which does " 432 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 433 "the Ifpack2 developers.");
434 mvgs.gather (*D_local, D, localRows);
435 RCP<MV_inv> X_scaled = rcp (
new MV_inv (Inverse_->getDomainMap (), numVecs));
436 X_scaled->elementWiseMultiply (STS::one (), *D_local, *X_local, STS::zero ());
444 if (beta == STS::zero ()) {
447 Y_temp = rcp (
new MV_inv (Inverse_->getRangeMap (), numVecs));
451 Inverse_->apply (*X_scaled, *Y_temp, mode);
457 Y_local->elementWiseMultiply (alpha, *D_local, *Y_temp, beta);
461 mvgs.scatter (Y, *Y_local, localRows);
465 template<
class MatrixType,
class InverseType>
468 Teuchos::FancyOStream fos(Teuchos::rcp(&os,
false));
469 fos.setOutputToRootOnly(0);
475 template<
class MatrixType,
class InverseType>
478 std::ostringstream oss;
479 oss << Teuchos::Describable::description();
482 oss <<
"{status = initialized, computed";
485 oss <<
"{status = initialized, not computed";
489 oss <<
"{status = not initialized, not computed";
497 template<
class MatrixType,
class InverseType>
501 if(verbLevel==Teuchos::VERB_NONE)
return;
502 os <<
"================================================================================" << endl;
503 os <<
"Ifpack2::SparseContainer" << endl;
504 os <<
"Number of rows = " << numRows_ << endl;
505 os <<
"isInitialized() = " << IsInitialized_ << endl;
506 os <<
"isComputed() = " << IsComputed_ << endl;
507 os <<
"================================================================================" << endl;
512 template<
class MatrixType,
class InverseType>
514 extract (
const Teuchos::RCP<const row_matrix_type>& globalMatrix)
516 using Teuchos::Array;
517 using Teuchos::ArrayView;
519 const size_t MatrixInNumRows = globalMatrix->getNodeNumRows ();
522 ArrayView<const local_ordinal_type> localRows = this->
getLocalRows ();
523 for (
size_t j = 0; j < numRows_; ++j) {
524 TEUCHOS_TEST_FOR_EXCEPTION(
526 static_cast<size_t> (localRows[j]) >= MatrixInNumRows,
527 std::runtime_error,
"Ifpack2::SparseContainer::extract: localRows[j=" 528 << j <<
"] = " << localRows[j] <<
", which is out of the valid range. " 529 "This probably means that compute() has not yet been called.");
532 const size_t maxNumEntriesInRow = globalMatrix->getNodeMaxNumRowEntries();
533 Array<scalar_type> Values;
534 Array<local_ordinal_type> Indices;
535 Array<InverseScalar> Values_insert;
536 Array<InverseGlobalOrdinal> Indices_insert;
538 Values.resize (maxNumEntriesInRow);
539 Indices.resize (maxNumEntriesInRow);
540 Values_insert.resize (maxNumEntriesInRow);
541 Indices_insert.resize (maxNumEntriesInRow);
543 const InverseLocalOrdinal INVALID =
544 Teuchos::OrdinalTraits<InverseLocalOrdinal>::invalid ();
545 for (
size_t j = 0; j < numRows_; ++j) {
546 const local_ordinal_type localRow = localRows[j];
548 globalMatrix->getLocalRowCopy (localRow, Indices (), Values (), numEntries);
550 size_t num_entries_found = 0;
551 for (
size_t k = 0; k < numEntries; ++k) {
552 const local_ordinal_type localCol = Indices[k];
563 if (static_cast<size_t> (localCol) >= MatrixInNumRows) {
568 InverseLocalOrdinal jj = INVALID;
569 for (
size_t kk = 0; kk < numRows_; ++kk) {
570 if (localRows[kk] == localCol) {
576 Indices_insert[num_entries_found] = localMap_->getGlobalElement (jj);
577 Values_insert[num_entries_found] = Values[k];
581 diagBlock_->insertGlobalValues (j, Indices_insert (0, num_entries_found),
582 Values_insert (0, num_entries_found));
588 diagBlock_->fillComplete ();
594 #include "Ifpack2_ILUT.hpp" 600 #define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \ 601 template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \ 602 Ifpack2::ILUT<Tpetra::RowMatrix<S,LO,GO,N> > >; 604 #endif // IFPACK2_SPARSECONTAINER_HPP virtual 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
Compute Y := alpha * M^{-1} X + beta*Y.
Definition: Ifpack2_SparseContainer_def.hpp:194
virtual bool isComputed() const
Whether the container has been successfully computed.
Definition: Ifpack2_SparseContainer_def.hpp:97
Store and solve a local sparse linear problem.
Definition: Ifpack2_SparseContainer_decl.hpp:130
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to the given FancyOStream.
Definition: Ifpack2_SparseContainer_def.hpp:498
virtual void weightedApply(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, const Tpetra::Vector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &D, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Compute Y := alpha * diag(D) * M^{-1} (diag(D) * X) + beta*Y.
Definition: Ifpack2_SparseContainer_def.hpp:304
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input (global) matrix.
Definition: Ifpack2_SparseContainer_decl.hpp:154
virtual void initialize()
Do all set-up operations that only require matrix structure.
Definition: Ifpack2_SparseContainer_def.hpp:111
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input (global) matrix.
Definition: Ifpack2_SparseContainer_decl.hpp:152
Teuchos::ArrayView< const local_ordinal_type > getLocalRows() const
Local indices of the rows of the input matrix that belong to this block.
Definition: Ifpack2_Container.hpp:169
virtual std::ostream & print(std::ostream &os) const
Print information about this object to the given output stream.
Definition: Ifpack2_SparseContainer_def.hpp:466
virtual size_t getNumRows() const
The number of rows in the local matrix on the calling process.
Definition: Ifpack2_SparseContainer_def.hpp:79
virtual std::string description() const
A one-line description of this object.
Definition: Ifpack2_SparseContainer_def.hpp:476
virtual bool isInitialized() const
Whether the container has been successfully initialized.
Definition: Ifpack2_SparseContainer_def.hpp:90
Definition: Ifpack2_Details_Amesos2Wrapper_decl.hpp:63
virtual ~SparseContainer()
Destructor (declared virtual for memory safety of derived classes).
Definition: Ifpack2_SparseContainer_def.hpp:74
Interface for creating and solving a local linear problem.
Definition: Ifpack2_Container.hpp:103
SparseContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::ArrayView< const local_ordinal_type > &localRows)
Constructor.
Definition: Ifpack2_SparseContainer_def.hpp:59
virtual void setParameters(const Teuchos::ParameterList &List)
Set all necessary parameters.
Definition: Ifpack2_SparseContainer_def.hpp:104
MatrixType::scalar_type scalar_type
The type of entries in the input (global) matrix.
Definition: Ifpack2_SparseContainer_decl.hpp:150
virtual void compute()
Extract the local diagonal block and prepare the solver.
Definition: Ifpack2_SparseContainer_def.hpp:145
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
Teuchos::RCP< const row_matrix_type > getMatrix() const
The input matrix to the constructor.
Definition: Ifpack2_Container.hpp:147
Ifpack2::SparseContainer class declaration.
Implementation detail of Ifpack2::Container subclasses.
Definition: Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85