Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_NodeFilter.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#ifdef IFPACK_NODE_AWARE_CODE
44
45#include "Ifpack_ConfigDefs.h"
46
47#include "Epetra_MultiVector.h"
48#include "Epetra_Vector.h"
49#include "Epetra_IntVector.h"
50#include "Epetra_RowMatrix.h"
51#include "Epetra_Map.h"
52#include "Epetra_BlockMap.h"
53#include "Ifpack_NodeFilter.h"
55#include "Epetra_Import.h"
56#include "Epetra_Export.h"
57#include "Epetra_CrsMatrix.h"
58#include "Epetra_BLAS_wrappers.h"
59
60using namespace Teuchos;
61
62#define OLD_AND_BUSTED
63
64//==============================================================================
65Ifpack_NodeFilter::Ifpack_NodeFilter(const RefCountPtr<const Epetra_RowMatrix>& Matrix,int nodeID) :
66 Matrix_(Matrix),
67 NumMyRows_(0),
68 NumMyNonzeros_(0),
69 NumGlobalRows_(0),
70 MaxNumEntries_(0),
71 MaxNumEntriesA_(0)
72{
73 sprintf(Label_,"%s","Ifpack_NodeFilter");
74
75 //ImportVector_=null;
76 //ExportVector_=null;
77 ImportVector_=0;
78 ExportVector_=0;
79
80 ovA_ = dynamic_cast<const Ifpack_OverlappingRowMatrix*>(&*Matrix_);
81 //assert(ovA_ != 0);
82
83 if (ovA_) {
84 Acrs_=dynamic_cast<const Epetra_CrsMatrix*>(&ovA_->A());
85 NumMyRowsA_ = ovA_->A().NumMyRows();
86 NumMyRowsB_ = ovA_->B().NumMyRows();
87 } else {
88 NumMyRowsA_ = Matrix->NumMyRows();
89 NumMyRowsB_ = 0;
90 }
91
92#ifdef HAVE_MPI
93 const Epetra_MpiComm *pComm = dynamic_cast<const Epetra_MpiComm*>( &(Matrix->Comm()) );
94 assert(pComm != NULL);
95 MPI_Comm_split(pComm->Comm(),nodeID,pComm->MyPID(),&nodeMPIComm_);
96 SubComm_ = rcp( new Epetra_MpiComm(nodeMPIComm_) );
97#else
98 SubComm_ = rcp( new Epetra_SerialComm );
99#endif
100
101 NumMyRows_ = Matrix->NumMyRows();
102 SubComm_->SumAll(&NumMyRows_,&NumGlobalRows_,1);
103 NumMyCols_ = Matrix->NumMyCols();
104
105 // build row map, based on the local communicator
106 try {
107 const Epetra_Map &globRowMap = Matrix->RowMatrixRowMap();
108 int *myGlobalElts = globRowMap.MyGlobalElements();
109 int numMyElts = globRowMap.NumMyElements();
110 Map_ = rcp( new Epetra_Map(-1,numMyElts,myGlobalElts,globRowMap.IndexBase(),*SubComm_) );
111 }
112 catch(...) {
113 printf("** * Ifpack_NodeFilter ctor: problem creating row map * **\n\n");
114 }
115
116
117 // build the column map, but don't use a copy constructor, b/c local communicator SubComm_ is
118 // different from that of Matrix.
119 try {
120 const Epetra_Map &globColMap = Matrix->RowMatrixColMap();
121 int *myGlobalElts = globColMap.MyGlobalElements();
122 int numMyElts = globColMap.NumMyElements();
123 colMap_ = rcp( new Epetra_Map(-1,numMyElts,myGlobalElts,globColMap.IndexBase(),*SubComm_) );
124 }
125 catch(...) {
126 printf("** * Ifpack_NodeFilter ctor: problem creating col map * **\n\n");
127 }
128
129 // NumEntries_ will contain the actual number of nonzeros
130 // for each localized row (that is, without external nodes,
131 // and always with the diagonal entry)
132 NumEntries_.resize(NumMyRows_);
133
134 // want to store the diagonal vector. FIXME: am I really useful?
135 Diagonal_ = rcp( new Epetra_Vector(*Map_) );
136 if (Diagonal_ == Teuchos::null) IFPACK_CHK_ERRV(-5);
137
138 // store this for future access to ExtractMyRowCopy().
139 // This is the # of nonzeros in the non-local matrix
140 MaxNumEntriesA_ = Matrix->MaxNumEntries();
141 // tentative value for MaxNumEntries. This is the number of
142 // nonzeros in the local matrix
143 MaxNumEntries_ = Matrix->MaxNumEntries();
144
145 // ExtractMyRowCopy() will use these vectors
146 Indices_.resize(MaxNumEntries_);
147 Values_.resize(MaxNumEntries_);
148
149 // now compute:
150 // - the number of nonzero per row
151 // - the total number of nonzeros
152 // - the diagonal entries
153
154
155 Ac_LIDMap_ = 0;
156 Bc_LIDMap_ = 0;
157 Ar_LIDMap_ = 0;
158 Br_LIDMap_ = 0;
159 tempX_ = 0;
160 tempY_ = 0;
161 // CMS: [A|B]-Local to Overlap-Local Column Indices
162 if(ovA_){
163 Ac_LIDMap_=new int[ovA_->A().NumMyCols()+1];
164 for(int i=0;i<ovA_->A().NumMyCols();i++) Ac_LIDMap_[i]=colMap_->LID(ovA_->A().RowMatrixColMap().GID(i));
165 Bc_LIDMap_=new int[ovA_->B().NumMyCols()+1];
166 for(int i=0;i<ovA_->B().NumMyCols();i++) Bc_LIDMap_[i]=colMap_->LID(ovA_->B().RowMatrixColMap().GID(i));
167
168 Ar_LIDMap_=new int[ovA_->A().NumMyRows()+1];
169 for(int i=0;i<ovA_->A().NumMyRows();i++) Ar_LIDMap_[i]=Map_->LID(ovA_->A().RowMatrixRowMap().GID(i));
170 Br_LIDMap_=new int[ovA_->B().NumMyRows()+1];
171 for(int i=0;i<ovA_->B().NumMyRows();i++) Br_LIDMap_[i]=Map_->LID(ovA_->B().RowMatrixRowMap().GID(i));
172
173
174#ifndef OLD_AND_BUSTED
175 NumMyColsA_=ovA_->A().NumMyCols();
176 tempX_=new double[NumMyColsA_];
177 tempY_=new double[NumMyRowsA_];
178#else
179 tempX_=0;
180 tempY_=0;
181#endif
182 }
183 // end CMS
184
185
186 // compute nonzeros (total and per-row), and store the
187 // diagonal entries (already modified)
188 int ActualMaxNumEntries = 0;
189
190 for (int i = 0 ; i < NumMyRows_ ; ++i) {
191
192 NumEntries_[i] = 0;
193 int Nnz, NewNnz = 0;
194 IFPACK_CHK_ERRV(ExtractMyRowCopy(i,MaxNumEntries_,Nnz,&Values_[0],&Indices_[0]));
195
196 for (int j = 0 ; j < Nnz ; ++j) {
197 NewNnz++;
198 if (Indices_[j] == i) (*Diagonal_)[i] = Values_[j];
199 }
200
201 if (NewNnz > ActualMaxNumEntries)
202 ActualMaxNumEntries = NewNnz;
203
204 NumMyNonzeros_ += NewNnz;
205 NumEntries_[i] = NewNnz;
206
207 }
208
209 SubComm_->SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
210 MaxNumEntries_ = ActualMaxNumEntries;
211
212 int gpid = Matrix->Comm().MyPID();
213 Exporter_ = null;
214 Importer_ = null;
215 // Check if non-trivial import/export operators
216 if (!(RowMatrixRowMap().SameAs(OperatorRangeMap()))) {
217 try{Exporter_ = rcp(new Epetra_Export(RowMatrixRowMap(), OperatorRangeMap()));}
218 catch(...) {
219 printf("** * gpid %d: Ifpack_NodeFilter ctor: problem creating Exporter_ * **\n\n",gpid);
220 }
221 }
222 if (!(*colMap_).SameAs(*Map_)) {
223 //TODO change this to RCP
224 try{Importer_ = rcp(new Epetra_Import(*colMap_, *Map_));}
225 catch(...) {
226 printf("** * gpid %d: Ifpack_NodeFilter ctor: problem creating Importer_ * **\n\n",gpid);
227 }
228 }
229
230} //Ifpack_NodeFilter() ctor
231
232//==============================================================================
233int Ifpack_NodeFilter::
234ExtractMyRowCopy(int MyRow, int Length, int & NumEntries,
235 double *Values, int * Indices) const
236{
237 if ((MyRow < 0) || (MyRow >= NumMyRows_)) {
238 IFPACK_CHK_ERR(-1); // range not valid
239 }
240
241 if (Length < NumEntries_[MyRow])
242 assert(1==0);
243 //return(-1);
244
245 //TODO will we have problems in the original use case?
246 // always extract using the object Values_ and Indices_.
247 // This is because I need more space than that given by
248 // the user (for the external nodes)
249
250 int ierr;
251 if (ovA_) {
252 int LocRow=ovA_->A().RowMatrixRowMap().LID(Map_->GID(MyRow));
253 if (LocRow != -1) {
254 ierr=ovA_->A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
255 for(int j=0;j<NumEntries;j++)
256 Indices[j]=Ac_LIDMap_[Indices[j]];
257
258 }
259 else {
260 LocRow=ovA_->B().RowMatrixRowMap().LID(Map_->GID(MyRow));
261 ierr=ovA_->B().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
262 for(int j=0;j<NumEntries;j++)
263 Indices[j]=Bc_LIDMap_[Indices[j]];
264 }
265 } else {
266 int Nnz;
267 ierr = Matrix_->ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz, &Values_[0],&Indices_[0]);
268 IFPACK_CHK_ERR(ierr);
269
270 NumEntries = 0;
271 for (int j = 0 ; j < Nnz ; ++j) {
272 // only local indices
273 if (Indices_[j] < NumMyRows_ ) {
274 Indices[NumEntries] = Indices_[j];
275 Values[NumEntries] = Values_[j];
276 ++NumEntries;
277 }
278 }
279 }
280
281 return(ierr);
282
283}
284
285//==============================================================================
286int Ifpack_NodeFilter::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
287{
288 if (!Diagonal.Map().SameAs(*Map_))
289 IFPACK_CHK_ERR(-1);
290 Diagonal = *Diagonal_;
291 return(0);
292}
293
294//==============================================================================
295/*
296//old Apply (no communication)
297int Ifpack_NodeFilter::Apply(const Epetra_MultiVector& X,
298 Epetra_MultiVector& Y) const
299{
300
301 // skip expensive checks, I suppose input data are ok
302
303 Y.PutScalar(0.0);
304 int NumVectors = Y.NumVectors();
305
306 double** X_ptr;
307 double** Y_ptr;
308 X.ExtractView(&X_ptr);
309 Y.ExtractView(&Y_ptr);
310
311 for (int i = 0 ; i < NumRows_ ; ++i) {
312
313 int Nnz;
314 int ierr = Matrix_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,&Values_[0],
315 &Indices_[0]);
316 IFPACK_CHK_ERR(ierr);
317
318 for (int j = 0 ; j < Nnz ; ++j) {
319 if (Indices_[j] < NumRows_ ) {
320 for (int k = 0 ; k < NumVectors ; ++k)
321 Y_ptr[k][i] += Values_[j] * X_ptr[k][Indices_[j]];
322 }
323 }
324 }
325
326 return(0);
327}
328*/
329int Ifpack_NodeFilter::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
330 //
331 // This function forms the product Y = A * X.
332 //
333
334 int NumEntries;
335
336 int NumVectors = X.NumVectors();
337 if (NumVectors!=Y.NumVectors()) {
338 EPETRA_CHK_ERR(-1); // Need same number of vectors in each MV
339 }
340
341 // Make sure Import and Export Vectors are compatible
342 UpdateImportVector(NumVectors);
343 UpdateExportVector(NumVectors);
344
345 double ** Xp = (double**) X.Pointers();
346 double ** Yp = (double**) Y.Pointers();
347
348
349 // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
350 if (Importer()!=0) {
351 EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
352 Xp = (double**)ImportVector_->Pointers();
353 }
354
355 // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
356 if (Exporter()!=0) {
357 Yp = (double**)ExportVector_->Pointers();
358 }
359
360 // Do actual computation
361 assert(ovA_!=0);
362 int *MyRows;
363 double *MyValues;
364 int *MyIndices;
365
366 if(Acrs_){
367 // A rows - CrsMatrix Case
368 IFPACK_CHK_ERR(Acrs_->ExtractCrsDataPointers(MyRows,MyIndices,MyValues));
369 //special case NumVectors==1
370 if (NumVectors==1) {
371#ifdef OLD_AND_BUSTED
372
373 for(int i=0;i<NumMyRowsA_;i++) {
374 int LocRow=Ar_LIDMap_[i];
375 double sum = 0.0;
376 for(int j = MyRows[i]; j < MyRows[i+1]; j++)
377 sum += MyValues[j]*Xp[0][Ac_LIDMap_[MyIndices[j]]];
378 Yp[0][LocRow] = sum;
379 }
380#else
381 int izero=0;
382 for(int i=0;i<NumMyColsA_;i++) tempX_[i]=Xp[0][Ac_LIDMap_[i]];
383 EPETRA_DCRSMV_F77(&izero,&NumMyRowsA_,&NumMyRowsA_,MyValues,MyIndices,MyRows,tempX_,tempY_);
384
385 /* for(int i=0;i<NumMyRowsA_;i++) {
386 double sum = 0.0;
387 for(int j = MyRows[i]; j < MyRows[i+1]; j++)
388 sum += MyValues[j]*tempX_[MyIndices[j]];
389 tempY_[i] = sum;
390 }*/
391
392 for(int i=0;i<NumMyRowsA_;i++) Yp[0][Ar_LIDMap_[i]]=tempY_[i];
393#endif
394
395 }
396 else {
397 for(int i=0;i<NumMyRowsA_;i++) {
398 int LocRow=Ar_LIDMap_[i];
399 for (int k=0; k<NumVectors; k++) {
400 double sum = 0.0;
401 for(int j = MyRows[i]; j < MyRows[i+1]; j++)
402 sum += MyValues[j]*Xp[k][Ac_LIDMap_[MyIndices[j]]];
403 Yp[k][LocRow] = sum;
404 }
405 }
406 }
407 }
408 else{
409 // A rows - RowMatrix Case
410 MyValues=&Values_[0];
411 MyIndices=&Indices_[0];
412 for(int i=0;i<NumMyRowsA_;i++) {
413 ovA_->A().ExtractMyRowCopy(i,MaxNumEntries_,NumEntries,MyValues,MyIndices);
414 int LocRow=Ar_LIDMap_[i];
415 for (int k=0; k<NumVectors; k++) {
416 double sum = 0.0;
417 for(int j = 0; j < NumEntries; j++)
418 sum += MyValues[j]*Xp[k][Ac_LIDMap_[MyIndices[j]]];
419 Yp[k][LocRow] = sum;
420 }
421 }
422 }
423
424 // B rows, always CrsMatrix
425 IFPACK_CHK_ERR(ovA_->B().ExtractCrsDataPointers(MyRows,MyIndices,MyValues));
426 //special case NumVectors==1
427 if (NumVectors==1) {
428 for(int i=0;i<NumMyRowsB_;i++) {
429 int LocRow=Br_LIDMap_[i];
430 double sum = 0.0;
431 for(int j = MyRows[i]; j < MyRows[i+1]; j++)
432 sum += MyValues[j]*Xp[0][Bc_LIDMap_[MyIndices[j]]];
433 Yp[0][LocRow] = sum;
434 }
435 } else {
436 for(int i=0;i<NumMyRowsB_;i++) {
437 int LocRow=Br_LIDMap_[i];
438 for (int k=0; k<NumVectors; k++) { //FIXME optimization, check for NumVectors=1
439 double sum = 0.0;
440 for(int j = MyRows[i]; j < MyRows[i+1]; j++)
441 sum += MyValues[j]*Xp[k][Bc_LIDMap_[MyIndices[j]]];
442 Yp[k][LocRow] = sum;
443 }
444 }
445 }
446
447 if (Exporter()!=0) {
448 Y.PutScalar(0.0); // Make sure target is zero
449 Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
450 }
451 // Handle case of rangemap being a local replicated map
452 if (!OperatorRangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
453
454 return(0);
455} //Apply
456
457//==============================================================================
458
459void Ifpack_NodeFilter::UpdateImportVector(int NumVectors) const {
460 if(Importer() != 0) {
461 if(ImportVector_ != 0) {
462 if(ImportVector_->NumVectors() != NumVectors) {
463 delete ImportVector_;
464 ImportVector_= 0;
465 }
466 }
467 if(ImportVector_ == 0)
468 ImportVector_ = new Epetra_MultiVector(Importer_->TargetMap(),NumVectors); // Create Import vector if needed
469 }
470 return;
471/*
472 if(ImportVector_ == null || ImportVector_->NumVectors() != NumVectors)
473 ImportVector_ = rcp(new Epetra_MultiVector(Importer_->TargetMap(),NumVectors));
474*/
475}
476
477//=======================================================================================================
478
479void Ifpack_NodeFilter::UpdateExportVector(int NumVectors) const {
480 if(Exporter() != 0) {
481 if(ExportVector_ != 0) {
482 if(ExportVector_->NumVectors() != NumVectors) {
483 delete ExportVector_;
484 ExportVector_= 0;
485 }
486 }
487 if(ExportVector_ == 0)
488 ExportVector_ = new Epetra_MultiVector(Exporter_->SourceMap(),NumVectors); // Create Export vector if needed
489 }
490 return;
491/*
492 if(ExportVector_ == null || ExportVector_->NumVectors() != NumVectors)
493 ExportVector_ = rcp(new Epetra_MultiVector(Exporter_->SourceMap(),NumVectors));
494*/
495}
496
497//=======================================================================================================
498int Ifpack_NodeFilter::ApplyInverse(const Epetra_MultiVector& X,
499 Epetra_MultiVector& Y) const
500{
501 IFPACK_CHK_ERR(-1); // not implemented
502}
503
504//==============================================================================
505const Epetra_BlockMap& Ifpack_NodeFilter::Map() const
506{
507 return(*Map_);
508
509}
510
511#endif //ifdef IFPACK_NODE_AWARE_CODE
void PREFIX EPETRA_DCRSMV_F77(const int *, const int *, const int *, const double *, const int *, const int *, double *, double *)
#define EPETRA_CHK_ERR(a)
#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
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
MPI_Comm Comm() const
int MyPID() const
int NumVectors() const
double ** Pointers() const
int PutScalar(double ScalarConstant)
Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.
const int NumVectors
Definition: performance.cpp:71