Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
60 using namespace Teuchos;
61 
62 #define OLD_AND_BUSTED
63 
64 //==============================================================================
65 Ifpack_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 //==============================================================================
233 int Ifpack_NodeFilter::
234 ExtractMyRowCopy(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 //==============================================================================
286 int 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)
297 int 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 */
329 int 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 
459 void 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 
479 void 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 //=======================================================================================================
498 int Ifpack_NodeFilter::ApplyInverse(const Epetra_MultiVector& X,
499  Epetra_MultiVector& Y) const
500 {
501  IFPACK_CHK_ERR(-1); // not implemented
502 }
503 
504 //==============================================================================
505 const Epetra_BlockMap& Ifpack_NodeFilter::Map() const
506 {
507  return(*Map_);
508 
509 }
510 
511 #endif //ifdef IFPACK_NODE_AWARE_CODE
bool SameAs(const Epetra_BlockMap &Map) const
int MyGlobalElements(int *MyGlobalElementList) const
#define EPETRA_CHK_ERR(a)
const int NumVectors
Definition: performance.cpp:71
int IndexBase() const
int NumMyElements() const
#define IFPACK_CHK_ERRV(ifpack_err)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual const Epetra_BlockMap & Map() const =0
int NumMyRows() const
Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.
MPI_Comm Comm() const
void PREFIX EPETRA_DCRSMV_F77(const int *, const int *, const int *, const double *, const int *, const int *, double *, double *)
#define IFPACK_CHK_ERR(ifpack_err)
int MyPID() const