43 #ifndef IFPACK2_CONTAINER_HPP 44 #define IFPACK2_CONTAINER_HPP 49 #include "Ifpack2_ConfigDefs.hpp" 50 #include "Tpetra_RowMatrix.hpp" 51 #include "Teuchos_Describable.hpp" 52 #include <Ifpack2_Partitioner.hpp> 53 #include <Tpetra_Map.hpp> 54 #include <Tpetra_Experimental_BlockCrsMatrix_decl.hpp> 55 #include <Teuchos_ParameterList.hpp> 56 #include <Teuchos_Time.hpp> 58 #include <type_traits> 60 #ifndef DOXYGEN_SHOULD_SKIP_THIS 65 #endif // DOXYGEN_SHOULD_SKIP_THIS 113 template<
class MatrixType>
118 typedef typename MatrixType::scalar_type scalar_type;
119 typedef typename MatrixType::local_ordinal_type local_ordinal_type;
120 typedef typename MatrixType::global_ordinal_type global_ordinal_type;
121 typedef typename MatrixType::node_type node_type;
122 typedef Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> mv_type;
123 typedef Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> vector_type;
124 typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
125 typedef Teuchos::ScalarTraits<scalar_type> STS;
126 typedef Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
128 typedef Tpetra::Experimental::BlockCrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> block_crs_matrix_type;
129 typedef Tpetra::RowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> row_matrix_type;
131 static_assert(std::is_same<MatrixType, row_matrix_type>::value,
132 "Ifpack2::Container: Please use MatrixType = Tpetra::RowMatrix.");
135 typedef typename Kokkos::Details::ArithTraits<scalar_type>::val_type
impl_scalar_type;
139 typedef typename mv_type::dual_view_type::t_host HostView;
149 Container (
const Teuchos::RCP<const row_matrix_type>& matrix,
150 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
151 const Teuchos::RCP<const import_type>& importer,
153 scalar_type DampingFactor) :
161 using Teuchos::Array;
162 using Teuchos::ArrayView;
170 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(
inputMatrix_);
190 Container (
const Teuchos::RCP<const row_matrix_type>& matrix,
191 const Teuchos::Array<local_ordinal_type>& localRows) :
207 const block_crs_matrix_type* global_bcrs =
208 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(
inputMatrix_).
get();
214 for(local_ordinal_type i = 0; i < localRows.size(); i++)
236 Teuchos::ArrayView<const local_ordinal_type>
getLocalRows(
int blockIndex)
const 238 return Teuchos::ArrayView<const local_ordinal_type>
253 void setBlockSizes(
const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions)
257 local_ordinal_type totalBlockRows = 0;
263 local_ordinal_type rowsInBlock = partitions[i].size();
266 totalBlockRows += rowsInBlock;
270 local_ordinal_type iter = 0;
280 void getMatDiag()
const 296 void DoJacobi(HostView& X, HostView& Y,
int stride)
const;
297 void DoOverlappingJacobi(HostView& X, HostView& Y, HostView& W,
int stride)
const;
298 void DoGaussSeidel(HostView& X, HostView& Y, HostView& Y2,
int stride)
const;
299 void DoSGS(HostView& X, HostView& Y, HostView& Y2,
int stride)
const;
302 virtual void setParameters (
const Teuchos::ParameterList& List) = 0;
327 Teuchos::ETransp mode = Teuchos::NO_TRANS,
328 scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one(),
329 scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero())
const = 0;
334 HostView XView = X.template getLocalView<Kokkos::HostSpace>();
335 HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
336 this->
apply (XView, YView, 0, X.getStride());
365 Teuchos::ETransp mode = Teuchos::NO_TRANS,
366 scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one(),
367 scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero())
const = 0;
374 HostView XView = X.template getLocalView<Kokkos::HostSpace>();
375 HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
376 HostView WView = W.template getLocalView<Kokkos::HostSpace>();
380 virtual void clearBlocks();
383 virtual std::ostream&
print (std::ostream& os)
const = 0;
405 mutable Teuchos::RCP<vector_type>
Diag_;
413 Teuchos::RCP<const Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type>>
Importer_;
427 template <
class MatrixType>
429 operator<< (std::ostream& os, const Ifpack2::Container<MatrixType>& obj)
431 return obj.print (os);
434 template <
class MatrixType>
435 void Container<MatrixType>::DoJacobi(HostView& X, HostView& Y,
int stride)
const 437 const scalar_type one = STS::one();
440 size_t numVecs = X.dimension_1();
442 for (local_ordinal_type i = 0; i < numBlocks_; i++)
445 if(blockRows_[i] != 1 || hasBlockCrs_)
447 if(blockRows_[i] == 0 )
449 apply(X, Y, i, stride, Teuchos::NO_TRANS, DampingFactor_, one);
453 local_ordinal_type LRID = partitions_[partitionIndices_[i]];
455 HostView diagView = Diag_->template getLocalView<Kokkos::HostSpace>();
456 impl_scalar_type d = (impl_scalar_type) one / diagView(LRID, 0);
457 for(
size_t nv = 0; nv < numVecs; nv++)
459 impl_scalar_type x = X(LRID, nv);
466 template <
class MatrixType>
467 void Container<MatrixType>::DoOverlappingJacobi(HostView& X, HostView& Y, HostView& W,
int stride)
const 470 for(local_ordinal_type i = 0; i < numBlocks_; i++)
473 if(blockRows_[i] == 0)
475 if(blockRows_[i] != 1)
476 weightedApply(X, Y, W, i, stride, Teuchos::NO_TRANS, DampingFactor_, STS::one());
480 template<
class MatrixType>
481 void Container<MatrixType>::DoGaussSeidel(HostView& X, HostView& Y, HostView& Y2,
int stride)
const 483 using Teuchos::Array;
484 using Teuchos::ArrayRCP;
485 using Teuchos::ArrayView;
489 using Teuchos::rcpFromRef;
491 const scalar_type one = STS::one();
492 const size_t Length = inputMatrix_->getNodeMaxNumRowEntries();
493 auto numVecs = X.dimension_1();
494 Array<scalar_type> Values;
495 Array<local_ordinal_type> Indices;
496 Indices.resize(Length);
498 Values.resize(bcrsBlockSize_ * bcrsBlockSize_ * Length);
503 HostView Resid(
"", X.dimension_0(), X.dimension_1());
504 for(local_ordinal_type i = 0; i < numBlocks_; i++)
506 if(blockRows_[i] > 1 || hasBlockCrs_)
508 if (blockRows_[i] == 0)
511 ArrayView<const local_ordinal_type> localRows = getLocalRows (i);
512 const size_t localNumRows = blockRows_[i];
513 for(
size_t j = 0; j < localNumRows; j++)
515 const local_ordinal_type LID = localRows[j];
517 inputMatrix_->getLocalRowCopy(LID, Indices(), Values(), NumEntries);
518 for(
size_t m = 0; m < numVecs; m++)
522 for (
int localR = 0; localR < bcrsBlockSize_; localR++)
523 Resid(LID * bcrsBlockSize_ + localR, m) = X(LID * bcrsBlockSize_ + localR, m);
524 for (
size_t k = 0; k < NumEntries; ++k)
526 const local_ordinal_type col = Indices[k];
527 for (
int localR = 0; localR < bcrsBlockSize_; localR++)
529 for(
int localC = 0; localC < bcrsBlockSize_; localC++)
531 Resid(LID * bcrsBlockSize_ + localR, m) -=
532 Values[k * bcrsBlockSize_ * bcrsBlockSize_ + localR + localC * bcrsBlockSize_]
533 * Y2(col * bcrsBlockSize_ + localC, m);
540 Resid(LID, m) = X(LID, m);
541 for (
size_t k = 0; k < NumEntries; ++k)
543 const local_ordinal_type col = Indices[k];
544 Resid(LID, m) -= Values[k] * Y2(col, m);
555 apply(Resid, Y2, i, stride, Teuchos::NO_TRANS, DampingFactor_, one);
557 else if(blockRows_[i] == 1)
561 local_ordinal_type LRID = partitionIndices_[i];
563 HostView diagView = Diag_->template getLocalView<Kokkos::HostSpace>();
564 impl_scalar_type d = (impl_scalar_type) one / diagView(LRID, 0);
565 for(
size_t m = 0; m < numVecs; m++)
567 impl_scalar_type x = X(LRID, m);
568 impl_scalar_type newy = x * d;
575 auto numMyRows = inputMatrix_->getNodeNumRows();
576 for (
size_t m = 0; m < numVecs; ++m)
578 for (
size_t i = 0; i < numMyRows * bcrsBlockSize_; ++i)
586 template<
class MatrixType>
587 void Container<MatrixType>::DoSGS(HostView& X, HostView& Y, HostView& Y2,
int stride)
const 589 using Teuchos::Array;
590 using Teuchos::ArrayRCP;
591 using Teuchos::ArrayView;
596 using Teuchos::rcpFromRef;
597 const scalar_type one = STS::one();
598 auto numVecs = X.dimension_1();
599 const size_t Length = inputMatrix_->getNodeMaxNumRowEntries();
600 Array<scalar_type> Values;
601 Array<local_ordinal_type> Indices(Length);
602 Values.resize(bcrsBlockSize_ * bcrsBlockSize_ * Length);
607 HostView Resid(
"", X.dimension_0(), X.dimension_1());
609 for(local_ordinal_type i = 0; i < numBlocks_; i++)
611 if(blockRows_[i] != 1 || hasBlockCrs_)
613 if(blockRows_[i] == 0)
616 ArrayView<const local_ordinal_type> localRows = getLocalRows(i);
617 for(local_ordinal_type j = 0; j < blockRows_[i]; j++)
619 const local_ordinal_type LID = localRows[j];
621 inputMatrix_->getLocalRowCopy(LID, Indices(), Values(), NumEntries);
623 for(
size_t m = 0; m < numVecs; m++)
627 for(
int localR = 0; localR < bcrsBlockSize_; localR++)
628 Resid(LID * bcrsBlockSize_ + localR, m) = X(LID * bcrsBlockSize_ + localR, m);
629 for(
size_t k = 0; k < NumEntries; ++k)
631 const local_ordinal_type col = Indices[k];
632 for (
int localR = 0; localR < bcrsBlockSize_; localR++)
634 for(
int localC = 0; localC < bcrsBlockSize_; localC++)
635 Resid(LID * bcrsBlockSize_ + localR, m) -=
636 Values[k * (bcrsBlockSize_ * bcrsBlockSize_) + (localR + localC * bcrsBlockSize_)]
637 * Y2(col * bcrsBlockSize_ + localC, m);
643 Resid(LID, m) = X(LID, m);
644 for(
size_t k = 0; k < NumEntries; k++)
646 local_ordinal_type col = Indices[k];
647 Resid(LID, m) -= Values[k] * Y2(col, m);
658 apply(Resid, Y2, i, stride, Teuchos::NO_TRANS, DampingFactor_, one);
663 local_ordinal_type LRID = partitions_[partitionIndices_[i]];
665 HostView diagView = Diag_->template getLocalView<Kokkos::HostSpace>();
666 impl_scalar_type d = (impl_scalar_type) one / diagView(LRID, 0);
667 for(
size_t m = 0; m < numVecs; m++)
669 impl_scalar_type x = X(LRID, m);
670 impl_scalar_type newy = x * d;
682 for(local_ordinal_type i = numBlocks_; i > 0; --i)
684 if(hasBlockCrs_ || blockRows_[i-1] != 1)
686 if(blockRows_[i - 1] == 0)
689 ArrayView<const local_ordinal_type> localRows = getLocalRows(i-1);
690 for(local_ordinal_type j = 0; j < blockRows_[i-1]; j++)
692 const local_ordinal_type LID = localRows[j];
694 inputMatrix_->getLocalRowCopy(LID, Indices(), Values(), NumEntries);
696 for (
size_t m = 0; m < numVecs; m++)
700 for(
int localR = 0; localR < bcrsBlockSize_; localR++)
701 Resid(LID * bcrsBlockSize_ + localR, m) = X(LID * bcrsBlockSize_ + localR, m);
702 for(
size_t k = 0; k < NumEntries; ++k)
704 const local_ordinal_type col = Indices[k];
705 for(
int localR = 0; localR < bcrsBlockSize_; localR++)
707 for(
int localC = 0; localC < bcrsBlockSize_; localC++)
708 Resid(LID*bcrsBlockSize_+localR, m) -=
709 Values[k * (bcrsBlockSize_ * bcrsBlockSize_) + (localR + localC * bcrsBlockSize_)]
710 * Y2(col * bcrsBlockSize_ + localC, m);
716 Resid(LID, m) = X(LID, m);
717 for(
size_t k = 0; k < NumEntries; ++k)
719 local_ordinal_type col = Indices[k];
720 Resid(LID, m) -= Values[k] * Y2(col, m);
731 apply(Resid, Y2, i - 1, stride, Teuchos::NO_TRANS, DampingFactor_, one);
738 auto numMyRows = inputMatrix_->getNodeNumRows();
739 for (
size_t m = 0; m < numVecs; ++m)
741 for (
size_t i = 0; i < numMyRows * bcrsBlockSize_; ++i)
749 template<
class MatrixType>
750 void Container<MatrixType>::clearBlocks()
755 partitionIndices_.clear();
756 Diag_ = Teuchos::null;
767 template<
class MatrixType>
771 static std::string name () {
772 return std::string (
"Ifpack2::Container<") +
773 TypeNameTraits<MatrixType>::name () +
" >";
776 static std::string concreteName (const ::Ifpack2::Container<MatrixType>&) {
783 #endif // IFPACK2_CONTAINER_HPP int OverlapLevel_
Number of rows of overlap for adjacent blocks.
Definition: Ifpack2_Container.hpp:409
global_ordinal_type NumGlobalNonzeros_
Number of nonzeros in input matrix.
Definition: Ifpack2_Container.hpp:419
virtual void initialize()=0
Do all set-up operations that only require matrix structure.
int numBlocks_
The number of blocks (partitions) in the container.
Definition: Ifpack2_Container.hpp:397
Teuchos::Array< local_ordinal_type > blockRows_
Number of rows in each block.
Definition: Ifpack2_Container.hpp:401
static std::string getName()
Definition: Ifpack2_Container.hpp:387
global_ordinal_type NumGlobalRows_
Number of global rows in input matrix.
Definition: Ifpack2_Container.hpp:417
bool hasBlockCrs_
Whether the input matrix is a BlockCRS matrix.
Definition: Ifpack2_Container.hpp:421
scalar_type DampingFactor_
Damping factor, passed to apply() as alpha.
Definition: Ifpack2_Container.hpp:411
Ifpack2::Partitioner:
Definition: Ifpack2_Partitioner.hpp:179
int bcrsBlockSize_
If hasBlockCrs_, the number of DOFs per vertex. Otherwise 1.
Definition: Ifpack2_Container.hpp:423
virtual bool isComputed() const =0
Return true if the container has been successfully computed.
local_ordinal_type NumLocalRows_
Number of local rows in input matrix.
Definition: Ifpack2_Container.hpp:415
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 =0
Compute Y := alpha * diag(D) * M^{-1} (diag(D) * X) + beta*Y.
Teuchos::RCP< const row_matrix_type > inputMatrix_
The input matrix to the constructor.
Definition: Ifpack2_Container.hpp:394
Teuchos::RCP< vector_type > Diag_
Diagonal elements.
Definition: Ifpack2_Container.hpp:405
Container(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_Container.hpp:149
Teuchos::RCP< const Tpetra::Import< local_ordinal_type, global_ordinal_type, node_type > > Importer_
Importer for importing off-process elements of MultiVectors.
Definition: Ifpack2_Container.hpp:413
bool IsParallel_
Whether the problem is distributed across multiple MPI processes.
Definition: Ifpack2_Container.hpp:407
void weightedApplyMV(mv_type &X, mv_type &Y, vector_type &W)
Wrapper for weightedApply with MVs, used in unit tests (never called by BlockRelaxation) ...
Definition: Ifpack2_Container.hpp:370
Definition: Ifpack2_Container.hpp:761
Teuchos::ArrayView< const local_ordinal_type > getLocalRows(int blockIndex) const
Local indices of the rows of the input matrix that belong to this block.
Definition: Ifpack2_Container.hpp:236
Teuchos::Array< local_ordinal_type > partitionIndices_
Starting index in partitions_ of local row indices for each block.
Definition: Ifpack2_Container.hpp:403
virtual void compute()=0
Extract the local diagonal block and prepare the solver.
Interface for creating and solving a local linear problem.
Definition: Ifpack2_Container.hpp:114
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 =0
Compute Y := alpha * M^{-1} X + beta*Y.
virtual void setParameters(const Teuchos::ParameterList &List)=0
Set parameters.
virtual std::ostream & print(std::ostream &os) const =0
Print basic information about the container to os.
virtual ~Container()
Destructor.
Definition: Ifpack2_Container.hpp:219
Teuchos::Array< local_ordinal_type > partitions_
Local indices of the rows of the input matrix that belong to this block.
Definition: Ifpack2_Container.hpp:399
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
virtual bool isInitialized() const =0
Return true if the container has been successfully initialized.
void setBlockSizes(const Teuchos::Array< Teuchos::Array< local_ordinal_type > > &partitions)
Initialize arrays with information about block sizes.
Definition: Ifpack2_Container.hpp:253
Kokkos::Details::ArithTraits< scalar_type >::val_type impl_scalar_type
Internal representation of Scalar in Kokkos::View.
Definition: Ifpack2_Container.hpp:132
Container(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< local_ordinal_type > &localRows)
Constructor for single block.
Definition: Ifpack2_Container.hpp:190
void applyMV(mv_type &X, mv_type &Y) const
Wrapper for apply with MVs, used in unit tests (never called by BlockRelaxation)
Definition: Ifpack2_Container.hpp:332