43#ifndef IFPACK2_SPARSECONTAINER_DEF_HPP
44#define IFPACK2_SPARSECONTAINER_DEF_HPP
49#include "Teuchos_DefaultMpiComm.hpp"
51#include "Teuchos_DefaultSerialComm.hpp"
53#include "Teuchos_TestForException.hpp"
58template<
class MatrixType,
class InverseType>
61 const Teuchos::Array<Teuchos::Array<LO> >& partitions,
62 const Teuchos::RCP<const import_type>&,
64 ContainerImpl<MatrixType, InverseScalar> (matrix, partitions, pointIndexed),
66 localComm_ (Teuchos::rcp (new Teuchos::MpiComm<int> (MPI_COMM_SELF)))
68 localComm_ (Teuchos::rcp (new Teuchos::SerialComm<int> ()))
73template<
class MatrixType,
class InverseType>
78template<
class MatrixType,
class InverseType>
85template<
class MatrixType,
class InverseType>
96 diagBlocks_.assign(this->numBlocks_, Teuchos::null);
97 Inverses_.assign(this->numBlocks_, Teuchos::null);
99 this->IsInitialized_ =
true;
100 this->IsComputed_ =
false;
104template<
class MatrixType,
class InverseType>
107 this->IsComputed_ =
false;
108 if (!this->isInitialized ()) {
118 for(
int i = 0; i < this->numBlocks_; i++)
120 Inverses_[i]->setParameters(List_);
121 Inverses_[i]->initialize ();
123 for(
int i = 0; i < this->numBlocks_; i++)
124 Inverses_[i]->compute ();
125 this->IsComputed_ =
true;
129template<
class MatrixType,
class InverseType>
132 for(
auto inv : Inverses_)
140template<
class MatrixType,
class InverseType>
145 Teuchos::ETransp mode,
147 InverseScalar beta)
const
149 TEUCHOS_TEST_FOR_EXCEPTION(
150 Inverses_[blockIndex]->getDomainMap()->getLocalNumElements() != X.getLocalLength(),
151 std::logic_error,
"Ifpack2::SparseContainer::apply: Inverse_ "
152 "operator and X have incompatible dimensions (" <<
153 Inverses_[blockIndex]->getDomainMap()->getLocalNumElements() <<
" resp. "
154 << X.getLocalLength() <<
"). Please report this bug to "
155 "the Ifpack2 developers.");
156 TEUCHOS_TEST_FOR_EXCEPTION(
157 Inverses_[blockIndex]->getRangeMap()->getLocalNumElements() != Y.getLocalLength(),
158 std::logic_error,
"Ifpack2::SparseContainer::apply: Inverse_ "
159 "operator and Y have incompatible dimensions (" <<
160 Inverses_[blockIndex]->getRangeMap()->getLocalNumElements() <<
" resp. "
161 << Y.getLocalLength() <<
"). Please report this bug to "
162 "the Ifpack2 developers.");
163 Inverses_[blockIndex]->apply(X, Y, mode, alpha, beta);
166template<
class MatrixType,
class InverseType>
168apply (ConstHostView X,
171 Teuchos::ETransp mode,
175 using Teuchos::ArrayView;
189 size_t numVecs = X.extent(1);
191 TEUCHOS_TEST_FOR_EXCEPTION(
192 ! this->IsComputed_, std::runtime_error,
"Ifpack2::SparseContainer::apply: "
193 "You must have called the compute() method before you may call apply(). "
194 "You may call the apply() method as many times as you want after calling "
195 "compute() once, but you must have called compute() at least once.");
196 TEUCHOS_TEST_FOR_EXCEPTION(
197 X.extent(1) != Y.extent(1), std::runtime_error,
198 "Ifpack2::SparseContainer::apply: X and Y have different numbers of "
199 "vectors. X has " << X.extent(1)
200 <<
", but Y has " << Y.extent(1) <<
".");
206 const LO numRows = this->blockSizes_[blockIndex];
235 for(LO i = 0; i < this->numBlocks_; i++)
236 invX.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
237 for(LO i = 0; i < this->numBlocks_; i++)
238 invY.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
240 inverse_mv_type& X_local = invX[blockIndex];
241 TEUCHOS_TEST_FOR_EXCEPTION(
242 X_local.getLocalLength() !=
size_t(numRows * this->scalarsPerRow_), std::logic_error,
243 "Ifpack2::SparseContainer::apply: "
244 "X_local has length " << X_local.getLocalLength() <<
", which does "
245 "not match numRows = " << numRows * this->scalarsPerRow_ <<
". Please report this bug to "
246 "the Ifpack2 developers.");
247 const ArrayView<const LO> blockRows = this->getBlockRows(blockIndex);
248 if(this->scalarsPerRow_ == 1)
249 mvgs.gatherMVtoView(X_local, X, blockRows);
251 mvgs.gatherMVtoViewBlock(X_local, X, blockRows, this->scalarsPerRow_);
258 inverse_mv_type& Y_local = invY[blockIndex];
259 TEUCHOS_TEST_FOR_EXCEPTION(
260 Y_local.getLocalLength () !=
size_t(numRows * this->scalarsPerRow_), std::logic_error,
261 "Ifpack2::SparseContainer::apply: "
262 "Y_local has length " << Y_local.getLocalLength () <<
", which does "
263 "not match numRows = " << numRows * this->scalarsPerRow_ <<
". Please report this bug to "
264 "the Ifpack2 developers.");
266 if(this->scalarsPerRow_ == 1)
267 mvgs.gatherMVtoView(Y_local, Y, blockRows);
269 mvgs.gatherMVtoViewBlock(Y_local, Y, blockRows, this->scalarsPerRow_);
273 this->solveBlockMV(X_local, Y_local, blockIndex, mode,
274 InverseScalar(alpha), InverseScalar(beta));
279 if(this->scalarsPerRow_ == 1)
280 mvgs.scatterMVtoView(Y, Y_local, blockRows);
282 mvgs.scatterMVtoViewBlock(Y, Y_local, blockRows, this->scalarsPerRow_);
286template<
class MatrixType,
class InverseType>
292 Teuchos::ETransp mode,
296 using Teuchos::ArrayView;
297 using Teuchos::Range1D;
300 typedef Teuchos::ScalarTraits<InverseScalar> STS;
312 typedef Tpetra::Vector<InverseScalar, InverseLocalOrdinal, InverseGlobalOrdinal, InverseNode> inverse_vector_type;
315 const size_t numVecs = X.extent(1);
317 TEUCHOS_TEST_FOR_EXCEPTION(
318 ! this->IsComputed_, std::runtime_error,
"Ifpack2::SparseContainer::"
319 "weightedApply: You must have called the compute() method before you may "
320 "call apply(). You may call the apply() method as many times as you want "
321 "after calling compute() once, but you must have called compute() at least "
323 TEUCHOS_TEST_FOR_EXCEPTION(
324 X.extent(1) != Y.extent(1), std::runtime_error,
325 "Ifpack2::SparseContainer::weightedApply: X and Y have different numbers "
326 "of vectors. X has " << X.extent(1) <<
", but Y has "
327 << Y.extent(1) <<
".");
330 TEUCHOS_TEST_FOR_EXCEPTION(
331 this->scalarsPerRow_ > 1, std::logic_error,
"Ifpack2::SparseContainer::weightedApply: "
332 "Use of block rows isn't allowed in overlapping Jacobi (the only method that uses weightedApply");
362 const LO numRows = this->blockSizes_[blockIndex];
366 for(LO i = 0; i < this->numBlocks_; i++)
367 invX.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
368 for(LO i = 0; i < this->numBlocks_; i++)
369 invY.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
371 inverse_mv_type& X_local = invX[blockIndex];
372 const ArrayView<const LO> blockRows = this->getBlockRows(blockIndex);
373 mvgs.gatherMVtoView(X_local, X, blockRows);
380 inverse_mv_type Y_local = invY[blockIndex];
381 TEUCHOS_TEST_FOR_EXCEPTION(
382 Y_local.getLocalLength() !=
size_t(numRows), std::logic_error,
383 "Ifpack2::SparseContainer::weightedApply: "
384 "Y_local has length " << X_local.getLocalLength() <<
", which does "
385 "not match numRows = " << numRows <<
". Please report this bug to "
386 "the Ifpack2 developers.");
387 mvgs.gatherMVtoView(Y_local, Y, blockRows);
399 inverse_vector_type D_local(Inverses_[blockIndex]->getDomainMap());
400 TEUCHOS_TEST_FOR_EXCEPTION(
401 D_local.getLocalLength() !=
size_t(this->blockSizes_[blockIndex]), std::logic_error,
402 "Ifpack2::SparseContainer::weightedApply: "
403 "D_local has length " << X_local.getLocalLength () <<
", which does "
404 "not match numRows = " << this->blockSizes_[blockIndex] <<
". Please report this bug to "
405 "the Ifpack2 developers.");
406 mvgs.gatherMVtoView(D_local, D, blockRows);
407 inverse_mv_type X_scaled(Inverses_[blockIndex]->getDomainMap(), numVecs);
408 X_scaled.elementWiseMultiply(STS::one(), D_local, X_local, STS::zero());
415 inverse_mv_type* Y_temp;
416 if (InverseScalar(beta) == STS::zero ()) {
419 Y_temp =
new inverse_mv_type(Inverses_[blockIndex]->getRangeMap(), numVecs);
422 Inverses_[blockIndex]->apply(X_scaled, *Y_temp, mode);
428 Y_local.elementWiseMultiply(alpha, D_local, *Y_temp, beta);
429 if(Y_temp != &Y_local)
433 mvgs.scatterMVtoView(Y, Y_local, blockRows);
437template<
class MatrixType,
class InverseType>
440 Teuchos::FancyOStream fos(Teuchos::rcp(&os,
false));
441 fos.setOutputToRootOnly(0);
447template<
class MatrixType,
class InverseType>
450 std::ostringstream oss;
451 oss <<
"\"Ifpack2::SparseContainer\": {";
452 if (this->isInitialized()) {
453 if (this->isComputed()) {
454 oss <<
"status = initialized, computed";
457 oss <<
"status = initialized, not computed";
461 oss <<
"status = not initialized, not computed";
463 for(
int i = 0; i < this->numBlocks_; i++)
465 oss <<
", Block Inverse " << i <<
": {";
466 oss << Inverses_[i]->description();
474template<
class MatrixType,
class InverseType>
478 if(verbLevel==Teuchos::VERB_NONE)
return;
479 os <<
"================================================================================" << endl;
480 os <<
"Ifpack2::SparseContainer" << endl;
481 for(
int i = 0; i < this->numBlocks_; i++)
483 os <<
"Block " << i <<
" rows: = " << this->blockSizes_[i] << endl;
485 os <<
"isInitialized() = " << this->IsInitialized_ << endl;
486 os <<
"isComputed() = " << this->IsComputed_ << endl;
487 os <<
"================================================================================" << endl;
492template<
class MatrixType,
class InverseType>
496 using Teuchos::Array;
497 using Teuchos::ArrayView;
499 const LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
506 Teuchos::Array<InverseGlobalOrdinal> indicesToInsert;
507 Teuchos::Array<InverseScalar> valuesToInsert;
508 if(this->scalarsPerRow_ > 1)
510 Array<LO> colToBlockOffset(this->inputBlockMatrix_->getLocalNumCols(), INVALID);
511 for(
int i = 0; i < this->numBlocks_; i++)
514 LO blockStart = this->blockOffsets_[i];
515 LO blockSize = this->blockSizes_[i];
516 LO blockPointSize = this->blockSizes_[i] * this->scalarsPerRow_;
517 LO blockEnd = blockStart + blockSize;
518 ArrayView<const LO> blockRows = this->getBlockRows(i);
522 for(LO j = 0; j < blockSize; j++)
524 LO localCol = this->translateRowToCol(blockRows[j]);
525 colToBlockOffset[localCol] = blockStart + j;
529 Array<size_t> rowEntryCounts(blockPointSize, 0);
532 using inds_type =
typename block_crs_matrix_type::local_inds_host_view_type;
533 using vals_type =
typename block_crs_matrix_type::values_host_view_type;
534 for(LO blockRow = 0; blockRow < blockSize; blockRow++)
539 LO inputRow = this->blockRows_[blockStart + blockRow];
540 this->inputBlockMatrix_->getLocalRowView(inputRow, indices, values);
541 LO numEntries = (LO) indices.size();
542 for(LO br = 0; br < this->bcrsBlockSize_; br++)
544 for(LO k = 0; k < numEntries; k++)
546 LO colOffset = colToBlockOffset[indices[k]];
547 if(blockStart <= colOffset && colOffset < blockEnd)
549 rowEntryCounts[blockRow * this->bcrsBlockSize_ + br] += this->bcrsBlockSize_;
555 RCP<InverseMap> tempMap(
new InverseMap(blockPointSize, 0, this->localComm_));
556 diagBlocks_[i] = rcp(
new InverseCrs(tempMap, rowEntryCounts));
557 Inverses_[i] = rcp(
new InverseType(diagBlocks_[i]));
559 for(LO blockRow = 0; blockRow < blockSize; blockRow++)
564 LO inputRow = this->blockRows_[blockStart + blockRow];
565 this->inputBlockMatrix_->getLocalRowView(inputRow, indices, values);
566 LO numEntries = (LO) indices.size();
567 for(LO br = 0; br < this->bcrsBlockSize_; br++)
569 indicesToInsert.clear();
570 valuesToInsert.clear();
571 for(LO k = 0; k < numEntries; k++)
573 LO colOffset = colToBlockOffset[indices[k]];
574 if(blockStart <= colOffset && colOffset < blockEnd)
576 LO blockCol = colOffset - blockStart;
578 for(LO bc = 0; bc < this->bcrsBlockSize_; bc++)
580 indicesToInsert.push_back(blockCol * this->bcrsBlockSize_ + bc);
581 valuesToInsert.push_back(values[k * this->bcrsBlockSize_ * this->bcrsBlockSize_ + bc * this->bcrsBlockSize_ + br]);
585 InverseGlobalOrdinal rowToInsert = blockRow * this->bcrsBlockSize_ + br;
586 if(indicesToInsert.size())
587 diagBlocks_[i]->insertGlobalValues(rowToInsert, indicesToInsert(), valuesToInsert());
590 diagBlocks_[i]->fillComplete();
597 Array<LO> colToBlockOffset(this->inputMatrix_->getLocalNumCols() * this->bcrsBlockSize_, INVALID);
598 for(
int i = 0; i < this->numBlocks_; i++)
601 LO blockStart = this->blockOffsets_[i];
602 LO blockSize = this->blockSizes_[i];
603 LO blockEnd = blockStart + blockSize;
604 ArrayView<const LO> blockRows = this->getBlockRows(i);
608 for(LO j = 0; j < blockSize; j++)
611 LO localCol = this->translateRowToCol(blockRows[j]);
612 colToBlockOffset[localCol] = blockStart + j;
614 Teuchos::Array<size_t> rowEntryCounts(blockSize, 0);
615 for(LO j = 0; j < blockSize; j++)
617 rowEntryCounts[j] = this->getInputRowView(this->blockRows_[blockStart + j]).size();
619 RCP<InverseMap> tempMap(
new InverseMap(blockSize, 0, this->localComm_));
620 diagBlocks_[i] = rcp(
new InverseCrs(tempMap, rowEntryCounts));
621 Inverses_[i] = rcp(
new InverseType(diagBlocks_[i]));
622 for(LO blockRow = 0; blockRow < blockSize; blockRow++)
624 valuesToInsert.clear();
625 indicesToInsert.clear();
627 LO inputSplitRow = this->blockRows_[blockStart + blockRow];
628 auto rowView = this->getInputRowView(inputSplitRow);
629 for(
size_t k = 0; k < rowView.size(); k++)
631 LO colOffset = colToBlockOffset[rowView.ind(k)];
632 if(blockStart <= colOffset && colOffset < blockEnd)
634 LO blockCol = colOffset - blockStart;
635 indicesToInsert.push_back(blockCol);
636 valuesToInsert.push_back(rowView.val(k));
639 if(indicesToInsert.size())
640 diagBlocks_[i]->insertGlobalValues(blockRow, indicesToInsert(), valuesToInsert());
642 diagBlocks_[i]->fillComplete ();
647template<
typename MatrixType,
typename InverseType>
651#ifdef HAVE_IFPACK2_AMESOS2
653 if(std::is_same<InverseType, ILUTInverse>::value)
657 else if(std::is_same<InverseType, AmesosInverse>::value)
659 return "SparseAmesos";
663 throw std::logic_error(
"InverseType for SparseContainer must be Ifpack2::ILUT or Details::Amesos2Wrapper");
668 constexpr bool inverseTypeIsILUT = std::is_same<InverseType, ILUTInverse>::value;
669 TEUCHOS_TEST_FOR_EXCEPTION(! inverseTypeIsILUT, std::logic_error,
670 "InverseType for SparseContainer must be Ifpack2::ILUT<ROW>");
678#include "Ifpack2_ILUT.hpp"
679#ifdef HAVE_IFPACK2_AMESOS2
680#include "Ifpack2_Details_Amesos2Wrapper.hpp"
687#ifdef HAVE_IFPACK2_AMESOS2
688# define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
689 template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \
690 Ifpack2::ILUT<Tpetra::RowMatrix<S,LO,GO,N> > >; \
691 template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \
692 Ifpack2::Details::Amesos2Wrapper<Tpetra::RowMatrix<S,LO,GO,N> > >;
694# define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
695 template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S,LO,GO,N>, \
696 Ifpack2::ILUT<Tpetra::RowMatrix<S, LO, GO, N> > >;
Ifpack2::SparseContainer class declaration.
The implementation of the numerical features of Container (Jacobi, Gauss-Seidel, SGS)....
Definition: Ifpack2_Container_decl.hpp:344
typename mv_type::dual_view_type::t_host HostView
Definition: Ifpack2_Container_decl.hpp:139
Wrapper class for direct solvers in Amesos2.
Definition: Ifpack2_Details_Amesos2Wrapper_decl.hpp:111
Implementation detail of Ifpack2::Container subclasses.
Definition: Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85
ILUT (incomplete LU factorization with threshold) of a Tpetra sparse matrix.
Definition: Ifpack2_ILUT_decl.hpp:100
Store and solve a local sparse linear problem.
Definition: Ifpack2_SparseContainer_decl.hpp:136
virtual void apply(ConstHostView X, HostView Y, int blockIndex, Teuchos::ETransp mode=Teuchos::NO_TRANS, SC alpha=Teuchos::ScalarTraits< SC >::one(), SC beta=Teuchos::ScalarTraits< SC >::zero()) const
Compute Y := alpha * M^{-1} X + beta*Y.
Definition: Ifpack2_SparseContainer_def.hpp:168
SparseContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< Teuchos::Array< LO > > &partitions, const Teuchos::RCP< const import_type > &importer, bool pointIndexed)
Constructor.
Definition: Ifpack2_SparseContainer_def.hpp:60
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:475
virtual void weightedApply(ConstHostView X, HostView Y, ConstHostView W, int blockIndex, Teuchos::ETransp mode=Teuchos::NO_TRANS, SC alpha=Teuchos::ScalarTraits< SC >::one(), SC beta=Teuchos::ScalarTraits< SC >::zero()) const
Compute Y := alpha * diag(D) * M^{-1} (diag(D) * X) + beta*Y.
Definition: Ifpack2_SparseContainer_def.hpp:288
static std::string getName()
Get the name of this container type for Details::constructContainer()
Definition: Ifpack2_SparseContainer_def.hpp:648
virtual void initialize()
Do all set-up operations that only require matrix structure.
Definition: Ifpack2_SparseContainer_def.hpp:86
void clearBlocks()
Definition: Ifpack2_SparseContainer_def.hpp:130
virtual void compute()
Initialize and compute all blocks.
Definition: Ifpack2_SparseContainer_def.hpp:105
virtual ~SparseContainer()
Destructor (declared virtual for memory safety of derived classes).
Definition: Ifpack2_SparseContainer_def.hpp:75
virtual std::ostream & print(std::ostream &os) const
Print information about this object to the given output stream.
Definition: Ifpack2_SparseContainer_def.hpp:438
virtual std::string description() const
A one-line description of this object.
Definition: Ifpack2_SparseContainer_def.hpp:448
virtual void setParameters(const Teuchos::ParameterList &List)
Set all necessary parameters.
Definition: Ifpack2_SparseContainer_def.hpp:79
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74