43 #ifdef IFPACK_NODE_AWARE_CODE 45 #include "Ifpack_ConfigDefs.h" 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" 54 #include "Ifpack_OverlappingRowMatrix.h" 55 #include "Epetra_Import.h" 56 #include "Epetra_Export.h" 57 #include "Epetra_CrsMatrix.h" 58 #include "Epetra_BLAS_wrappers.h" 62 #define OLD_AND_BUSTED 65 Ifpack_NodeFilter::Ifpack_NodeFilter(
const RefCountPtr<const Epetra_RowMatrix>& Matrix,
int nodeID) :
73 sprintf(Label_,
"%s",
"Ifpack_NodeFilter");
84 Acrs_=
dynamic_cast<const Epetra_CrsMatrix*
>(&ovA_->A());
85 NumMyRowsA_ = ovA_->A().NumMyRows();
86 NumMyRowsB_ = ovA_->B().NumMyRows();
88 NumMyRowsA_ = Matrix->NumMyRows();
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_) );
98 SubComm_ = rcp(
new Epetra_SerialComm );
101 NumMyRows_ = Matrix->NumMyRows();
102 SubComm_->SumAll(&NumMyRows_,&NumGlobalRows_,1);
103 NumMyCols_ = Matrix->NumMyCols();
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_) );
113 printf(
"** * Ifpack_NodeFilter ctor: problem creating row map * **\n\n");
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_) );
126 printf(
"** * Ifpack_NodeFilter ctor: problem creating col map * **\n\n");
132 NumEntries_.resize(NumMyRows_);
135 Diagonal_ = rcp(
new Epetra_Vector(*Map_) );
136 if (Diagonal_ == Teuchos::null) IFPACK_CHK_ERRV(-5);
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_)) {
224 try{Importer_ = rcp(
new Epetra_Import(*colMap_, *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]);
268 IFPACK_CHK_ERR(ierr);
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 288 if (!Diagonal.Map().SameAs(*Map_))
290 Diagonal = *Diagonal_;
329 int Ifpack_NodeFilter::Apply(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const {
336 int NumVectors = X.NumVectors();
337 if (NumVectors!=Y.NumVectors()) {
342 UpdateImportVector(NumVectors);
343 UpdateExportVector(NumVectors);
345 double ** Xp = (
double**) X.Pointers();
346 double ** Yp = (
double**) Y.Pointers();
351 EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
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];
399 for (
int k=0; k<NumVectors; k++) {
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];
415 for (
int k=0; k<NumVectors; k++) {
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];
438 for (
int k=0; k<NumVectors; k++) {
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)
468 ImportVector_ =
new Epetra_MultiVector(Importer_->TargetMap(),NumVectors);
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)
488 ExportVector_ =
new Epetra_MultiVector(Exporter_->SourceMap(),NumVectors);
498 int Ifpack_NodeFilter::ApplyInverse(
const Epetra_MultiVector& X,
499 Epetra_MultiVector& Y)
const 505 const Epetra_BlockMap& Ifpack_NodeFilter::Map()
const 511 #endif //ifdef IFPACK_NODE_AWARE_CODE Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.