43 #ifdef IFPACK_NODE_AWARE_CODE
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"
55 #include "Epetra_Import.h"
56 #include "Epetra_Export.h"
57 #include "Epetra_CrsMatrix.h"
58 #include "Epetra_BLAS_wrappers.h"
60 using namespace Teuchos;
62 #define OLD_AND_BUSTED
65 Ifpack_NodeFilter::Ifpack_NodeFilter(
const RefCountPtr<const Epetra_RowMatrix>& Matrix,
int nodeID) :
73 sprintf(Label_,
"%s",
"Ifpack_NodeFilter");
86 NumMyRowsB_ = ovA_->B().NumMyRows();
88 NumMyRowsA_ = Matrix->NumMyRows();
94 assert(pComm != NULL);
95 MPI_Comm_split(pComm->
Comm(),nodeID,pComm->
MyPID(),&nodeMPIComm_);
101 NumMyRows_ = Matrix->NumMyRows();
102 SubComm_->SumAll(&NumMyRows_,&NumGlobalRows_,1);
103 NumMyCols_ = Matrix->NumMyCols();
107 const Epetra_Map &globRowMap = Matrix->RowMatrixRowMap();
113 printf(
"** * Ifpack_NodeFilter ctor: problem creating row map * **\n\n");
120 const Epetra_Map &globColMap = Matrix->RowMatrixColMap();
126 printf(
"** * Ifpack_NodeFilter ctor: problem creating col map * **\n\n");
132 NumEntries_.resize(NumMyRows_);
140 MaxNumEntriesA_ = Matrix->MaxNumEntries();
143 MaxNumEntries_ = Matrix->MaxNumEntries();
146 Indices_.resize(MaxNumEntries_);
147 Values_.resize(MaxNumEntries_);
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));
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));
174 #ifndef OLD_AND_BUSTED
175 NumMyColsA_=ovA_->A().NumMyCols();
176 tempX_=
new double[NumMyColsA_];
177 tempY_=
new double[NumMyRowsA_];
188 int ActualMaxNumEntries = 0;
190 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
194 IFPACK_CHK_ERRV(ExtractMyRowCopy(i,MaxNumEntries_,Nnz,&Values_[0],&Indices_[0]));
196 for (
int j = 0 ; j < Nnz ; ++j) {
198 if (Indices_[j] == i) (*Diagonal_)[i] = Values_[j];
201 if (NewNnz > ActualMaxNumEntries)
202 ActualMaxNumEntries = NewNnz;
204 NumMyNonzeros_ += NewNnz;
205 NumEntries_[i] = NewNnz;
209 SubComm_->SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
210 MaxNumEntries_ = ActualMaxNumEntries;
212 int gpid = Matrix->Comm().MyPID();
216 if (!(RowMatrixRowMap().SameAs(OperatorRangeMap()))) {
217 try{Exporter_ =
rcp(
new Epetra_Export(RowMatrixRowMap(), OperatorRangeMap()));}
219 printf(
"** * gpid %d: Ifpack_NodeFilter ctor: problem creating Exporter_ * **\n\n",gpid);
222 if (!(*colMap_).SameAs(*Map_)) {
226 printf(
"** * gpid %d: Ifpack_NodeFilter ctor: problem creating Importer_ * **\n\n",gpid);
233 int Ifpack_NodeFilter::
234 ExtractMyRowCopy(
int MyRow,
int Length,
int & NumEntries,
235 double *Values,
int * Indices)
const
237 if ((MyRow < 0) || (MyRow >= NumMyRows_)) {
241 if (Length < NumEntries_[MyRow])
252 int LocRow=ovA_->A().RowMatrixRowMap().LID(Map_->GID(MyRow));
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]];
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]];
267 ierr = Matrix_->ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz, &Values_[0],&Indices_[0]);
271 for (
int j = 0 ; j < Nnz ; ++j) {
273 if (Indices_[j] < NumMyRows_ ) {
274 Indices[NumEntries] = Indices_[j];
275 Values[NumEntries] = Values_[j];
286 int Ifpack_NodeFilter::ExtractDiagonalCopy(
Epetra_Vector & Diagonal)
const
290 Diagonal = *Diagonal_;
337 if (NumVectors!=Y.NumVectors()) {
342 UpdateImportVector(NumVectors);
343 UpdateExportVector(NumVectors);
345 double ** Xp = (
double**) X.Pointers();
346 double ** Yp = (
double**) Y.Pointers();
352 Xp = (
double**)ImportVector_->Pointers();
357 Yp = (
double**)ExportVector_->Pointers();
368 IFPACK_CHK_ERR(Acrs_->ExtractCrsDataPointers(MyRows,MyIndices,MyValues));
371 #ifdef OLD_AND_BUSTED
373 for(
int i=0;i<NumMyRowsA_;i++) {
374 int LocRow=Ar_LIDMap_[i];
376 for(
int j = MyRows[i]; j < MyRows[i+1]; j++)
377 sum += MyValues[j]*Xp[0][Ac_LIDMap_[MyIndices[j]]];
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_);
392 for(
int i=0;i<NumMyRowsA_;i++) Yp[0][Ar_LIDMap_[i]]=tempY_[i];
397 for(
int i=0;i<NumMyRowsA_;i++) {
398 int LocRow=Ar_LIDMap_[i];
401 for(
int j = MyRows[i]; j < MyRows[i+1]; j++)
402 sum += MyValues[j]*Xp[k][Ac_LIDMap_[MyIndices[j]]];
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];
417 for(
int j = 0; j < NumEntries; j++)
418 sum += MyValues[j]*Xp[k][Ac_LIDMap_[MyIndices[j]]];
425 IFPACK_CHK_ERR(ovA_->B().ExtractCrsDataPointers(MyRows,MyIndices,MyValues));
428 for(
int i=0;i<NumMyRowsB_;i++) {
429 int LocRow=Br_LIDMap_[i];
431 for(
int j = MyRows[i]; j < MyRows[i+1]; j++)
432 sum += MyValues[j]*Xp[0][Bc_LIDMap_[MyIndices[j]]];
436 for(
int i=0;i<NumMyRowsB_;i++) {
437 int LocRow=Br_LIDMap_[i];
440 for(
int j = MyRows[i]; j < MyRows[i+1]; j++)
441 sum += MyValues[j]*Xp[k][Bc_LIDMap_[MyIndices[j]]];
449 Y.Export(*ExportVector_, *Exporter(), Add);
452 if (!OperatorRangeMap().DistributedGlobal() &&
Comm().NumProc()>1)
EPETRA_CHK_ERR(Y.Reduce());
459 void Ifpack_NodeFilter::UpdateImportVector(
int NumVectors)
const {
460 if(Importer() != 0) {
461 if(ImportVector_ != 0) {
462 if(ImportVector_->NumVectors() !=
NumVectors) {
463 delete ImportVector_;
467 if(ImportVector_ == 0)
479 void Ifpack_NodeFilter::UpdateExportVector(
int NumVectors)
const {
480 if(Exporter() != 0) {
481 if(ExportVector_ != 0) {
482 if(ExportVector_->NumVectors() !=
NumVectors) {
483 delete ExportVector_;
487 if(ExportVector_ == 0)
511 #endif //ifdef IFPACK_NODE_AWARE_CODE
bool SameAs(const Epetra_BlockMap &Map) const
int MyGlobalElements(int *MyGlobalElementList) const
#define EPETRA_CHK_ERR(a)
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
Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.
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)