Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
44 #include "Ifpack_SparsityFilter.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 //==============================================================================
53 Ifpack_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()))
78  IFPACK_CHK_ERRV(-1);
79 
80  NumRows_ = A_->NumMyRows();
81  MaxNumEntriesA_ = A_->MaxNumEntries();
82  Indices_.resize(MaxNumEntriesA_);
83  Values_.resize(MaxNumEntriesA_);
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 //==============================================================================
114 ExtractMyRowCopy(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 //==============================================================================
177 {
178  int ierr = A_->ExtractDiagonalCopy(Diagonal);
179  IFPACK_RETURN(ierr);
180 }
181 
182 //==============================================================================
184 Multiply(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 //==============================================================================
225 Solve(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 //==============================================================================
241 ApplyInverse(const Epetra_MultiVector& /* X */, Epetra_MultiVector& /* Y */) const
242 {
243  IFPACK_CHK_ERR(-98);
244 }
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int AllowedEntries_
Maximum allowed entries per row.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
std::vector< int > Indices_
Used in ExtractMyRowCopy, to avoid allocation each time.
virtual int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
int MaxNumEntries_
Maximum entries in each row.
const int NumVectors
Definition: performance.cpp:71
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
std::vector< double > Values_
Used in ExtractMyRowCopy, to avoid allocation each time.
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
#define IFPACK_CHK_ERRV(ifpack_err)
Ifpack_SparsityFilter(const Teuchos::RefCountPtr< Epetra_RowMatrix > &Matrix, int AllowedNumEntries, int AllowedBandwidth=-1)
#define IFPACK_RETURN(ifpack_err)
#define IFPACK_ABS(x)
std::vector< int > NumEntries_
int NumNonzeros_
Number of nonzeros for the dropped matrix.
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Pointer to the matrix to be preconditioned.
#define IFPACK_CHK_ERR(ifpack_err)
int AllowedBandwidth_
Maximum allowed bandwidth.
virtual int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const