Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_SparsityFilter.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"
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_SparsityFilter::Ifpack_SparsityFilter(const Teuchos::RefCountPtr<Epetra_RowMatrix>& Matrix,
54 int AllowedEntries,
55 int AllowedBandwidth) :
56 A_(Matrix),
57 MaxNumEntries_(0),
58 MaxNumEntriesA_(0),
59 AllowedBandwidth_(AllowedBandwidth),
60 AllowedEntries_(AllowedEntries),
61 NumNonzeros_(0),
62 NumRows_(0)
63{
64 using std::cerr;
65 using std::endl;
66
67 // use this filter only on serial matrices
68 if (A_->Comm().NumProc() != 1) {
69 cerr << "Ifpack_SparsityFilter can be used with Comm().NumProc() == 1" << endl;
70 cerr << "only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
71 cerr << "and it is not meant to be used otherwise." << endl;
72 exit(EXIT_FAILURE);
73 }
74
75 // only square serial matrices
76 if ((A_->NumMyRows() != A_->NumMyCols()) ||
77 (A_->NumMyRows() != A_->NumGlobalRows64()))
79
80 NumRows_ = A_->NumMyRows();
81 MaxNumEntriesA_ = A_->MaxNumEntries();
84
85 // default value is to not consider bandwidth
86 if (AllowedBandwidth_ == -1)
88
89 // computes the number of nonzero elements per row in the
90 // dropped matrix. Stores this number in NumEntries_.
91 // Also, computes the global number of nonzeros.
92 std::vector<int> Ind(MaxNumEntriesA_);
93 std::vector<double> Val(MaxNumEntriesA_);
94
95 NumEntries_.resize(NumRows_);
96 for (int i = 0 ; i < NumRows_ ; ++i)
98
99 for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
100 int Nnz;
102 &Val[0], &Ind[0]));
103
104 NumEntries_[i] = Nnz;
105 NumNonzeros_ += Nnz;
106 if (Nnz > MaxNumEntries_)
107 MaxNumEntries_ = Nnz;
108
109 }
110}
111
112//==============================================================================
114ExtractMyRowCopy(int MyRow, int Length, int & NumEntries,
115 double *Values, int * Indices) const
116{
117 if (Length < NumEntries_[MyRow])
118 IFPACK_CHK_ERR(-1);
119
120 int Nnz;
121 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz,
122 &Values_[0],&Indices_[0]));
123
124 double Threshold = 0.0;
125
126 // this `if' is to define the cut-off value
127 if (Nnz > AllowedEntries_) {
128
129 std::vector<double> Values2(Nnz);
130 int count = 0;
131 for (int i = 0 ; i < Nnz ; ++i) {
132 // skip diagonal entry (which is always inserted)
133 if (Indices_[i] == MyRow)
134 continue;
135 // put absolute value
136 Values2[count] = IFPACK_ABS(Values_[i]);
137 count++;
138 }
139
140 for (int i = count ; i < Nnz ; ++i)
141 Values2[i] = 0.0;
142
143 // sort in descending order
144 std::sort(Values2.rbegin(),Values2.rend());
145 // get the cut-off value
146 Threshold = Values2[AllowedEntries_ - 1];
147
148 }
149
150 // loop over all nonzero elements of row MyRow,
151 // and drop elements below specified threshold.
152 // Also, add value to zero diagonal
153 NumEntries = 0;
154
155 for (int i = 0 ; i < Nnz ; ++i) {
156
157 if (IFPACK_ABS(Indices_[i] - MyRow) > AllowedBandwidth_)
158 continue;
159
160 if ((Indices_[i] != MyRow) && (IFPACK_ABS(Values_[i]) < Threshold))
161 continue;
162
163 Values[NumEntries] = Values_[i];
164 Indices[NumEntries] = Indices_[i];
165
166 NumEntries++;
167 if (NumEntries > AllowedEntries_)
168 break;
169 }
170
171 return(0);
172}
173
174//==============================================================================
176ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
177{
178 int ierr = A_->ExtractDiagonalCopy(Diagonal);
179 IFPACK_RETURN(ierr);
180}
181
182//==============================================================================
184Multiply(bool TransA, const Epetra_MultiVector& X,
185 Epetra_MultiVector& Y) const
186{
187
188 int NumVectors = X.NumVectors();
189 if (NumVectors != Y.NumVectors())
190 IFPACK_CHK_ERR(-1);
191
192 Y.PutScalar(0.0);
193
194 std::vector<int> Indices(MaxNumEntries_);
195 std::vector<double> Values(MaxNumEntries_);
196
197 for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
198
199 int Nnz;
201 &Values[0], &Indices[0]);
202 if (!TransA) {
203 // no transpose first
204 for (int j = 0 ; j < NumVectors ; ++j) {
205 for (int k = 0 ; k < Nnz ; ++k) {
206 Y[j][i] += Values[k] * X[j][Indices[k]];
207 }
208 }
209 }
210 else {
211 // transpose here
212 for (int j = 0 ; j < NumVectors ; ++j) {
213 for (int k = 0 ; k < Nnz ; ++k) {
214 Y[j][Indices[k]] += Values[k] * X[j][i];
215 }
216 }
217 }
218 }
219
220 return(0);
221}
222
223//==============================================================================
225Solve(bool /* Upper */, bool /* Trans */, bool /* UnitDiagonal */,
226 const Epetra_MultiVector& /* X */, Epetra_MultiVector& /* Y */) const
227{
228 IFPACK_CHK_ERR(-98);
229}
230
231//==============================================================================
234{
235 int ierr = Multiply(UseTranspose(),X,Y);
236 IFPACK_RETURN(ierr);
237}
238
239//==============================================================================
241ApplyInverse(const Epetra_MultiVector& /* X */, Epetra_MultiVector& /* Y */) const
242{
243 IFPACK_CHK_ERR(-98);
244}
#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)
int MaxNumEntries_
Maximum entries in each row.
Ifpack_SparsityFilter(const Teuchos::RefCountPtr< Epetra_RowMatrix > &Matrix, int AllowedNumEntries, int AllowedBandwidth=-1)
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Pointer to the matrix to be preconditioned.
int NumNonzeros_
Number of nonzeros for the dropped matrix.
virtual int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
std::vector< double > Values_
Used in ExtractMyRowCopy, to avoid allocation each time.
std::vector< int > NumEntries_
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
virtual int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
int AllowedEntries_
Maximum allowed entries per row.
std::vector< int > Indices_
Used in ExtractMyRowCopy, to avoid allocation each time.
int AllowedBandwidth_
Maximum allowed bandwidth.
const int NumVectors
Definition: performance.cpp:71