43 #ifndef IFPACK2_SPARSECONTAINER_DEF_HPP
44 #define IFPACK2_SPARSECONTAINER_DEF_HPP
51 #include "Teuchos_DefaultSerialComm.hpp"
58 template<
class MatrixType,
class InverseType>
64 scalar_type DampingFactor) :
65 Container<MatrixType> (matrix, partitions, importer, OverlapLevel,
67 IsInitialized_ (false),
70 localComm_ (Teuchos::
rcp (new Teuchos::MpiComm<int> (MPI_COMM_SELF)))
72 localComm_ (Teuchos::
rcp (new Teuchos::SerialComm<int> ()))
75 global_ordinal_type indexBase = 0;
77 localMaps_.emplace_back(this->blockRows_[i], indexBase, localComm_);
81 template<
class MatrixType,
class InverseType>
85 Container<MatrixType> (matrix, localRows),
86 IsInitialized_(false),
89 localComm_ (Teuchos::
rcp(new Teuchos::MpiComm<int>(MPI_COMM_SELF)))
91 localComm_ (Teuchos::
rcp(new Teuchos::SerialComm<int>()))
94 global_ordinal_type indexBase = 0;
95 localMaps_.emplace_back(this->
blockRows_[0], indexBase, localComm_);
99 template<
class MatrixType,
class InverseType>
102 for(
auto inv : Inverses_)
107 template<
class MatrixType,
class InverseType>
110 return IsInitialized_;
114 template<
class MatrixType,
class InverseType>
121 template<
class MatrixType,
class InverseType>
128 template<
class MatrixType,
class InverseType>
137 IsInitialized_ =
false;
142 diagBlocks_.reserve(this->numBlocks_);
143 Inverses_.reserve(this->numBlocks_);
144 for(
int i = 0; i < this->numBlocks_; i++)
149 RCP<InverseMap> tempMap(
new InverseMap(this->blockRows_[i], 0, localComm_));
150 diagBlocks_.emplace_back(
new InverseCrs(tempMap, 0));
154 Inverses_.push_back(ptr(
new InverseType(diagBlocks_[i])));
156 IsInitialized_ =
true;
160 template<
class MatrixType,
class InverseType>
164 if (! this->isInitialized ()) {
174 for(
int i = 0; i < this->numBlocks_; i++)
175 Inverses_[i]->initialize ();
176 for(
int i = 0; i < this->numBlocks_; i++)
177 Inverses_[i]->compute ();
182 template<
class MatrixType,
class InverseType>
185 for(
auto inv : Inverses_)
193 template<
class MatrixType,
class InverseType>
201 InverseScalar beta)
const
204 Inverses_[blockIndex]->getDomainMap()->getNodeNumElements() != X.getLocalLength(),
205 std::logic_error,
"Ifpack2::SparseContainer::apply: Inverse_ "
206 "operator and X have incompatible dimensions (" <<
207 Inverses_[blockIndex]->getDomainMap()->getNodeNumElements() <<
" resp. "
208 << X.getLocalLength() <<
"). Please report this bug to "
209 "the Ifpack2 developers.");
211 Inverses_[blockIndex]->getRangeMap()->getNodeNumElements() != Y.getLocalLength(),
212 std::logic_error,
"Ifpack2::SparseContainer::apply: Inverse_ "
213 "operator and Y have incompatible dimensions (" <<
214 Inverses_[blockIndex]->getRangeMap()->getNodeNumElements() <<
" resp. "
215 << Y.getLocalLength() <<
"). Please report this bug to "
216 "the Ifpack2 developers.");
217 Inverses_[blockIndex]->apply(X, Y, mode, alpha, beta);
220 template<
class MatrixType,
class InverseType>
228 scalar_type beta)
const
246 size_t numVecs = X.extent(1);
249 ! IsComputed_, std::runtime_error,
"Ifpack2::SparseContainer::apply: "
250 "You must have called the compute() method before you may call apply(). "
251 "You may call the apply() method as many times as you want after calling "
252 "compute() once, but you must have called compute() at least once.");
254 X.extent(1) != Y.extent(1), std::runtime_error,
255 "Ifpack2::SparseContainer::apply: X and Y have different numbers of "
256 "vectors. X has " << X.extent(1)
257 <<
", but Y has " << Y.extent(1) <<
".");
263 const local_ordinal_type numRows_ = this->blockRows_[blockIndex];
292 for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
293 invX.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
294 for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
295 invY.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
297 inverse_mv_type& X_local = invX[blockIndex];
299 X_local.getLocalLength() != (size_t) numRows_, std::logic_error,
300 "Ifpack2::SparseContainer::apply: "
301 "X_local has length " << X_local.getLocalLength() <<
", which does "
302 "not match numRows_ = " << numRows_ <<
". Please report this bug to "
303 "the Ifpack2 developers.");
304 const ArrayView<const local_ordinal_type> localRows = this->getLocalRows(blockIndex);
305 mvgs.gatherMVtoView(X_local, X, localRows);
312 inverse_mv_type& Y_local = invY[blockIndex];
314 Y_local.getLocalLength () != (size_t) numRows_, std::logic_error,
315 "Ifpack2::SparseContainer::apply: "
316 "Y_local has length " << Y_local.getLocalLength () <<
", which does "
317 "not match numRows_ = " << numRows_ <<
". Please report this bug to "
318 "the Ifpack2 developers.");
319 mvgs.gatherMVtoView(Y_local, Y, localRows);
323 this->applyImpl(X_local, Y_local, blockIndex, stride, mode, as<InverseScalar>(alpha),
324 as<InverseScalar>(beta));
328 mvgs.scatterMVtoView(Y, Y_local, localRows);
332 template<
class MatrixType,
class InverseType>
341 scalar_type beta)
const
350 using Teuchos::rcp_const_cast;
365 typedef Tpetra::Vector<InverseScalar, InverseLocalOrdinal, InverseGlobalOrdinal, InverseNode> inverse_vector_type;
368 const size_t numVecs = X.extent(1);
371 ! IsComputed_, std::runtime_error,
"Ifpack2::SparseContainer::"
372 "weightedApply: You must have called the compute() method before you may "
373 "call apply(). You may call the apply() method as many times as you want "
374 "after calling compute() once, but you must have called compute() at least "
377 X.extent(1) != Y.extent(1), std::runtime_error,
378 "Ifpack2::SparseContainer::weightedApply: X and Y have different numbers "
379 "of vectors. X has " << X.extent(1) <<
", but Y has "
380 << Y.extent(1) <<
".");
410 const local_ordinal_type numRows = this->blockRows_[blockIndex];
414 for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
415 invX.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
416 for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
417 invY.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
419 inverse_mv_type& X_local = invX[blockIndex];
420 const ArrayView<const local_ordinal_type> localRows = this->getLocalRows(blockIndex);
421 mvgs.gatherMVtoView(X_local, X, localRows);
428 inverse_mv_type Y_local = invY[blockIndex];
430 Y_local.getLocalLength() != (size_t) numRows, std::logic_error,
431 "Ifpack2::SparseContainer::weightedApply: "
432 "Y_local has length " << X_local.getLocalLength() <<
", which does "
433 "not match numRows_ = " << numRows <<
". Please report this bug to "
434 "the Ifpack2 developers.");
435 mvgs.gatherMVtoView(Y_local, Y, localRows);
447 inverse_vector_type D_local(Inverses_[blockIndex]->getDomainMap());
449 D_local.getLocalLength() != (size_t) this->blockRows_[blockIndex], std::logic_error,
450 "Ifpack2::SparseContainer::weightedApply: "
451 "D_local has length " << X_local.getLocalLength () <<
", which does "
452 "not match numRows_ = " << this->blockRows_[blockIndex] <<
". Please report this bug to "
453 "the Ifpack2 developers.");
454 mvgs.gatherMVtoView(D_local, D, localRows);
455 inverse_mv_type X_scaled(Inverses_[blockIndex]->getDomainMap(), numVecs);
456 X_scaled.elementWiseMultiply(STS::one(), D_local, X_local, STS::zero());
463 Ptr<inverse_mv_type> Y_temp;
464 bool deleteYT =
false;
465 if (beta == STS::zero ()) {
466 Y_temp = ptr(&Y_local);
468 Y_temp = ptr(
new inverse_mv_type(Inverses_[blockIndex]->getRangeMap(), numVecs));
472 Inverses_[blockIndex]->apply(X_scaled, *Y_temp, mode);
478 Y_local.elementWiseMultiply(alpha, D_local, *Y_temp, beta);
483 mvgs.scatterMVtoView(Y, Y_local, localRows);
487 template<
class MatrixType,
class InverseType>
497 template<
class MatrixType,
class InverseType>
500 std::ostringstream oss;
501 oss <<
"\"Ifpack2::SparseContainer\": {";
502 if (isInitialized()) {
504 oss <<
"status = initialized, computed";
507 oss <<
"status = initialized, not computed";
511 oss <<
"status = not initialized, not computed";
513 for(
int i = 0; i < this->numBlocks_; i++)
515 oss <<
", Block Inverse " << i <<
": {";
516 oss << Inverses_[i]->description();
524 template<
class MatrixType,
class InverseType>
529 os <<
"================================================================================" << endl;
530 os <<
"Ifpack2::SparseContainer" << endl;
531 for(
int i = 0; i < this->numBlocks_; i++)
533 os <<
"Block " << i <<
" rows: = " << this->blockRows_[i] << endl;
535 os <<
"isInitialized() = " << IsInitialized_ << endl;
536 os <<
"isComputed() = " << IsComputed_ << endl;
537 os <<
"================================================================================" << endl;
542 template<
class MatrixType,
class InverseType>
549 auto& A = *this->inputMatrix_;
550 const size_t MatrixInNumRows = A.getNodeNumRows ();
553 for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
555 const ArrayView<const local_ordinal_type> localRows = this->getLocalRows(i);
556 for (local_ordinal_type j = 0; j < localRows.size(); j++)
560 static_cast<size_t> (localRows[j]) >= MatrixInNumRows,
561 std::runtime_error,
"Ifpack2::SparseContainer::extract: localRows[j="
562 << j <<
"] = " << localRows[j] <<
", which is out of the valid range. "
563 "This probably means that compute() has not yet been called.");
567 const size_t maxNumEntriesInRow = A.getNodeMaxNumRowEntries();
568 Array<scalar_type> Values;
569 Array<local_ordinal_type> Indices;
570 Array<InverseScalar> Values_insert;
571 Array<InverseGlobalOrdinal> Indices_insert;
573 Values.resize (maxNumEntriesInRow);
574 Indices.resize (maxNumEntriesInRow);
575 Values_insert.resize (maxNumEntriesInRow);
576 Indices_insert.resize (maxNumEntriesInRow);
579 for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
581 const local_ordinal_type numRows_ = this->blockRows_[i];
582 const ArrayView<const local_ordinal_type> localRows = this->getLocalRows(i);
583 for (local_ordinal_type j = 0; j < numRows_; j++)
585 const local_ordinal_type localRow = this->partitions_[this->partitionIndices_[i] + j];
587 A.getLocalRowCopy(localRow, Indices(), Values(), numEntries);
588 size_t num_entries_found = 0;
589 for(
size_t k = 0; k < numEntries; ++k)
591 const local_ordinal_type localCol = Indices[k];
601 if(static_cast<size_t> (localCol) >= MatrixInNumRows)
605 InverseLocalOrdinal jj = INVALID;
606 for(local_ordinal_type kk = 0; kk < numRows_; kk++)
608 if(localRows[kk] == localCol)
613 Indices_insert[num_entries_found] = localMaps_[i].getGlobalElement(jj);
614 Values_insert[num_entries_found] = Values[k];
618 diagBlocks_[i]->insertGlobalValues(j, Indices_insert (0, num_entries_found),
619 Values_insert (0, num_entries_found));
624 diagBlocks_[i]->fillComplete ();
628 template<
typename MatrixType,
typename InverseType>
632 #ifdef HAVE_IFPACK2_AMESOS2
634 if(std::is_same<InverseType, ILUTInverse>::value)
638 else if(std::is_same<InverseType, AmesosInverse>::value)
640 return "SparseAmesos";
644 throw std::logic_error(
"InverseType for SparseContainer must be Ifpack2::ILUT or Details::Amesos2Wrapper");
649 constexpr
bool inverseTypeIsILUT = std::is_same<InverseType, ILUTInverse>::value;
651 "InverseType for SparseContainer must be Ifpack2::ILUT<ROW>");
659 #include "Ifpack2_ILUT.hpp"
660 #ifdef HAVE_IFPACK2_AMESOS2
661 #include "Ifpack2_Details_Amesos2Wrapper.hpp"
668 #ifdef HAVE_IFPACK2_AMESOS2
669 # define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
670 template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \
671 Ifpack2::ILUT<Tpetra::RowMatrix<S,LO,GO,N> > >; \
672 template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \
673 Ifpack2::Details::Amesos2Wrapper<Tpetra::RowMatrix<S,LO,GO,N> > >;
675 # define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
676 template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S,LO,GO,N>, \
677 Ifpack2::ILUT<Tpetra::RowMatrix<S, LO, GO, N> > >;
679 #endif // IFPACK2_SPARSECONTAINER_DEF_HPP
static std::string getName()
Get the name of this container type for Details::constructContainer()
Definition: Ifpack2_SparseContainer_def.hpp:629
virtual bool isComputed() const
Whether the container has been successfully computed.
Definition: Ifpack2_SparseContainer_def.hpp:115
Store and solve a local sparse linear problem.
Definition: Ifpack2_SparseContainer_decl.hpp:134
int numBlocks_
The number of blocks (partitions) in the container.
Definition: Ifpack2_Container.hpp:421
Teuchos::Array< local_ordinal_type > blockRows_
Number of rows in each block.
Definition: Ifpack2_Container.hpp:425
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:525
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void initialize()
Do all set-up operations that only require matrix structure.
Definition: Ifpack2_SparseContainer_def.hpp:129
ILUT (incomplete LU factorization with threshold) of a Tpetra sparse matrix.
Definition: Ifpack2_ILUT_decl.hpp:91
virtual std::ostream & print(std::ostream &os) const
Print information about this object to the given output stream.
Definition: Ifpack2_SparseContainer_def.hpp:488
virtual void weightedApply(HostView &X, HostView &Y, HostView &W, int blockIndex, int stride, 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:334
virtual std::string description() const
A one-line description of this object.
Definition: Ifpack2_SparseContainer_def.hpp:498
virtual bool isInitialized() const
Whether the container has been successfully initialized.
Definition: Ifpack2_SparseContainer_def.hpp:108
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Wrapper class for direct solvers in Amesos2.
Definition: Ifpack2_Details_Amesos2Wrapper_decl.hpp:102
SparseContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< Teuchos::Array< local_ordinal_type > > &partitions, const Teuchos::RCP< const import_type > &importer, int OverlapLevel, scalar_type DampingFactor)
Constructor.
Definition: Ifpack2_SparseContainer_def.hpp:60
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
virtual void apply(HostView &X, HostView &Y, int blockIndex, int stride, 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:222
virtual ~SparseContainer()
Destructor (declared virtual for memory safety of derived classes).
Definition: Ifpack2_SparseContainer_def.hpp:100
TypeTo as(const TypeFrom &t)
Interface for creating and solving a local linear problem.
Definition: Ifpack2_Container.hpp:114
virtual void setParameters(const Teuchos::ParameterList &List)
Set all necessary parameters.
Definition: Ifpack2_SparseContainer_def.hpp:122
virtual void compute()
Initialize and compute all blocks.
Definition: Ifpack2_SparseContainer_def.hpp:161
Ifpack2::SparseContainer class declaration.
void clearBlocks()
Definition: Ifpack2_SparseContainer_def.hpp:183
Implementation detail of Ifpack2::Container subclasses.
Definition: Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85