Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_DropFilter.cpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) 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#include "Ifpack_ConfigDefs.h"
44#include "Ifpack_DropFilter.h"
45#include "Epetra_ConfigDefs.h"
46#include "Epetra_RowMatrix.h"
47#include "Epetra_Comm.h"
48#include "Epetra_Map.h"
49#include "Epetra_MultiVector.h"
50#include "Epetra_Vector.h"
51
52//==============================================================================
53Ifpack_DropFilter::Ifpack_DropFilter(const Teuchos::RefCountPtr<Epetra_RowMatrix>& Matrix,
54 double DropTol) :
55 A_(Matrix),
56 DropTol_(DropTol),
57 MaxNumEntries_(0),
58 MaxNumEntriesA_(0),
59 NumNonzeros_(0)
60{
61 using std::cerr;
62 using std::endl;
63
64 // use this filter only on serial matrices
65 if (A_->Comm().NumProc() != 1) {
66 cerr << "Ifpack_DropFilter can be used with Comm().NumProc() == 1" << endl;
67 cerr << "only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
68 cerr << "and it is not meant to be used otherwise." << endl;
69 exit(EXIT_FAILURE);
70 }
71
72 if ((A_->NumMyRows() != A_->NumGlobalRows64()) ||
73 (A_->NumMyRows() != A_->NumMyCols()))
75
76 NumRows_ = A_->NumMyRows();
77 MaxNumEntriesA_ = A_->MaxNumEntries();
78
79 NumEntries_.resize(NumRows_);
82
83 std::vector<int> Ind(MaxNumEntriesA_);
84 std::vector<double> Val(MaxNumEntriesA_);
85
86 for (int i = 0 ; i < NumRows_ ; ++i) {
88 int Nnz;
90 &Val[0], &Ind[0]));
91
92 NumEntries_[i] = Nnz;
93 NumNonzeros_ += Nnz;
94 if (Nnz > MaxNumEntries_)
95 MaxNumEntries_ = Nnz;
96 }
97
98}
99
100//==============================================================================
102ExtractMyRowCopy(int MyRow, int Length, int & NumEntries,
103 double *Values, int * Indices) const
104{
105 if (Length < NumEntries_[MyRow])
106 IFPACK_CHK_ERR(-1);
107
108 int Nnz;
109
110 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz,
111 &Values_[0],&Indices_[0]));
112
113 // loop over all nonzero elements of row MyRow,
114 // and drop elements below specified threshold.
115 // Also, add value to zero diagonal
116 int count = 0;
117 for (int i = 0 ; i < Nnz ; ++i) {
118
119 // if element is above specified tol, add to the
120 // user's defined arrays. Check that we are not
121 // exceeding allocated space. Do not drop any diagonal entry.
122 if ((Indices_[i] == MyRow) || (IFPACK_ABS(Values_[i]) >= DropTol_)) {
123 if (count == Length)
124 IFPACK_CHK_ERR(-1);
125 Values[count] = Values_[i];
126 Indices[count] = Indices_[i];
127 count++;
128 }
129 }
130
131 NumEntries = count;
132 return(0);
133}
134
135//==============================================================================
137ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
138{
139 IFPACK_CHK_ERR(A_->ExtractDiagonalCopy(Diagonal));
140 return(0);
141}
142
143//==============================================================================
145Multiply(bool TransA, const Epetra_MultiVector& X,
146 Epetra_MultiVector& Y) const
147{
148 // NOTE: I suppose that the matrix has been localized,
149 // hence all maps are trivial.
150 int NumVectors = X.NumVectors();
151 if (NumVectors != Y.NumVectors())
152 IFPACK_CHK_ERR(-1);
153
154 Y.PutScalar(0.0);
155
156 std::vector<int> Indices(MaxNumEntries_);
157 std::vector<double> Values(MaxNumEntries_);
158
159 for (int i = 0 ; i < NumRows_ ; ++i) {
160
161 int Nnz;
163 &Values[0], &Indices[0]);
164 if (!TransA) {
165 // no transpose first
166 for (int j = 0 ; j < NumVectors ; ++j) {
167 for (int k = 0 ; k < Nnz ; ++k) {
168 Y[j][i] += Values[k] * X[j][Indices[k]];
169 }
170 }
171 }
172 else {
173 // transpose here
174 for (int j = 0 ; j < NumVectors ; ++j) {
175 for (int k = 0 ; k < Nnz ; ++k) {
176 Y[j][Indices[k]] += Values[k] * X[j][i];
177 }
178 }
179 }
180 }
181
182 return(0);
183}
184
185//==============================================================================
187Solve(bool /* Upper */, bool /* Trans */, bool /* UnitDiagonal */,
188 const Epetra_MultiVector& /* X */, Epetra_MultiVector& /* Y */) const
189{
190 IFPACK_CHK_ERR(-99);
191}
192
193//==============================================================================
196{
197 int ierr = Multiply(UseTranspose(),X,Y);
198 IFPACK_RETURN(ierr);
199}
200
201//==============================================================================
203ApplyInverse(const Epetra_MultiVector& /* X */, Epetra_MultiVector& /* Y */) const
204{
205 IFPACK_CHK_ERR(-99);
206}
207
208//==============================================================================
210{
211 IFPACK_CHK_ERR(-1);
212}
213
214//==============================================================================
216{
217 IFPACK_CHK_ERR(-1);
218}
#define IFPACK_RETURN(ifpack_err)
#define IFPACK_CHK_ERR(ifpack_err)
#define IFPACK_ABS(x)
#define IFPACK_CHK_ERRV(ifpack_err)
int NumVectors() const
int PutScalar(double ScalarConstant)
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int InvRowSums(Epetra_Vector &x) const
int NumNonzeros_
Number of nonzeros for the dropped matrix.
virtual int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
std::vector< int > Indices_
Used in ExtractMyRowCopy, to avoid allocation each time.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
double DropTol_
Drop tolerance.
int MaxNumEntries_
Maximum entries in each row.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int InvColSums(Epetra_Vector &x) const
std::vector< double > Values_
Used in ExtractMyRowCopy, to avoid allocation each time.
virtual int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Pointer to the matrix to be preconditioned.
Ifpack_DropFilter(const Teuchos::RefCountPtr< Epetra_RowMatrix > &Matrix, double DropTol)
Constructor.
std::vector< int > NumEntries_
bool UseTranspose() const
const int NumVectors
Definition: performance.cpp:71