40 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS 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" 62 Ifpack_SubdomainFilter::Ifpack_SubdomainFilter(
const RefCountPtr<const Epetra_RowMatrix>& Matrix,
int subdomainId) :
71 sprintf(
Label_,
"%s",
"Ifpack_SubdomainFilter");
79 Acrs_=
dynamic_cast<const Epetra_CrsMatrix*
>(&ovA_->A());
80 NumMyRowsA_ = ovA_->
A().NumMyRows();
81 NumMyRowsB_ = ovA_->B().NumMyRows();
83 NumMyRowsA_ = Matrix->NumMyRows();
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_) );
93 SubComm_ = rcp(
new Epetra_SerialComm );
96 NumMyRows_ = Matrix->NumMyRows();
97 SubComm_->SumAll(&NumMyRows_,&NumGlobalRows_,1);
100 const Epetra_Map &globRowMap = Matrix->RowMatrixRowMap();
101 std::vector<int> myRowsGID(NumMyRows_);
102 globRowMap.MyGlobalElements(&myRowsGID[0]);
115 int MaxLocalRows = 0;
116 SubComm_->MaxAll(&NumMyRows_, &MaxLocalRows, 1);
118 std::vector<int> SubdomainRowGID(SubComm_->NumProc() * MaxLocalRows);
119 myRowsGID.resize(MaxLocalRows, -1);
121 SubComm_->GatherAll(&myRowsGID[0],
125 std::sort(SubdomainRowGID.begin(), SubdomainRowGID.end());
127 int PaddingSize = SubdomainRowGID.size() - NumGlobalRows_;
128 SubdomainRowGID.erase(SubdomainRowGID.begin(),
129 SubdomainRowGID.begin() + PaddingSize);
132 const Epetra_Map &globColMap = Matrix->RowMatrixColMap();
133 NumMyCols_ = globColMap.NumMyElements();
134 std::vector<int> myGlobalCols(NumMyCols_);
135 globColMap.MyGlobalElements(&myGlobalCols[0]);
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]);
145 NumMyCols_ = filteredColGID.size();
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_;
156 Diagonal_ = rcp(
new Epetra_Vector(*Map_) );
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));
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));
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));
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));
186 int ActualMaxNumEntries = 0;
188 for (
int i = 0 ; i < NumMyRows_; ++i) {
193 for (
int j = 0 ; j < Nnz; ++j) {
199 if (Nnz > ActualMaxNumEntries) {
200 ActualMaxNumEntries = Nnz;
203 NumMyNonzeros_ += Nnz;
207 SubComm_->SumAll(&NumMyNonzeros_, &NumGlobalNonzeros_, 1);
210 int gpid = Matrix->Comm().MyPID();
219 printf(
"** * gpid %d: Ifpack_SubdomainFilter ctor: problem creating Exporter_ * **\n\n",gpid);
223 if (!(*colMap_).SameAs(*Map_)) {
225 Importer_ = rcp(
new Epetra_Import(*colMap_, *Map_));
228 printf(
"** * gpid %d: Ifpack_SubdomainFilter ctor: problem creating Importer_ * **\n\n",gpid);
235 Ifpack_SubdomainFilter::~Ifpack_SubdomainFilter()
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_;
242 if(ImportVector_)
delete ImportVector_;
246 int Ifpack_SubdomainFilter:: ExtractMyRowCopy(
int MyRow,
int Length,
int & NumEntries,
247 double *Values,
int * Indices)
const 249 if ((MyRow < 0) || (MyRow >= NumMyRows_)) {
253 assert(Length >= NumEntries_[MyRow]);
257 int LocRow = ovA_->A().RowMatrixRowMap().LID(Map_->GID(MyRow));
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]];
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]];
273 ierr = Matrix_->ExtractMyRowCopy(MyRow, MaxNumEntriesA_, Nnz, &Values_[0], &Indices_[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];
291 int Ifpack_SubdomainFilter::ExtractDiagonalCopy(Epetra_Vector & Diagonal)
const 293 if (!Diagonal.Map().SameAs(*Map_)) {
297 Diagonal = *Diagonal_;
302 int Ifpack_SubdomainFilter::Apply(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 304 std::cout <<
"Ifpack_SubdomainFilter::Apply not implemented.\n";
309 void Ifpack_SubdomainFilter::UpdateImportVector(
int NumVectors)
const 314 void Ifpack_SubdomainFilter::UpdateExportVector(
int NumVectors)
const 319 int Ifpack_SubdomainFilter::ApplyInverse(
const Epetra_MultiVector& X,
320 Epetra_MultiVector& Y)
const 322 std::cout <<
"Ifpack_SubdomainFilter::ApplyInverse not implemented.\n";
327 const Epetra_BlockMap& Ifpack_SubdomainFilter::Map()
const 332 #endif //ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
char Label_[80]
Label for this object.
std::vector< int > Indices_
Used in ExtractMyRowCopy, to avoid allocation each time.
virtual const Epetra_Map & RowMatrixRowMap() const
int MaxNumEntries_
Maximum entries in each row.
std::vector< double > Values_
Used in ExtractMyRowCopy, to avoid allocation each time.
#define IFPACK_CHK_ERRV(ifpack_err)
Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.
const Epetra_RowMatrix & A() const
std::vector< int > NumEntries_
#define IFPACK_CHK_ERR(ifpack_err)
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
static bool find(Parser_dh p, char *option, OptionsNode **ptr)
const Epetra_Map & OperatorRangeMap() const