Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_SparsityFilter_def.hpp
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#ifndef IFPACK2_SPARSITYFILTER_DEF_HPP
44#define IFPACK2_SPARSITYFILTER_DEF_HPP
45
46#include "Tpetra_Map.hpp"
47#include "Tpetra_MultiVector.hpp"
48#include "Tpetra_Vector.hpp"
49
50#include <algorithm>
51#include <vector>
52
53namespace Ifpack2 {
54
55//==========================================================================
56template<class MatrixType>
57SparsityFilter<MatrixType>::SparsityFilter(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Matrix,
58 size_t AllowedNumEntries,
59 LocalOrdinal AllowedBandwidth):
60 A_(Matrix),
61 AllowedNumEntries_(AllowedNumEntries),
62 AllowedBandwidth_(AllowedBandwidth),
63 NumRows_(0),
64 NumNonzeros_(0),
65 MaxNumEntries_(0),
66 MaxNumEntriesA_(0)
67{
68
69 // use this filter only on serial matrices
70 TEUCHOS_TEST_FOR_EXCEPTION(
71 A_->getComm()->getSize() != 1 || A_->getLocalNumRows() != A_->getGlobalNumRows(),
72 std::runtime_error, "Ifpack2::SparsityFilter: "
73 "This class may only be used when A.getComm()->getSize() == 1.");
74
75 // localized matrix has all the local rows of Matrix
76 NumRows_ = A_->getLocalNumRows();
77
78 // NodeNumEntries_ will contain the actual number of nonzeros
79 // for each localized row (that is, without external nodes,
80 // and always with the diagonal entry)
81 NumEntries_.resize(NumRows_);
82
83 // tentative value for MaxNumEntries. This is the number of
84 // nonzeros in the local matrix
85 MaxNumEntries_ = A_->getLocalMaxNumRowEntries();
86 MaxNumEntriesA_ = A_->getLocalMaxNumRowEntries();
87
88 // ExtractMyRowCopy() will use these vectors
89 Kokkos::resize(Indices_,MaxNumEntries_);
90 Kokkos::resize(Values_,MaxNumEntries_);
91
92 size_t ActualMaxNumEntries = 0;
93 for (size_t i = 0 ; i < NumRows_ ; ++i) {
94 NumEntries_[i] = MaxNumEntriesA_;
95 size_t Nnz;
96 A_->getLocalRowCopy(i,Indices_,Values_,Nnz);
97
98 NumNonzeros_ += Nnz;
99 NumEntries_[i] = Nnz;
100 if (Nnz > ActualMaxNumEntries)
101 ActualMaxNumEntries = Nnz;
102 }
103
104 MaxNumEntries_ = ActualMaxNumEntries;
105}
106
107//=========================================================================
108template<class MatrixType>
110
111//==========================================================================
112template<class MatrixType>
113Teuchos::RCP<const Teuchos::Comm<int> > SparsityFilter<MatrixType>::getComm() const
114{
115 return A_->getComm();
116}
117
118
119//==========================================================================
120template<class MatrixType>
121Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
122 typename MatrixType::global_ordinal_type,
123 typename MatrixType::node_type> >
125{
126 return A_->getRowMap();
127}
128
129//==========================================================================
130template<class MatrixType>
131Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
132 typename MatrixType::global_ordinal_type,
133 typename MatrixType::node_type> >
135{
136 return A_->getColMap();
137}
138
139//==========================================================================
140template<class MatrixType>
141Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
142 typename MatrixType::global_ordinal_type,
143 typename MatrixType::node_type> >
145{
146 return A_->getDomainMap();
147}
148
149//==========================================================================
150template<class MatrixType>
151Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
152 typename MatrixType::global_ordinal_type,
153 typename MatrixType::node_type> >
155{
156 return A_->getRangeMap();
157}
158
159//==========================================================================
160template<class MatrixType>
161Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type,
162 typename MatrixType::global_ordinal_type,
163 typename MatrixType::node_type> >
165{
166 throw std::runtime_error("Ifpack2::SparsityFilter: does not support getGraph.");
167}
168
169//==========================================================================
170template<class MatrixType>
172{
173 return NumRows_;
174}
175
176//==========================================================================
177template<class MatrixType>
179{
180 return NumRows_;
181}
182
183//==========================================================================
184template<class MatrixType>
186{
187 return NumRows_;
188}
189
190//==========================================================================
191
192template<class MatrixType>
194{
195 return NumRows_;
196}
197
198//==========================================================================
199template<class MatrixType>
200typename MatrixType::global_ordinal_type SparsityFilter<MatrixType>::getIndexBase() const
201{
202 return A_->getIndexBase();
203}
204
205//==========================================================================
206template<class MatrixType>
208{
209 return NumNonzeros_;
210}
211
212//==========================================================================
213template<class MatrixType>
215{
216 return NumNonzeros_;
217}
218
219//==========================================================================
220template<class MatrixType>
221size_t SparsityFilter<MatrixType>::getNumEntriesInGlobalRow(GlobalOrdinal /* globalRow */) const
222{
223 throw std::runtime_error("Ifpack2::SparsityFilter does not implement getNumEntriesInGlobalRow.");
224}
225
226//==========================================================================
227template<class MatrixType>
229{
230 return NumEntries_[localRow];
231}
232
233//==========================================================================
234template<class MatrixType>
236{
237 return MaxNumEntries_;
238}
239
240//==========================================================================
241template<class MatrixType>
243{
244 return MaxNumEntries_;
245}
246
247//==========================================================================
248template<class MatrixType>
250{
251 return true;
252}
253
254//==========================================================================
255template<class MatrixType>
257{
258 return A_->isLocallyIndexed();
259}
260
261//==========================================================================
262template<class MatrixType>
264{
265 return A_->isGloballyIndexed();
266}
267
268//==========================================================================
269template<class MatrixType>
271{
272 return A_->isFillComplete();
273}
274
275//==========================================================================
276template<class MatrixType>
278getGlobalRowCopy (GlobalOrdinal /*GlobalRow*/,
279 nonconst_global_inds_host_view_type &/*Indices*/,
280 nonconst_values_host_view_type &/*Values*/,
281 size_t& /*NumEntries*/) const {
282 throw std::runtime_error("Ifpack2::SparsityFilter does not implement getGlobalRowCopy.");
283}
284
285//==========================================================================
286template<class MatrixType>
288 getLocalRowCopy (LocalOrdinal LocalRow,
289 nonconst_local_inds_host_view_type &Indices,
290 nonconst_values_host_view_type &Values,
291 size_t& NumEntries) const
292{
293TEUCHOS_TEST_FOR_EXCEPTION((LocalRow < 0 || (size_t) LocalRow >= NumRows_ || (size_t) Indices.size() < NumEntries_[LocalRow]), std::runtime_error, "Ifpack2::SparsityFilter::getLocalRowCopy invalid row or array size.");
294
295 // Note: This function will work correctly if called by apply, say, with Indices_ and Values_ as
296 // parameters. The structure of the loop below should make that obvious.
297
298 // always extract using the object Values_ and Indices_.
299 // This is because I need more space than that given by
300 // the user (for the external nodes)
301 size_t A_NumEntries=0;
302 A_->getLocalRowCopy(LocalRow,Indices_,Values_,A_NumEntries);
303 magnitudeType Threshold = Teuchos::ScalarTraits<magnitudeType>::zero();
304 std::vector<magnitudeType> Values2(A_NumEntries,Teuchos::ScalarTraits<magnitudeType>::zero());
305
306 if (A_NumEntries > AllowedNumEntries_) {
307 size_t count = 0;
308 for (size_t i = 0 ; i < A_NumEntries ; ++i) {
309 // skip diagonal entry (which is always inserted)
310 if (Indices_[i] == LocalRow)
311 continue;
312 // put magnitude
313 Values2[count] = Teuchos::ScalarTraits<Scalar>::magnitude(Values_[i]);
314 count++;
315 }
316
317 // sort in descending order
318 std::sort(Values2.rbegin(),Values2.rend());
319 // get the cut-off value
320 Threshold = Values2[AllowedNumEntries_ - 1];
321 }
322
323
324 // loop over all nonzero elements of row MyRow,
325 // and drop elements below specified threshold.
326 // Also, add value to zero diagonal
327 NumEntries = 0;
328 for (size_t i = 0 ; i < A_NumEntries ; ++i) {
329 if (std::abs(Indices_[i] - LocalRow) > AllowedBandwidth_)
330 continue;
331
332 if ((Indices_[i] != LocalRow) && (Teuchos::ScalarTraits<Scalar>::magnitude(Values_[i]) < Threshold))
333 continue;
334
335 Values[NumEntries] = Values_[i];
336 Indices[NumEntries] = Indices_[i];
337
338 NumEntries++;
339 if (NumEntries > AllowedNumEntries_)
340 break;
341 }
342
343
344}
345
346//==========================================================================
347template<class MatrixType>
348void SparsityFilter<MatrixType>::getGlobalRowView(GlobalOrdinal /* GlobalRow */,
349 global_inds_host_view_type &/*indices*/,
350 values_host_view_type &/*values*/) const
351{
352 throw std::runtime_error("Ifpack2::SparsityFilter: does not support getGlobalRowView.");
353}
354
355//==========================================================================
356template<class MatrixType>
357void SparsityFilter<MatrixType>::getLocalRowView(LocalOrdinal /* LocalRow */,
358 local_inds_host_view_type & /*indices*/,
359 values_host_view_type & /*values*/) const
360{
361 throw std::runtime_error("Ifpack2::SparsityFilter: does not support getLocalRowView.");
362}
363
364//==========================================================================
365template<class MatrixType>
366void SparsityFilter<MatrixType>::getLocalDiagCopy(Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const
367{
368 // This is somewhat dubious as to how the maps match.
369 return A_->getLocalDiagCopy(diag);
370}
371
372//==========================================================================
373template<class MatrixType>
374void SparsityFilter<MatrixType>::leftScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* x */)
375{
376 throw std::runtime_error("Ifpack2::SparsityFilter does not support leftScale.");
377}
378
379//==========================================================================
380template<class MatrixType>
381void SparsityFilter<MatrixType>::rightScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* x */)
382{
383 throw std::runtime_error("Ifpack2::SparsityFilter does not support rightScale.");
384}
385
386//==========================================================================
387template<class MatrixType>
388void SparsityFilter<MatrixType>::apply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
389 Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
390 Teuchos::ETransp mode,
391 Scalar /* alpha */,
392 Scalar /* beta */) const
393{
394 // Note: This isn't AztecOO compliant. But neither was Ifpack's version.
395 // Note: The localized maps mean the matvec is trivial (and has no import)
396 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
397 "Ifpack2::SparsityFilter::apply ERROR: X.getNumVectors() != Y.getNumVectors().");
398
399 Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
400 Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > x_ptr = X.get2dView();
401 Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > y_ptr = Y.get2dViewNonConst();
402
403 Y.putScalar(zero);
404 size_t NumVectors = Y.getNumVectors();
405
406 for (size_t i = 0 ; i < NumRows_ ; ++i) {
407 size_t Nnz;
408 // Use this class's getrow to make the below code simpler
409 getLocalRowCopy(i,Indices_,Values_,Nnz);
410 Scalar* Values = reinterpret_cast<Scalar*>(Values_.data());
411 if (mode==Teuchos::NO_TRANS){
412 for (size_t j = 0 ; j < Nnz ; ++j)
413 for (size_t k = 0 ; k < NumVectors ; ++k)
414 y_ptr[k][i] += Values[j] * x_ptr[k][Indices_[j]];
415 }
416 else if (mode==Teuchos::TRANS){
417 for (size_t j = 0 ; j < Nnz ; ++j)
418 for (size_t k = 0 ; k < NumVectors ; ++k)
419 y_ptr[k][Indices_[j]] += Values[j] * x_ptr[k][i];
420 }
421 else { //mode==Teuchos::CONJ_TRANS
422 for (size_t j = 0 ; j < Nnz ; ++j)
423 for (size_t k = 0 ; k < NumVectors ; ++k)
424 y_ptr[k][Indices_[j]] += Teuchos::ScalarTraits<Scalar>::conjugate(Values[j]) * x_ptr[k][i];
425 }
426 }
427}
428
429
430//==========================================================================
431template<class MatrixType>
433{
434 return true;
435}
436
437//==========================================================================
438template<class MatrixType>
440{
441 return false;
442}
443
444//==========================================================================
445template<class MatrixType>
446typename SparsityFilter<MatrixType>::mag_type SparsityFilter<MatrixType>::getFrobeniusNorm() const
447{
448 throw std::runtime_error("Ifpack2::SparsityFilter does not implement getFrobeniusNorm.");
449}
450
451} // namespace Ifpack2
452
453#define IFPACK2_SPARSITYFILTER_INSTANT(S,LO,GO,N) \
454 template class Ifpack2::SparsityFilter< Tpetra::RowMatrix<S, LO, GO, N> >;
455
456#endif
virtual global_size_t getGlobalNumRows() const
Returns the number of global rows in this matrix.
Definition Ifpack2_SparsityFilter_def.hpp:171
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
Definition Ifpack2_SparsityFilter_def.hpp:228
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_SparsityFilter_def.hpp:357
virtual bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
Definition Ifpack2_SparsityFilter_def.hpp:432
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
Returns the Map that describes the domain distribution in this matrix.
Definition Ifpack2_SparsityFilter_def.hpp:144
virtual GlobalOrdinal getIndexBase() const
Returns the index base for global indices for this matrix.
Definition Ifpack2_SparsityFilter_def.hpp:200
virtual mag_type getFrobeniusNorm() const
Returns the Frobenius norm of the matrix.
Definition Ifpack2_SparsityFilter_def.hpp:446
virtual global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
Definition Ifpack2_SparsityFilter_def.hpp:207
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the communicator.
Definition Ifpack2_SparsityFilter_def.hpp:113
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition Ifpack2_SparsityFilter_def.hpp:124
virtual size_t getLocalNumCols() const
Returns the number of columns needed to apply the forward operator on this node, i....
Definition Ifpack2_SparsityFilter_def.hpp:193
virtual size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
Returns the current number of entries on this node in the specified global row.
Definition Ifpack2_SparsityFilter_def.hpp:221
virtual size_t getLocalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on this node.
Definition Ifpack2_SparsityFilter_def.hpp:242
virtual ~SparsityFilter()
Destructor.
Definition Ifpack2_SparsityFilter_def.hpp:109
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_SparsityFilter_def.hpp:348
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_SparsityFilter_def.hpp:366
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_SparsityFilter_def.hpp:278
virtual bool isLocallyIndexed() const
If matrix indices are in the local range, this function returns true. Otherwise, this function return...
Definition Ifpack2_SparsityFilter_def.hpp:256
virtual global_size_t getGlobalNumCols() const
Returns the number of global columns in this matrix.
Definition Ifpack2_SparsityFilter_def.hpp:178
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
Returns the Map that describes the range distribution in this matrix.
Definition Ifpack2_SparsityFilter_def.hpp:154
virtual bool hasColMap() const
Indicates whether this matrix has a well-defined column map.
Definition Ifpack2_SparsityFilter_def.hpp:249
virtual void rightScale(const Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &x)
Scales the RowMatrix on the right with the Vector x.
Definition Ifpack2_SparsityFilter_def.hpp:381
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_SparsityFilter_def.hpp:388
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getColMap() const
Returns the Map that describes the column distribution in this matrix.
Definition Ifpack2_SparsityFilter_def.hpp:134
virtual bool supportsRowViews() const
Returns true if RowViews are supported.
Definition Ifpack2_SparsityFilter_def.hpp:439
virtual bool isFillComplete() const
Returns true if fillComplete() has been called.
Definition Ifpack2_SparsityFilter_def.hpp:270
virtual void leftScale(const Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &x)
Scales the RowMatrix on the left with the Vector x.
Definition Ifpack2_SparsityFilter_def.hpp:374
virtual Teuchos::RCP< const Tpetra::RowGraph< LocalOrdinal, GlobalOrdinal, Node > > getGraph() const
Returns the RowGraph associated with this matrix.
Definition Ifpack2_SparsityFilter_def.hpp:164
virtual size_t getLocalNumRows() const
Returns the number of rows owned on the calling node.
Definition Ifpack2_SparsityFilter_def.hpp:185
SparsityFilter(const Teuchos::RCP< const row_matrix_type > &Matrix, size_t AllowedNumEntries, LocalOrdinal AllowedBandwidth=-Teuchos::ScalarTraits< LocalOrdinal >::one())
Constructor.
Definition Ifpack2_SparsityFilter_def.hpp:57
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_SparsityFilter_def.hpp:288
virtual size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
Definition Ifpack2_SparsityFilter_def.hpp:235
virtual size_t getLocalNumEntries() const
Returns the local number of entries in this matrix.
Definition Ifpack2_SparsityFilter_def.hpp:214
virtual bool isGloballyIndexed() const
If matrix indices are in the global range, this function returns true. Otherwise, this function retur...
Definition Ifpack2_SparsityFilter_def.hpp:263
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74