Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 //==============================================================================
53 Ifpack_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()))
74  IFPACK_CHK_ERRV(-2);
75 
76  NumRows_ = A_->NumMyRows();
77  MaxNumEntriesA_ = A_->MaxNumEntries();
78 
79  NumEntries_.resize(NumRows_);
80  Indices_.resize(MaxNumEntriesA_);
81  Values_.resize(MaxNumEntriesA_);
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 //==============================================================================
102 ExtractMyRowCopy(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 //==============================================================================
138 {
139  IFPACK_CHK_ERR(A_->ExtractDiagonalCopy(Diagonal));
140  return(0);
141 }
142 
143 //==============================================================================
145 Multiply(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 //==============================================================================
187 Solve(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 //==============================================================================
203 ApplyInverse(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 }
virtual int InvColSums(Epetra_Vector &x) const
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
const int NumVectors
Definition: performance.cpp:71
double DropTol_
Drop tolerance.
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Pointer to the matrix to be preconditioned.
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
#define IFPACK_CHK_ERRV(ifpack_err)
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int InvRowSums(Epetra_Vector &x) const
int NumNonzeros_
Number of nonzeros for the dropped matrix.
int MaxNumEntries_
Maximum entries in each row.
#define IFPACK_RETURN(ifpack_err)
std::vector< double > Values_
Used in ExtractMyRowCopy, to avoid allocation each time.
#define IFPACK_ABS(x)
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
#define IFPACK_CHK_ERR(ifpack_err)
std::vector< int > NumEntries_
Ifpack_DropFilter(const Teuchos::RefCountPtr< Epetra_RowMatrix > &Matrix, double DropTol)
Constructor.
std::vector< int > Indices_
Used in ExtractMyRowCopy, to avoid allocation each time.
bool UseTranspose() const