Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_SubdomainFilter.cpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Ifpack_SubdomainFilter
5// Author: Radu Popescu <radu.popescu@epfl.ch>
6// Copyright 2011 EPFL - CMCS
7//
8// Redistribution and use in source and binary forms, with or without
9// modification, are permitted provided that the following conditions are
10// met:
11//
12// 1. Redistributions of source code must retain the above copyright
13// notice, this list of conditions and the following disclaimer.
14//
15// 2. Redistributions in binary form must reproduce the above copyright
16// notice, this list of conditions and the following disclaimer in the
17// documentation and/or other materials provided with the distribution.
18//
19// 3. Neither the name of the copyright holder nor the names of the
20// contributors may be used to endorse or promote products derived from
21// this software without specific prior written permission.
22//
23// THIS SOFTWARE IS PROVIDED BY EPFL - CMCS "AS IS" AND ANY EXPRESS OR
24// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
25// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26// ARE DISCLAIMED. IN NO EVENT SHALL EPFL - CMCS OR THE CONTRIBUTORS
27// BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
28// OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
29// OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
30// BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
31// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
32// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
33// EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34//
35// ************************************************************************
36//@HEADER
37
38#include "Ifpack_ConfigDefs.h"
39
40#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
41
42#include <vector>
43#include <algorithm>
44
45#include "Epetra_MultiVector.h"
46#include "Epetra_Vector.h"
47#include "Epetra_IntVector.h"
48#include "Epetra_RowMatrix.h"
49#include "Epetra_Map.h"
50#include "Epetra_BlockMap.h"
52#include "Epetra_Import.h"
53#include "Epetra_Export.h"
54#include "Epetra_CrsMatrix.h"
55#include "Epetra_BLAS_wrappers.h"
56
58
59using namespace Teuchos;
60
61//==============================================================================
62Ifpack_SubdomainFilter::Ifpack_SubdomainFilter(const RefCountPtr<const Epetra_RowMatrix>& Matrix,int subdomainId) :
63 Matrix_(Matrix),
64 NumMyRows_(0),
65 NumMyNonzeros_(0),
66 NumGlobalRows_(0),
67 NumGlobalCols_(0),
68 MaxNumEntries_(0),
69 MaxNumEntriesA_(0)
70{
71 sprintf(Label_,"%s","Ifpack_SubdomainFilter");
72
73 ImportVector_ = 0;
74 ExportVector_ = 0;
75
76 ovA_ = dynamic_cast<const Ifpack_OverlappingRowMatrix*>(&*Matrix_);
77
78 if (ovA_) {
79 Acrs_=dynamic_cast<const Epetra_CrsMatrix*>(&ovA_->A());
80 NumMyRowsA_ = ovA_->A().NumMyRows();
81 NumMyRowsB_ = ovA_->B().NumMyRows();
82 } else {
83 NumMyRowsA_ = Matrix->NumMyRows();
84 NumMyRowsB_ = 0;
85 }
86
87#ifdef HAVE_MPI
88 const Epetra_MpiComm *pComm = dynamic_cast<const Epetra_MpiComm*>( &(Matrix->Comm()) );
89 assert(pComm != NULL);
90 MPI_Comm_split(pComm->Comm(), subdomainId, pComm->MyPID(), &subdomainMPIComm_);
91 SubComm_ = rcp( new Epetra_MpiComm(subdomainMPIComm_) );
92#else
93 SubComm_ = rcp( new Epetra_SerialComm );
94#endif
95
96 NumMyRows_ = Matrix->NumMyRows();
97 SubComm_->SumAll(&NumMyRows_,&NumGlobalRows_,1);
98
99 // Get the row GID corresponding to the process
100 const Epetra_Map &globRowMap = Matrix->RowMatrixRowMap();
101 std::vector<int> myRowsGID(NumMyRows_);
102 globRowMap.MyGlobalElements(&myRowsGID[0]);
103
104 // Gather the GID of all rows in the subdomain
105 // Get the maximum number of local rows on each processor in the subdomain
106 // Allocate a vector large enough (Proc cout * max number local rows
107 // Gather all the local row indices. If a process has less than max num
108 // local rows, pad the local vector with -1;
109 // After the gather, eliminate the -1 elements from the target vector
110 //
111 // TODO: this section should be rewritten to avoid the GatherAll call.
112 // Right now it's not memory efficient. It could be moved to
113 // Epetra_Map since it a map operation
114
115 int MaxLocalRows = 0;
116 SubComm_->MaxAll(&NumMyRows_, &MaxLocalRows, 1);
117
118 std::vector<int> SubdomainRowGID(SubComm_->NumProc() * MaxLocalRows);
119 myRowsGID.resize(MaxLocalRows, -1);
120
121 SubComm_->GatherAll(&myRowsGID[0],
122 &SubdomainRowGID[0],
123 MaxLocalRows);
124
125 std::sort(SubdomainRowGID.begin(), SubdomainRowGID.end());
126
127 int PaddingSize = SubdomainRowGID.size() - NumGlobalRows_;
128 SubdomainRowGID.erase(SubdomainRowGID.begin(),
129 SubdomainRowGID.begin() + PaddingSize);
130
131 // Get the col GID corresponding to the process
132 const Epetra_Map &globColMap = Matrix->RowMatrixColMap();
133 NumMyCols_ = globColMap.NumMyElements();
134 std::vector<int> myGlobalCols(NumMyCols_);
135 globColMap.MyGlobalElements(&myGlobalCols[0]);
136
137 // Eliminate column GID that are outside the subdomain
138 std::vector<int> filteredColGID;
139 for (int i = 0; i < NumMyCols_; ++i) {
140 if (std::find(SubdomainRowGID.begin(), SubdomainRowGID.end(),
141 myGlobalCols[i]) != SubdomainRowGID.end()) {
142 filteredColGID.push_back(myGlobalCols[i]);
143 }
144 }
145 NumMyCols_ = filteredColGID.size();
146
147 // Create the maps with the reindexed GIDs
148 Map_ = rcp( new Epetra_Map(-1, NumMyRows_, &myRowsGID[0], globRowMap.IndexBase(), *SubComm_) );
149 colMap_ = rcp( new Epetra_Map(-1, NumMyCols_, &filteredColGID[0], globColMap.IndexBase(), *SubComm_) );
150 NumGlobalCols_ = NumGlobalRows_;
151
152 NumEntries_.resize(NumMyRows_);
153 MaxNumEntriesA_ = Matrix->MaxNumEntries();
154 MaxNumEntries_ = Matrix->MaxNumEntries();
155
156 Diagonal_ = rcp( new Epetra_Vector(*Map_) );
157 if (Diagonal_ == Teuchos::null) IFPACK_CHK_ERRV(-5);
158
159 Indices_.resize(MaxNumEntries_);
160 Values_.resize(MaxNumEntries_);
161
162 Ac_LIDMap_ = 0;
163 Bc_LIDMap_ = 0;
164 Ar_LIDMap_ = 0;
165 Br_LIDMap_ = 0;
166
167 if(ovA_){
168 Ac_LIDMap_ = new int[ovA_->A().NumMyCols() + 1];
169 for(int i = 0; i < ovA_->A().NumMyCols(); i++) {
170 Ac_LIDMap_[i] = colMap_->LID(ovA_->A().RowMatrixColMap().GID(i));
171 }
172 Bc_LIDMap_ = new int[ovA_->B().NumMyCols() + 1];
173 for(int i = 0; i < ovA_->B().NumMyCols(); i++) {
174 Bc_LIDMap_[i] = colMap_->LID(ovA_->B().RowMatrixColMap().GID(i));
175 }
176 Ar_LIDMap_ = new int[ovA_->A().NumMyRows() + 1];
177 for(int i = 0; i < ovA_->A().NumMyRows(); i++) {
178 Ar_LIDMap_[i] = Map_->LID(ovA_->A().RowMatrixRowMap().GID(i));
179 }
180 Br_LIDMap_ = new int[ovA_->B().NumMyRows() + 1];
181 for(int i = 0; i < ovA_->B().NumMyRows(); i++) {
182 Br_LIDMap_[i] = Map_->LID(ovA_->B().RowMatrixRowMap().GID(i));
183 }
184 }
185
186 int ActualMaxNumEntries = 0;
187
188 for (int i = 0 ; i < NumMyRows_; ++i) {
189 NumEntries_[i] = 0;
190 int Nnz;
191 IFPACK_CHK_ERRV(ExtractMyRowCopy(i, MaxNumEntries_, Nnz, &Values_[0], &Indices_[0]));
192
193 for (int j = 0 ; j < Nnz; ++j) {
194 if (Indices_[j] == i) {
195 (*Diagonal_)[i] = Values_[j];
196 }
197 }
198
199 if (Nnz > ActualMaxNumEntries) {
200 ActualMaxNumEntries = Nnz;
201 }
202
203 NumMyNonzeros_ += Nnz;
204 NumEntries_[i] = Nnz;
205 }
206
207 SubComm_->SumAll(&NumMyNonzeros_, &NumGlobalNonzeros_, 1);
208 MaxNumEntries_ = ActualMaxNumEntries;
209
210 int gpid = Matrix->Comm().MyPID();
211 Exporter_ = null;
212 Importer_ = null;
213
214 if (!(RowMatrixRowMap().SameAs(OperatorRangeMap()))) {
215 try{
216 Exporter_ = rcp(new Epetra_Export(RowMatrixRowMap(), OperatorRangeMap()));
217 }
218 catch(...) {
219 printf("** * gpid %d: Ifpack_SubdomainFilter ctor: problem creating Exporter_ * **\n\n",gpid);
220 }
221 }
222
223 if (!(*colMap_).SameAs(*Map_)) {
224 try{
225 Importer_ = rcp(new Epetra_Import(*colMap_, *Map_));
226 }
227 catch(...) {
228 printf("** * gpid %d: Ifpack_SubdomainFilter ctor: problem creating Importer_ * **\n\n",gpid);
229 }
230 }
231
232} //Ifpack_SubdomainFilter() ctor
233
234//==============================================================================
235Ifpack_SubdomainFilter::~Ifpack_SubdomainFilter()
236{
237 if(Ac_LIDMap_) delete [] Ac_LIDMap_;
238 if(Bc_LIDMap_) delete [] Bc_LIDMap_;
239 if(Ar_LIDMap_) delete [] Ar_LIDMap_;
240 if(Br_LIDMap_) delete [] Br_LIDMap_;
241
242 if(ImportVector_) delete ImportVector_;
243} //Ifpack_SubdomainFilter destructor
244
245//==============================================================================
246int Ifpack_SubdomainFilter:: ExtractMyRowCopy(int MyRow, int Length, int & NumEntries,
247 double *Values, int * Indices) const
248{
249 if ((MyRow < 0) || (MyRow >= NumMyRows_)) {
250 IFPACK_CHK_ERR(-1);
251 }
252
253 assert(Length >= NumEntries_[MyRow]);
254
255 int ierr;
256 if (ovA_) {
257 int LocRow = ovA_->A().RowMatrixRowMap().LID(Map_->GID(MyRow));
258 if (LocRow != -1) {
259 ierr = ovA_->A().ExtractMyRowCopy(LocRow, Length, NumEntries, Values, Indices);
260 for(int j = 0;j < NumEntries; j++) {
261 Indices[j] = Ac_LIDMap_[Indices[j]];
262 }
263 }
264 else {
265 LocRow = ovA_->B().RowMatrixRowMap().LID(Map_->GID(MyRow));
266 ierr = ovA_->B().ExtractMyRowCopy(LocRow, Length, NumEntries, Values, Indices);
267 for(int j = 0; j < NumEntries; j++) {
268 Indices[j] = Bc_LIDMap_[Indices[j]];
269 }
270 }
271 } else {
272 int Nnz = 0;
273 ierr = Matrix_->ExtractMyRowCopy(MyRow, MaxNumEntriesA_, Nnz, &Values_[0], &Indices_[0]);
274 IFPACK_CHK_ERR(ierr);
275
276 NumEntries = 0;
277 for (int j = 0 ; j < Nnz ; ++j) {
278 int subdomainLocalIndex = colMap_->LID(Matrix_->RowMatrixColMap().GID(Indices_[j]));
279 if (subdomainLocalIndex != -1) {
280 Indices[NumEntries] = subdomainLocalIndex;
281 Values[NumEntries] = Values_[j];
282 NumEntries++;
283 }
284 }
285 }
286
287 return(ierr);
288}
289
290//==============================================================================
291int Ifpack_SubdomainFilter::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
292{
293 if (!Diagonal.Map().SameAs(*Map_)) {
294 IFPACK_CHK_ERR(-1);
295 }
296
297 Diagonal = *Diagonal_;
298 return(0);
299}
300
301//==============================================================================
302int Ifpack_SubdomainFilter::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
303{
304 std::cout << "Ifpack_SubdomainFilter::Apply not implemented.\n";
305 IFPACK_CHK_ERR(-1); // not implemented
306}
307
308//==============================================================================
309void Ifpack_SubdomainFilter::UpdateImportVector(int NumVectors) const
310{
311}
312
313//==============================================================================
314void Ifpack_SubdomainFilter::UpdateExportVector(int NumVectors) const
315{
316}
317
318//==============================================================================
319int Ifpack_SubdomainFilter::ApplyInverse(const Epetra_MultiVector& X,
320 Epetra_MultiVector& Y) const
321{
322 std::cout << "Ifpack_SubdomainFilter::ApplyInverse not implemented.\n";
323 IFPACK_CHK_ERR(-1); // not implemented
324}
325
326//==============================================================================
327const Epetra_BlockMap& Ifpack_SubdomainFilter::Map() const
328{
329 return(*Map_);
330}
331
332#endif //ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
#define IFPACK_CHK_ERR(ifpack_err)
#define IFPACK_CHK_ERRV(ifpack_err)
int MyGlobalElements(int *MyGlobalElementList) const
int IndexBase() const
bool SameAs(const Epetra_BlockMap &Map) const
int NumMyElements() const
int NumMyRows() const
const Epetra_BlockMap & Map() const
MPI_Comm Comm() const
int MyPID() const
Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.
const int NumVectors
Definition: performance.cpp:71