IFPACK Development
Loading...
Searching...
No Matches
Ifpack_LocalFilter.cpp
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
45#include "Epetra_MultiVector.h"
46#include "Epetra_Vector.h"
47#include "Epetra_RowMatrix.h"
48#include "Epetra_Map.h"
49#include "Epetra_BlockMap.h"
50#include "Ifpack_LocalFilter.h"
51
52//==============================================================================
53Ifpack_LocalFilter::Ifpack_LocalFilter(const Teuchos::RefCountPtr<const Epetra_RowMatrix>& Matrix) :
54 Matrix_(Matrix),
55 NumRows_(0),
56 NumNonzeros_(0),
57 MaxNumEntries_(0),
58 MaxNumEntriesA_(0)
59{
60 sprintf(Label_,"%s","Ifpack_LocalFilter");
61
62#ifdef HAVE_MPI
63 SerialComm_ = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_SELF) );
64#else
65 SerialComm_ = Teuchos::rcp( new Epetra_SerialComm );
66#endif
67
68 // localized matrix has all the local rows of Matrix
69 NumRows_ = Matrix->NumMyRows();
70
71#if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
72 // build a linear map, based on the serial communicator
73 Map_ = Teuchos::rcp( new Epetra_Map(NumRows_,0,*SerialComm_) );
74#endif
75
76 // NumEntries_ will contain the actual number of nonzeros
77 // for each localized row (that is, without external nodes,
78 // and always with the diagonal entry)
79 NumEntries_.resize(NumRows_);
80
81 // want to store the diagonal vector. FIXME: am I really useful?
82 Diagonal_ = Teuchos::rcp( new Epetra_Vector(*Map_) );
83 if (Diagonal_ == Teuchos::null) IFPACK_CHK_ERRV(-5);
84
85 // store this for future access to ExtractMyRowCopy().
86 // This is the # of nonzeros in the non-local matrix
87 MaxNumEntriesA_ = Matrix->MaxNumEntries();
88 // tentative value for MaxNumEntries. This is the number of
89 // nonzeros in the local matrix
90 MaxNumEntries_ = Matrix->MaxNumEntries();
91
92 // ExtractMyRowCopy() will use these vectors
93 Indices_.resize(MaxNumEntries_);
94 Values_.resize(MaxNumEntries_);
95
96 // now compute:
97 // - the number of nonzero per row
98 // - the total number of nonzeros
99 // - the diagonal entries
100
101 // compute nonzeros (total and per-row), and store the
102 // diagonal entries (already modified)
103 int ActualMaxNumEntries = 0;
104
105 for (int i = 0 ; i < NumRows_ ; ++i) {
106
107 NumEntries_[i] = 0;
108 int Nnz, NewNnz = 0;
109 IFPACK_CHK_ERRV(ExtractMyRowCopy(i,MaxNumEntries_,Nnz,&Values_[0],&Indices_[0]));
110
111 for (int j = 0 ; j < Nnz ; ++j) {
112 if (Indices_[j] < NumRows_ ) ++NewNnz;
113
114 if (Indices_[j] == i)
115 (*Diagonal_)[i] = Values_[j];
116 }
117
118 if (NewNnz > ActualMaxNumEntries)
119 ActualMaxNumEntries = NewNnz;
120
121 NumNonzeros_ += NewNnz;
122 NumEntries_[i] = NewNnz;
123
124 }
125
126 MaxNumEntries_ = ActualMaxNumEntries;
127}
128
129//==============================================================================
131ExtractMyRowCopy(int MyRow, int Length, int & NumEntries,
132 double *Values, int * Indices) const
133{
134 if ((MyRow < 0) || (MyRow >= NumRows_)) {
135 IFPACK_CHK_ERR(-1); // range not valid
136 }
137
138 if (Length < NumEntries_[MyRow])
139 return(-1);
140
141 // always extract using the object Values_ and Indices_.
142 // This is because I need more space than that given by
143 // the user (for the external nodes)
144 int Nnz;
145 int ierr = Matrix_->ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz,
146 &Values_[0],&Indices_[0]);
147
148 IFPACK_CHK_ERR(ierr);
149
150 // populate the user's vectors, add diagonal if not found
151 NumEntries = 0;
152
153 for (int j = 0 ; j < Nnz ; ++j) {
154 // only local indices
155 if (Indices_[j] < NumRows_ ) {
156 Indices[NumEntries] = Indices_[j];
157 Values[NumEntries] = Values_[j];
158 ++NumEntries;
159 }
160 }
161
162 return(0);
163
164}
165
166//==============================================================================
168{
169 if (!Diagonal.Map().SameAs(*Map_))
170 IFPACK_CHK_ERR(-1);
171 Diagonal = *Diagonal_;
172 return(0);
173}
174
175//==============================================================================
176int Ifpack_LocalFilter::Apply(const Epetra_MultiVector& X,
177 Epetra_MultiVector& Y) const
178{
179
180 // skip expensive checks, I suppose input data are ok
181
182 Y.PutScalar(0.0);
183 int NumVectors = Y.NumVectors();
184
185 double** X_ptr;
186 double** Y_ptr;
187 X.ExtractView(&X_ptr);
188 Y.ExtractView(&Y_ptr);
189
190 for (int i = 0 ; i < NumRows_ ; ++i) {
191
192 int Nnz;
193 int ierr = Matrix_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,&Values_[0],
194 &Indices_[0]);
195 IFPACK_CHK_ERR(ierr);
196
197 for (int j = 0 ; j < Nnz ; ++j) {
198 if (Indices_[j] < NumRows_ ) {
199 for (int k = 0 ; k < NumVectors ; ++k)
200 Y_ptr[k][i] += Values_[j] * X_ptr[k][Indices_[j]];
201 }
202 }
203 }
204
205 return(0);
206}
207
208//==============================================================================
209int Ifpack_LocalFilter::ApplyInverse(const Epetra_MultiVector& /* X */,
210 Epetra_MultiVector& /* Y */) const
211{
212 IFPACK_CHK_ERR(-1); // not implemented
213}
214
215//==============================================================================
216const Epetra_BlockMap& Ifpack_LocalFilter::Map() const
217{
218 return(*Map_);
219}
bool SameAs(const Epetra_BlockMap &Map) const
const Epetra_BlockMap & Map() const
int NumVectors() const
int ExtractView(double **A, int *MyLDA) const
int PutScalar(double ScalarConstant)
Ifpack_LocalFilter(const Teuchos::RefCountPtr< const Epetra_RowMatrix > &Matrix)
Constructor.
virtual int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.