10 #ifndef IFPACK2_SINGLETONFILTER_DEF_HPP
11 #define IFPACK2_SINGLETONFILTER_DEF_HPP
12 #include "Ifpack2_SingletonFilter_decl.hpp"
14 #include "Tpetra_ConfigDefs.hpp"
15 #include "Tpetra_RowMatrix.hpp"
16 #include "Tpetra_Map.hpp"
17 #include "Tpetra_MultiVector.hpp"
18 #include "Tpetra_Vector.hpp"
22 template <
class MatrixType>
29 , MaxNumEntriesA_(0) {
31 if (A_->getComm()->getSize() != 1 || A_->getLocalNumRows() != A_->getGlobalNumRows()) {
32 throw std::runtime_error(
"Ifpack2::SingeltonFilter can be used with Comm().getSize() == 1 only. This class is a tool for Ifpack2_AdditiveSchwarz, and it is not meant to be used otherwise.");
36 size_t NumRowsA_ = A_->getLocalNumRows();
40 MaxNumEntriesA_ = A_->getLocalMaxNumRowEntries();
43 Kokkos::resize(Indices_, MaxNumEntriesA_);
44 Kokkos::resize(Values_, MaxNumEntriesA_);
47 Reorder_.resize(NumRowsA_);
48 Reorder_.assign(Reorder_.size(), -1);
52 for (
size_t i = 0; i < NumRowsA_; ++i) {
54 A_->getLocalRowCopy(i, Indices_, Values_, Nnz);
56 Reorder_[i] = NumRows_++;
63 InvReorder_.resize(NumRows_);
64 for (
size_t i = 0; i < NumRowsA_; ++i) {
67 InvReorder_[Reorder_[i]] = i;
69 NumEntries_.resize(NumRows_);
70 SingletonIndex_.resize(NumSingletons_);
74 for (
size_t i = 0; i < NumRowsA_; ++i) {
76 A_->getLocalRowCopy(i, Indices_, Values_, Nnz);
77 LocalOrdinal ii = Reorder_[i];
79 NumEntries_[ii] = Nnz;
81 if (Nnz > MaxNumEntries_)
84 SingletonIndex_[count] = i;
90 ReducedMap_ =
Teuchos::rcp(
new Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>(NumRows_, 0, A_->getComm()));
93 Diagonal_ =
Teuchos::rcp(
new Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(ReducedMap_));
95 Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> DiagonalA(A_->getRowMap());
96 A_->getLocalDiagCopy(DiagonalA);
98 for (
size_t i = 0; i < NumRows_; ++i) {
99 LocalOrdinal ii = InvReorder_[i];
100 Diagonal_->replaceLocalValue((LocalOrdinal)i, DiagonalAview[ii]);
104 template <
class MatrixType>
107 template <
class MatrixType>
110 return A_->getComm();
113 template <
class MatrixType>
114 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
115 typename MatrixType::global_ordinal_type,
116 typename MatrixType::node_type> >
121 template <
class MatrixType>
122 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
123 typename MatrixType::global_ordinal_type,
124 typename MatrixType::node_type> >
129 template <
class MatrixType>
130 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
131 typename MatrixType::global_ordinal_type,
132 typename MatrixType::node_type> >
137 template <
class MatrixType>
138 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
139 typename MatrixType::global_ordinal_type,
140 typename MatrixType::node_type> >
145 template <
class MatrixType>
146 Teuchos::RCP<
const Tpetra::RowGraph<
typename MatrixType::local_ordinal_type,
147 typename MatrixType::global_ordinal_type,
148 typename MatrixType::node_type> >
150 throw std::runtime_error(
"Ifpack2::SingletonFilter: does not support getGraph.");
153 template <
class MatrixType>
158 template <
class MatrixType>
163 template <
class MatrixType>
168 template <
class MatrixType>
173 template <
class MatrixType>
175 return A_->getIndexBase();
178 template <
class MatrixType>
183 template <
class MatrixType>
188 template <
class MatrixType>
190 throw std::runtime_error(
"Ifpack2::SingletonFilter does not implement getNumEntriesInGlobalRow.");
193 template <
class MatrixType>
195 return NumEntries_[localRow];
198 template <
class MatrixType>
200 return MaxNumEntries_;
203 template <
class MatrixType>
205 return MaxNumEntries_;
208 template <
class MatrixType>
210 return A_->getBlockSize();
213 template <
class MatrixType>
218 template <
class MatrixType>
220 return A_->isLocallyIndexed();
223 template <
class MatrixType>
225 return A_->isGloballyIndexed();
228 template <
class MatrixType>
230 return A_->isFillComplete();
233 template <
class MatrixType>
236 nonconst_global_inds_host_view_type& ,
237 nonconst_values_host_view_type& ,
239 throw std::runtime_error(
"Ifpack2::SingletonFilter does not implement getGlobalRowCopy.");
242 template <
class MatrixType>
245 nonconst_local_inds_host_view_type& Indices,
246 nonconst_values_host_view_type& Values,
247 size_t& NumEntries)
const {
248 TEUCHOS_TEST_FOR_EXCEPTION((LocalRow < 0 || (
size_t)LocalRow >= NumRows_ || (
size_t)Indices.size() < NumEntries_[LocalRow]), std::runtime_error,
"Ifpack2::SingletonFilter::getLocalRowCopy invalid row or array size.");
251 LocalOrdinal ARow = InvReorder_[LocalRow];
252 A_->getLocalRowCopy(ARow, Indices_, Values_, Nnz);
256 for (
size_t i = 0; i < Nnz; ++i) {
257 LocalOrdinal ii = Reorder_[Indices_[i]];
259 Indices[NumEntries] = ii;
260 Values[NumEntries] = Values_[i];
266 template <
class MatrixType>
268 global_inds_host_view_type& ,
269 values_host_view_type& )
const {
270 throw std::runtime_error(
"Ifpack2::SingletonFilter: does not support getGlobalRowView.");
273 template <
class MatrixType>
275 local_inds_host_view_type& ,
276 values_host_view_type& )
const {
277 throw std::runtime_error(
"Ifpack2::SingletonFilter: does not support getLocalRowView.");
280 template <
class MatrixType>
283 return A_->getLocalDiagCopy(diag);
286 template <
class MatrixType>
288 throw std::runtime_error(
"Ifpack2::SingletonFilter does not support leftScale.");
291 template <
class MatrixType>
293 throw std::runtime_error(
"Ifpack2::SingletonFilter does not support rightScale.");
296 template <
class MatrixType>
298 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Y,
302 typedef Scalar DomainScalar;
303 typedef Scalar RangeScalar;
308 "Ifpack2::SingletonFilter::apply ERROR: X.getNumVectors() != Y.getNumVectors().");
315 size_t NumVectors = Y.getNumVectors();
317 for (
size_t i = 0; i < NumRows_; ++i) {
320 getLocalRowCopy(i, Indices_, Values_, Nnz);
322 for (
size_t j = 0; j < Nnz; ++j)
323 for (
size_t k = 0; k < NumVectors; ++k)
324 y_ptr[k][i] += (RangeScalar)Values_[j] * (RangeScalar)x_ptr[k][Indices_[j]];
326 for (
size_t j = 0; j < Nnz; ++j)
327 for (
size_t k = 0; k < NumVectors; ++k)
328 y_ptr[k][Indices_[j]] += (RangeScalar)Values_[j] * (RangeScalar)x_ptr[k][i];
330 for (
size_t j = 0; j < Nnz; ++j)
331 for (
size_t k = 0; k < NumVectors; ++k)
337 template <
class MatrixType>
342 template <
class MatrixType>
347 template <
class MatrixType>
349 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& LHS) {
350 this->
template SolveSingletonsTempl<Scalar, Scalar>(RHS, LHS);
353 template <
class MatrixType>
354 template <
class DomainScalar,
class RangeScalar>
356 Tpetra::MultiVector<RangeScalar, LocalOrdinal, GlobalOrdinal, Node>& LHS) {
360 for (
size_t i = 0; i < NumSingletons_; ++i) {
361 LocalOrdinal ii = SingletonIndex_[i];
364 A_->getLocalRowCopy(ii, Indices_, Values_, Nnz);
365 for (
size_t j = 0; j < Nnz; ++j) {
366 if (Indices_[j] == ii) {
367 for (
size_t k = 0; k < LHS.getNumVectors(); ++k)
368 LHS_ptr[k][ii] = (RangeScalar)RHS_ptr[k][ii] / (RangeScalar)Values_[j];
374 template <
class MatrixType>
376 const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& RHS,
377 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& ReducedRHS) {
378 this->
template CreateReducedRHSTempl<Scalar, Scalar>(LHS, RHS, ReducedRHS);
381 template <
class MatrixType>
382 template <
class DomainScalar,
class RangeScalar>
384 const Tpetra::MultiVector<RangeScalar, LocalOrdinal, GlobalOrdinal, Node>& RHS,
385 Tpetra::MultiVector<RangeScalar, LocalOrdinal, GlobalOrdinal, Node>& ReducedRHS) {
390 size_t NumVectors = LHS.getNumVectors();
392 for (
size_t i = 0; i < NumRows_; ++i)
393 for (
size_t k = 0; k < NumVectors; ++k)
394 ReducedRHS_ptr[k][i] = RHS_ptr[k][InvReorder_[i]];
396 for (
size_t i = 0; i < NumRows_; ++i) {
397 LocalOrdinal ii = InvReorder_[i];
399 A_->getLocalRowCopy(ii, Indices_, Values_, Nnz);
401 for (
size_t j = 0; j < Nnz; ++j) {
402 if (Reorder_[Indices_[j]] == -1) {
403 for (
size_t k = 0; k < NumVectors; ++k)
404 ReducedRHS_ptr[k][i] -= (RangeScalar)Values_[j] * (RangeScalar)LHS_ptr[k][Indices_[j]];
410 template <
class MatrixType>
412 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& LHS) {
413 this->
template UpdateLHSTempl<Scalar, Scalar>(ReducedLHS, LHS);
416 template <
class MatrixType>
417 template <
class DomainScalar,
class RangeScalar>
419 Tpetra::MultiVector<RangeScalar, LocalOrdinal, GlobalOrdinal, Node>& LHS) {
423 for (
size_t i = 0; i < NumRows_; ++i)
424 for (
size_t k = 0; k < LHS.getNumVectors(); ++k)
425 LHS_ptr[k][InvReorder_[i]] = (RangeScalar)ReducedLHS_ptr[k][i];
428 template <
class MatrixType>
430 throw std::runtime_error(
"Ifpack2::SingletonFilter does not implement getFrobeniusNorm.");
435 #define IFPACK2_SINGLETONFILTER_INSTANT(S, LO, GO, N) \
436 template class Ifpack2::SingletonFilter<Tpetra::RowMatrix<S, LO, GO, N> >;
virtual bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
Definition: Ifpack2_SingletonFilter_def.hpp:338
virtual mag_type getFrobeniusNorm() const
Returns the Frobenius norm of the matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:429
virtual void rightScale(const Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &x)
Scales the RowMatrix on the right with the Vector x.
Definition: Ifpack2_SingletonFilter_def.hpp:292
virtual ~SingletonFilter()
Destructor.
Definition: Ifpack2_SingletonFilter_def.hpp:105
virtual Teuchos::RCP< const Tpetra::RowGraph< LocalOrdinal, GlobalOrdinal, Node > > getGraph() const
Returns the RowGraph associated with this matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:149
virtual void getLocalRowCopy(LocalOrdinal LocalRow, nonconst_local_inds_host_view_type &Indices, nonconst_values_host_view_type &Values, size_t &NumEntries) const
Extract a list of entries in a specified local row of the graph. Put into storage allocated by callin...
Definition: Ifpack2_SingletonFilter_def.hpp:244
virtual void UpdateLHS(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &ReducedLHS, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &LHS)
Updates a full LHS from a reduces LHS.
Definition: Ifpack2_SingletonFilter_def.hpp:411
virtual LocalOrdinal getBlockSize() const
The number of degrees of freedom per mesh point.
Definition: Ifpack2_SingletonFilter_def.hpp:209
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:117
virtual void CreateReducedRHS(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &LHS, const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &ReducedRHS)
Creates a RHS for the reduced singleton-free system.
Definition: Ifpack2_SingletonFilter_def.hpp:375
virtual size_t getLocalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on this node.
Definition: Ifpack2_SingletonFilter_def.hpp:204
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
Returns the Map that describes the domain distribution in this matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:133
virtual bool supportsRowViews() const
Returns true if RowViews are supported.
Definition: Ifpack2_SingletonFilter_def.hpp:343
virtual global_size_t getGlobalNumRows() const
Returns the number of global rows in this matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:154
virtual void getGlobalRowCopy(GlobalOrdinal GlobalRow, nonconst_global_inds_host_view_type &Indices, nonconst_values_host_view_type &Values, size_t &NumEntries) const
Extract a list of entries in a specified global row of this matrix. Put into pre-allocated storage...
Definition: Ifpack2_SingletonFilter_def.hpp:235
virtual GlobalOrdinal getIndexBase() const
Returns the index base for global indices for this matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:174
virtual bool isFillComplete() const
Returns true if fillComplete() has been called.
Definition: Ifpack2_SingletonFilter_def.hpp:229
virtual global_size_t getGlobalNumCols() const
Returns the number of global columns in this matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:159
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
Returns the Map that describes the range distribution in this matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:141
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual bool isGloballyIndexed() const
If matrix indices are in the global range, this function returns true. Otherwise, this function retur...
Definition: Ifpack2_SingletonFilter_def.hpp:224
virtual global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:179
virtual bool hasColMap() const
Indicates whether this matrix has a well-defined column map.
Definition: Ifpack2_SingletonFilter_def.hpp:214
virtual bool isLocallyIndexed() const
If matrix indices are in the local range, this function returns true. Otherwise, this function return...
Definition: Ifpack2_SingletonFilter_def.hpp:219
virtual void getLocalRowView(LocalOrdinal LocalRow, local_inds_host_view_type &indices, values_host_view_type &values) const
Extract a const, non-persisting view of local indices in a specified row of the matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:274
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
Definition: Ifpack2_SingletonFilter_def.hpp:194
virtual size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
Returns the current number of entries on this node in the specified global row.
Definition: Ifpack2_SingletonFilter_def.hpp:189
virtual size_t getLocalNumCols() const
Returns the number of columns needed to apply the forward operator on this node, i.e., the number of elements listed in the column map.
Definition: Ifpack2_SingletonFilter_def.hpp:169
virtual size_t getLocalNumEntries() const
Returns the local number of entries in this matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:184
virtual size_t getLocalNumRows() const
Returns the number of rows owned on the calling node.
Definition: Ifpack2_SingletonFilter_def.hpp:164
virtual size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
Definition: Ifpack2_SingletonFilter_def.hpp:199
virtual void getLocalDiagCopy(Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &diag) const
Get a copy of the diagonal entries owned by this node, with local row indices.
Definition: Ifpack2_SingletonFilter_def.hpp:281
virtual void getGlobalRowView(GlobalOrdinal GlobalRow, global_inds_host_view_type &indices, values_host_view_type &values) const
Extract a const, non-persisting view of global indices in a specified row of the matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:267
SingletonFilter(const Teuchos::RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &Matrix)
Constructor.
Definition: Ifpack2_SingletonFilter_def.hpp:23
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the communicator.
Definition: Ifpack2_SingletonFilter_def.hpp:109
virtual void apply(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const
Computes the operator-multivector application.
Definition: Ifpack2_SingletonFilter_def.hpp:297
virtual void SolveSingletons(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &LHS)
Solve the singleton components of the linear system.
Definition: Ifpack2_SingletonFilter_def.hpp:348
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getColMap() const
Returns the Map that describes the column distribution in this matrix.
Definition: Ifpack2_SingletonFilter_def.hpp:125
Filter based on matrix entries.
Definition: Ifpack2_SingletonFilter_decl.hpp:30
virtual void leftScale(const Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &x)
Scales the RowMatrix on the left with the Vector x.
Definition: Ifpack2_SingletonFilter_def.hpp:287