Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_SingletonFilter.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#include "Epetra_Import.h"
52
53//==============================================================================
54Ifpack_SingletonFilter::Ifpack_SingletonFilter(const Teuchos::RefCountPtr<Epetra_RowMatrix>& Matrix) :
55 A_(Matrix),
56 NumSingletons_(0),
57 NumRows_(0),
58 NumRowsA_(0),
59 MaxNumEntries_(0),
60 MaxNumEntriesA_(0),
61 NumNonzeros_(0)
62{
63 using std::cerr;
64 using std::endl;
65
66 // use this filter only on serial matrices
67 if (A_->Comm().NumProc() != 1) {
68 cerr << "Ifpack_SingletonFilter can be used with Comm().NumProc() == 1" << endl;
69 cerr << "only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
70 cerr << "and it is not meant to be used otherwise." << endl;
71 exit(EXIT_FAILURE);
72 }
73
74 if ((A_->NumMyRows() != A_->NumGlobalRows64()) ||
75 (A_->NumMyRows() != A_->NumMyCols()))
77
78 NumRowsA_ = (A_->NumMyRows());
79 MaxNumEntriesA_ = A_->MaxNumEntries();
80
83 Reorder_.resize(A_->NumMyRows());
84
85 for (int i = 0 ; i < NumRowsA_ ; ++i)
86 Reorder_[i] = -1;
87
88 // first check how may singletons I do have
89 for (int i = 0 ; i < NumRowsA_ ; ++i) {
90 int Nnz;
91 IFPACK_CHK_ERRV(A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
92 &Values_[0], &Indices_[0]));
93 if (Nnz != 1) {
94 Reorder_[i] = NumRows_++;
95 }
96 else {
98 }
99 }
100
101 InvReorder_.resize(NumRows_);
102 for (int i = 0 ; i < NumRowsA_ ; ++i) {
103 if (Reorder_[i] < 0)
104 continue;
105 InvReorder_[Reorder_[i]] = i;
106 }
107
108 NumEntries_.resize(NumRows_);
110
111 // now compute the nonzeros per row
112 int count = 0;
113 for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
114
115 int Nnz;
116 IFPACK_CHK_ERRV(A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
117 &Values_[0], &Indices_[0]));
118 int ii = Reorder_[i];
119 if (ii >= 0) {
120 assert (Nnz != 1);
121
122 NumEntries_[ii] = Nnz;
123 NumNonzeros_ += Nnz;
124 if (Nnz > MaxNumEntries_)
125 MaxNumEntries_ = Nnz;
126 }
127 else {
128 SingletonIndex_[count] = i;
129 count++;
130 }
131 }
132
133#if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
134 Map_ = Teuchos::rcp( new Epetra_Map(NumRows_,0,Comm()) );
135#endif
136
137 // and finish up with the diagonal entry
138 Diagonal_ = Teuchos::rcp( new Epetra_Vector(*Map_) );
139
140 Epetra_Vector Diagonal(A_->Map());
141 A_->ExtractDiagonalCopy(Diagonal);
142 for (int i = 0 ; i < NumRows_ ; ++i) {
143 int ii = InvReorder_[i];
144 (*Diagonal_)[i] = Diagonal[ii];
145 }
146
147}
148
149//==============================================================================
151ExtractMyRowCopy(int MyRow, int Length, int & NumEntries,
152 double *Values, int * Indices) const
153{
154 int Nnz;
155
156 if (Length < NumEntries_[MyRow])
157 IFPACK_CHK_ERR(-1);
158
159 int Row = InvReorder_[MyRow];
160 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(Row,MaxNumEntriesA_,Nnz,
161 &Values_[0],&Indices_[0]));
162 NumEntries = 0;
163 for (int i = 0 ; i < Nnz ; ++i) {
164 int ii = Reorder_[Indices_[i]];
165 if ( ii >= 0) {
166 Indices[NumEntries] = ii;
167 Values[NumEntries] = Values_[i];
168 NumEntries++;
169 }
170 }
171 return(0);
172}
173
174//==============================================================================
176ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
177{
178 Diagonal = *Diagonal_;
179 return(0);
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 // cycle over all rows of the original matrix
198 for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
199
200 if (Reorder_[i] < 0)
201 continue; // skip singleton rows
202
203 int Nnz;
204 A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
205 &Values[0], &Indices[0]);
206 if (!TransA) {
207 // no transpose first
208 for (int j = 0 ; j < NumVectors ; ++j) {
209 for (int k = 0 ; k < Nnz ; ++k) {
210 if (Reorder_[i] >= 0)
211 Y[j][i] += Values[k] * X[j][Reorder_[Indices[k]]];
212 }
213 }
214 }
215 else {
216 // transpose here
217 for (int j = 0 ; j < NumVectors ; ++j) {
218 for (int k = 0 ; k < Nnz ; ++k) {
219 if (Reorder_[i] >= 0)
220 Y[j][Reorder_[Indices[k]]] += Values[k] * X[j][i];
221 }
222 }
223 }
224 }
225
226 return(0);
227}
228
229//==============================================================================
231Solve(bool /* Upper */, bool /* Trans */, bool /* UnitDiagonal */,
232 const Epetra_MultiVector& /* X */, Epetra_MultiVector& /* Y */) const
233{
234 return(-1);
235}
236
237//==============================================================================
240{
241 IFPACK_CHK_ERR(Multiply(false,X,Y));
242 return(0);
243}
244
245//==============================================================================
247ApplyInverse(const Epetra_MultiVector& /* X */, Epetra_MultiVector& /* Y */) const
248{
249 return(-1); // NOT IMPLEMENTED AT THIS STAGE
250}
251
252//==============================================================================
256{
257 for (int i = 0 ; i < NumSingletons_ ; ++i) {
258 int ii = SingletonIndex_[i];
259 // get the diagonal value for the singleton
260 int Nnz;
261 A_->ExtractMyRowCopy(ii,MaxNumEntriesA_,Nnz,
262 &Values_[0], &Indices_[0]);
263 for (int j = 0 ; j < Nnz ; ++j) {
264 if (Indices_[j] == ii) {
265 for (int k = 0 ; k < LHS.NumVectors() ; ++k)
266 LHS[k][ii] = RHS[k][ii] / Values_[j];
267 }
268 }
269 }
270
271 return(0);
272}
273
274//==============================================================================
277 const Epetra_MultiVector& RHS,
278 Epetra_MultiVector& ReducedRHS)
279{
280 int NumVectors = LHS.NumVectors();
281
282 for (int i = 0 ; i < NumRows_ ; ++i)
283 for (int k = 0 ; k < NumVectors ; ++k)
284 ReducedRHS[k][i] = RHS[k][InvReorder_[i]];
285
286 for (int i = 0 ; i < NumRows_ ; ++i) {
287 int ii = InvReorder_[i];
288 int Nnz;
289 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(ii,MaxNumEntriesA_,Nnz,
290 &Values_[0], &Indices_[0]));
291
292 for (int j = 0 ; j < Nnz ; ++j) {
293 if (Reorder_[Indices_[j]] == -1) {
294 for (int k = 0 ; k < NumVectors ; ++k)
295 ReducedRHS[k][i] -= Values_[j] * LHS[k][Indices_[j]];
296 }
297 }
298 }
299 return(0);
300}
301
302//==============================================================================
304UpdateLHS(const Epetra_MultiVector& ReducedLHS,
306{
307 for (int i = 0 ; i < NumRows_ ; ++i)
308 for (int k = 0 ; k < LHS.NumVectors() ; ++k)
309 LHS[k][InvReorder_[i]] = ReducedLHS[k][i];
310
311 return(0);
312}
#define IFPACK_CHK_ERR(ifpack_err)
#define IFPACK_CHK_ERRV(ifpack_err)
#define RHS(a)
Definition MatGenFD.c:60
int NumVectors() const
std::vector< int > Indices_
Used in ExtractMyRowCopy, to avoid allocation each time.
Ifpack_SingletonFilter(const Teuchos::RefCountPtr< Epetra_RowMatrix > &Matrix)
Constructor.
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Teuchos::RefCountPtr< Epetra_Map > Map_
std::vector< int > SingletonIndex_
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int SolveSingletons(const Epetra_MultiVector &RHS, Epetra_MultiVector &LHS)
int UpdateLHS(const Epetra_MultiVector &ReducedLHS, Epetra_MultiVector &LHS)
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Pointer to the matrix to be preconditioned.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
int CreateReducedRHS(const Epetra_MultiVector &LHS, const Epetra_MultiVector &RHS, Epetra_MultiVector &ReducedRHS)
Teuchos::RefCountPtr< Epetra_Vector > Diagonal_
virtual int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
std::vector< double > Values_
Used in ExtractMyRowCopy, to avoid allocation each time.
const Epetra_Comm & Comm() const
const int NumVectors