Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_DropFilter_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_DROPFILTER_DEF_HPP
44#define IFPACK2_DROPFILTER_DEF_HPP
45#include "Ifpack2_DropFilter_decl.hpp"
46#include <vector>
47
48#include "Tpetra_ConfigDefs.hpp"
49#include "Tpetra_RowMatrix.hpp"
50#include "Tpetra_Map.hpp"
51#include "Tpetra_MultiVector.hpp"
52#include "Tpetra_Vector.hpp"
53
54namespace Ifpack2 {
55
56//==========================================================================
57template<class MatrixType>
58DropFilter<MatrixType>::DropFilter(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Matrix,
59 magnitudeType DropTol):
60 A_(Matrix),
61 DropTol_(DropTol),
62 NumRows_(0),
63 NumNonzeros_(0),
64 MaxNumEntries_(0),
65 MaxNumEntriesA_(0)
66{
67
68 // use this filter only on serial matrices
69 if (A_->getComm()->getSize() != 1 || A_->getLocalNumRows() != A_->getGlobalNumRows()) {
70 throw std::runtime_error("Ifpack2::DropFilter 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.");
71 }
72
73
74 // localized matrix has all the local rows of Matrix
75 NumRows_ = A_->getLocalNumRows();
76
77 // NodeNumEntries_ will contain the actual number of nonzeros
78 // for each localized row (that is, without external nodes,
79 // and always with the diagonal entry)
80 NumEntries_.resize(NumRows_);
81
82 // tentative value for MaxNumEntries. This is the number of
83 // nonzeros in the local matrix
84 MaxNumEntries_ = A_->getLocalMaxNumRowEntries();
85 MaxNumEntriesA_ = A_->getLocalMaxNumRowEntries();
86
87 // ExtractMyRowCopy() will use these vectors
88 Kokkos::resize(Indices_,MaxNumEntries_);
89 Kokkos::resize(Values_,MaxNumEntries_);
90
91 size_t ActualMaxNumEntries = 0;
92 for (size_t i = 0 ; i < NumRows_ ; ++i) {
93 NumEntries_[i] = MaxNumEntriesA_;
94 size_t Nnz, NewNnz=0;
95 A_->getLocalRowCopy(i,Indices_,Values_,Nnz);
96 for (size_t j = 0 ; j < Nnz ; ++j) {
97 if (((size_t)Indices_[j] == i) || (Teuchos::ScalarTraits<Scalar>::magnitude(Values_[j]) >= DropTol_))
98 NewNnz++;
99 }
100
101 NumNonzeros_ += NewNnz;
102 NumEntries_[i] = NewNnz;
103 if (NewNnz > ActualMaxNumEntries)
104 ActualMaxNumEntries = NewNnz;
105 }
106
107 MaxNumEntries_ = ActualMaxNumEntries;
108}
109
110//=========================================================================
111template<class MatrixType>
113
114//==========================================================================
115template<class MatrixType>
116Teuchos::RCP<const Teuchos::Comm<int> >
118{
119 return A_->getComm ();
120}
121
122
123//==========================================================================
124template<class MatrixType>
125Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
126 typename MatrixType::global_ordinal_type,
127 typename MatrixType::node_type> >
129{
130 return A_->getRowMap();
131}
132
133//==========================================================================
134template<class MatrixType>
135Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
136 typename MatrixType::global_ordinal_type,
137 typename MatrixType::node_type> >
139{
140 return A_->getColMap();
141}
142
143//==========================================================================
144template<class MatrixType>
145Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
146 typename MatrixType::global_ordinal_type,
147 typename MatrixType::node_type> >
149{
150 return A_->getDomainMap();
151}
152
153//==========================================================================
154template<class MatrixType>
155Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
156 typename MatrixType::global_ordinal_type,
157 typename MatrixType::node_type> >
159{
160 return A_->getRangeMap();
161}
162
163//==========================================================================
164template<class MatrixType>
165Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type,
166 typename MatrixType::global_ordinal_type,
167 typename MatrixType::node_type> >
169{
170 throw std::runtime_error("Ifpack2::DropFilter: does not support getGraph.");
171}
172
173//==========================================================================
174template<class MatrixType>
176{
177 return NumRows_;
178}
179
180//==========================================================================
181template<class MatrixType>
183{
184 return NumRows_;
185}
186
187//==========================================================================
188template<class MatrixType>
190{
191 return NumRows_;
192}
193
194//==========================================================================
195
196template<class MatrixType>
198{
199 return NumRows_;
200}
201
202//==========================================================================
203template<class MatrixType>
204typename MatrixType::global_ordinal_type DropFilter<MatrixType>::getIndexBase() const
205{
206 return A_->getIndexBase();
207}
208
209//==========================================================================
210template<class MatrixType>
212{
213 return NumNonzeros_;
214}
215
216//==========================================================================
217template<class MatrixType>
219{
220 return NumNonzeros_;
221}
222
223//==========================================================================
224template<class MatrixType>
225size_t DropFilter<MatrixType>::getNumEntriesInGlobalRow(GlobalOrdinal /* globalRow */) const
226{
227 throw std::runtime_error("Ifpack2::DropFilter does not implement getNumEntriesInGlobalRow.");
228}
229
230//==========================================================================
231template<class MatrixType>
232size_t DropFilter<MatrixType>::getNumEntriesInLocalRow(LocalOrdinal localRow) const
233{
234 return NumEntries_[localRow];
235}
236
237//==========================================================================
238template<class MatrixType>
240{
241 return MaxNumEntries_;
242}
243
244//==========================================================================
245template<class MatrixType>
247{
248 return MaxNumEntries_;
249}
250
251//==========================================================================
252template<class MatrixType>
254{
255 return true;
256}
257
258//==========================================================================
259template<class MatrixType>
261{
262 return A_->isLocallyIndexed();
263}
264
265//==========================================================================
266template<class MatrixType>
268{
269 return A_->isGloballyIndexed();
270}
271
272//==========================================================================
273template<class MatrixType>
275{
276 return A_->isFillComplete();
277}
278
279//==========================================================================
280template<class MatrixType>
282getGlobalRowCopy (GlobalOrdinal /*GlobalRow*/,
283 nonconst_global_inds_host_view_type &/*Indices*/,
284 nonconst_values_host_view_type &/*Values*/,
285 size_t& /*NumEntries*/) const
286{
287 throw std::runtime_error("Ifpack2::DropFilter does not implement getGlobalRowCopy.");
288}
289
290//==========================================================================
291template<class MatrixType>
293 getLocalRowCopy (LocalOrdinal LocalRow,
294 nonconst_local_inds_host_view_type &Indices,
295 nonconst_values_host_view_type &Values,
296 size_t& NumEntries) const
297{
298 TEUCHOS_TEST_FOR_EXCEPTION((LocalRow < 0 || (size_t) LocalRow >= NumRows_ || (size_t) Indices.size() < NumEntries_[LocalRow]), std::runtime_error, "Ifpack2::DropFilter::getLocalRowCopy invalid row or array size.");
299
300 // Note: This function will work correctly if called by apply, say, with Indices_ and Values_ as
301 // parameters. The structure of the loop below should make that obvious.
302
303 // always extract using the object Values_ and Indices_.
304 // This is because I need more space than that given by
305 // the user (for the external nodes)
306 size_t A_NumEntries=0;
307 A_->getLocalRowCopy(LocalRow,Indices_,Values_,A_NumEntries);
308
309 // loop over all nonzero elements of row MyRow,
310 // and drop elements below specified threshold.
311 // Also, add value to zero diagonal
312 NumEntries = 0;
313 for (size_t i = 0 ; i < A_NumEntries ; ++i) {
314 // if element is above specified tol, add to the
315 // user's defined arrays. Check that we are not
316 // exceeding allocated space. Do not drop any diagonal entry.
317 if ((Indices_[i] == LocalRow) || (Teuchos::ScalarTraits<Scalar>::magnitude(Values_[i]) >= DropTol_)) {
318 Values[NumEntries] = Values_[i];
319 Indices[NumEntries] = Indices_[i];
320 NumEntries++;
321 }
322 }
323}
324
325//==========================================================================
326template<class MatrixType>
327void DropFilter<MatrixType>::getGlobalRowView(GlobalOrdinal /* GlobalRow */,
328 global_inds_host_view_type &/*indices*/,
329 values_host_view_type &/*values*/) const
330{
331 throw std::runtime_error("Ifpack2::DropFilter: does not support getGlobalRowView.");
332}
333
334//==========================================================================
335template<class MatrixType>
336void DropFilter<MatrixType>::getLocalRowView(LocalOrdinal /* LocalRow */,
337 local_inds_host_view_type & /*indices*/,
338 values_host_view_type & /*values*/) const
339{
340 throw std::runtime_error("Ifpack2::DropFilter: does not support getLocalRowView.");
341}
342
343//==========================================================================
344template<class MatrixType>
345void DropFilter<MatrixType>::getLocalDiagCopy(Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const
346{
347 // This is somewhat dubious as to how the maps match.
348 return A_->getLocalDiagCopy(diag);
349}
350
351//==========================================================================
352template<class MatrixType>
353void DropFilter<MatrixType>::leftScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* x */)
354{
355 throw std::runtime_error("Ifpack2::DropFilter does not support leftScale.");
356}
357
358//==========================================================================
359template<class MatrixType>
360void DropFilter<MatrixType>::rightScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* x */)
361{
362 throw std::runtime_error("Ifpack2::DropFilter does not support rightScale.");
363}
364
365//==========================================================================
366template<class MatrixType>
367void DropFilter<MatrixType>::apply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
368 Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
369 Teuchos::ETransp mode,
370 Scalar /* alpha */,
371 Scalar /* beta */) const
372{
373 // Note: This isn't AztecOO compliant. But neither was Ifpack's version.
374 // Note: The localized maps mean the matvec is trivial (and has no import)
375 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
376 "Ifpack2::DropFilter::apply ERROR: X.getNumVectors() != Y.getNumVectors().");
377
378 Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
379 Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > x_ptr = X.get2dView();
380 Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > y_ptr = Y.get2dViewNonConst();
381
382 Y.putScalar(zero);
383 size_t NumVectors = Y.getNumVectors();
384
385 for (size_t i = 0 ; i < NumRows_ ; ++i) {
386 size_t Nnz;
387 // Use this class's getrow to make the below code simpler
388 getLocalRowCopy(i,Indices_,Values_,Nnz);
389 Scalar* Values = reinterpret_cast<Scalar*>(Values_.data());
390 if (mode==Teuchos::NO_TRANS){
391 for (size_t j = 0 ; j < Nnz ; ++j)
392 for (size_t k = 0 ; k < NumVectors ; ++k)
393 y_ptr[k][i] += Values[j] * x_ptr[k][Indices_[j]];
394 }
395 else if (mode==Teuchos::TRANS){
396 for (size_t j = 0 ; j < Nnz ; ++j)
397 for (size_t k = 0 ; k < NumVectors ; ++k)
398 y_ptr[k][Indices_[j]] += Values[j] * x_ptr[k][i];
399 }
400 else { //mode==Teuchos::CONJ_TRANS
401 for (size_t j = 0 ; j < Nnz ; ++j)
402 for (size_t k = 0 ; k < NumVectors ; ++k)
403 y_ptr[k][Indices_[j]] += Teuchos::ScalarTraits<Scalar>::conjugate(Values[j]) * x_ptr[k][i];
404 }
405 }
406}
407
408
409//==========================================================================
410template<class MatrixType>
412{
413 return true;
414}
415
416//==========================================================================
417template<class MatrixType>
419{
420 return false;
421}
422
423//==========================================================================
424template<class MatrixType>
425typename DropFilter<MatrixType>::mag_type DropFilter<MatrixType>::getFrobeniusNorm() const
426{
427 throw std::runtime_error("Ifpack2::DropFilter does not implement getFrobeniusNorm.");
428}
429
430
431#define IFPACK2_DROPFILTER_INSTANT(S,LO,GO,N) \
432 template class Ifpack2::DropFilter< Tpetra::RowMatrix<S, LO, GO, N> >;
433
434
435} // namespace Ifpack2
436
437#endif
virtual size_t getLocalNumRows() const
Returns the number of rows owned on the calling node.
Definition Ifpack2_DropFilter_def.hpp:189
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getColMap() const
Returns the Map that describes the column distribution in this matrix.
Definition Ifpack2_DropFilter_def.hpp:138
virtual void rightScale(const Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &x)
Scales the RowMatrix on the right with the Vector x.
Definition Ifpack2_DropFilter_def.hpp:360
virtual bool hasColMap() const
Indicates whether this matrix has a well-defined column map.
Definition Ifpack2_DropFilter_def.hpp:253
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_DropFilter_def.hpp:336
virtual global_size_t getGlobalNumRows() const
Returns the number of global rows in this matrix.
Definition Ifpack2_DropFilter_def.hpp:175
virtual bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
Definition Ifpack2_DropFilter_def.hpp:411
virtual bool supportsRowViews() const
Returns true if RowViews are supported.
Definition Ifpack2_DropFilter_def.hpp:418
virtual void leftScale(const Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &x)
Scales the RowMatrix on the left with the Vector x.
Definition Ifpack2_DropFilter_def.hpp:353
virtual GlobalOrdinal getIndexBase() const
Returns the index base for global indices for this matrix.
Definition Ifpack2_DropFilter_def.hpp:204
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
Definition Ifpack2_DropFilter_def.hpp:232
virtual Teuchos::RCP< const Tpetra::RowGraph< LocalOrdinal, GlobalOrdinal, Node > > getGraph() const
Returns the RowGraph associated with this matrix.
Definition Ifpack2_DropFilter_def.hpp:168
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition Ifpack2_DropFilter_def.hpp:128
virtual bool isGloballyIndexed() const
If matrix indices are in the global range, this function returns true. Otherwise, this function retur...
Definition Ifpack2_DropFilter_def.hpp:267
virtual size_t getLocalNumEntries() const
Returns the local number of entries in this matrix.
Definition Ifpack2_DropFilter_def.hpp:218
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_DropFilter_def.hpp:282
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
Returns the Map that describes the range distribution in this matrix.
Definition Ifpack2_DropFilter_def.hpp:158
DropFilter(const Teuchos::RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &Matrix, magnitudeType DropTol)
Constructor.
Definition Ifpack2_DropFilter_def.hpp:58
virtual bool isLocallyIndexed() const
If matrix indices are in the local range, this function returns true. Otherwise, this function return...
Definition Ifpack2_DropFilter_def.hpp:260
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
Returns the Map that describes the domain distribution in this matrix.
Definition Ifpack2_DropFilter_def.hpp:148
virtual mag_type getFrobeniusNorm() const
Returns the Frobenius norm of the matrix.
Definition Ifpack2_DropFilter_def.hpp:425
virtual bool isFillComplete() const
Returns true if fillComplete() has been called.
Definition Ifpack2_DropFilter_def.hpp:274
virtual size_t getLocalNumCols() const
Returns the number of columns needed to apply the forward operator on this node, i....
Definition Ifpack2_DropFilter_def.hpp:197
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_DropFilter_def.hpp:327
virtual size_t getLocalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on this node.
Definition Ifpack2_DropFilter_def.hpp:246
virtual void getLocalRowCopy(LocalOrdinal DropRow, 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_DropFilter_def.hpp:293
virtual ~DropFilter()
Destructor.
Definition Ifpack2_DropFilter_def.hpp:112
virtual size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
Returns the current number of entries on this node in the specified global row.
Definition Ifpack2_DropFilter_def.hpp:225
virtual global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
Definition Ifpack2_DropFilter_def.hpp:211
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_DropFilter_def.hpp:367
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_DropFilter_def.hpp:345
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the communicator.
Definition Ifpack2_DropFilter_def.hpp:117
virtual global_size_t getGlobalNumCols() const
Returns the number of global columns in this matrix.
Definition Ifpack2_DropFilter_def.hpp:182
virtual size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
Definition Ifpack2_DropFilter_def.hpp:239
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74