38 #include "Ifpack_ConfigDefs.h"
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"
51 #include "Ifpack_OverlappingRowMatrix.h"
52 #include "Epetra_Import.h"
53 #include "Epetra_Export.h"
54 #include "Epetra_CrsMatrix.h"
55 #include "Epetra_BLAS_wrappers.h"
57 #include "Ifpack_SubdomainFilter.h"
59 using namespace Teuchos;
62 Ifpack_SubdomainFilter::Ifpack_SubdomainFilter(
const RefCountPtr<const Epetra_RowMatrix>& Matrix,
int subdomainId) :
71 sprintf(Label_,
"%s",
"Ifpack_SubdomainFilter");
81 NumMyRowsB_ = ovA_->B().NumMyRows();
83 NumMyRowsA_ = Matrix->NumMyRows();
89 assert(pComm != NULL);
90 MPI_Comm_split(pComm->
Comm(), subdomainId, pComm->
MyPID(), &subdomainMPIComm_);
96 NumMyRows_ = Matrix->NumMyRows();
97 SubComm_->SumAll(&NumMyRows_,&NumGlobalRows_,1);
100 const Epetra_Map &globRowMap = Matrix->RowMatrixRowMap();
101 std::vector<int> myRowsGID(NumMyRows_);
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();
134 std::vector<int> myGlobalCols(NumMyCols_);
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_;
152 NumEntries_.resize(NumMyRows_);
153 MaxNumEntriesA_ = Matrix->MaxNumEntries();
154 MaxNumEntries_ = Matrix->MaxNumEntries();
157 if (Diagonal_ == Teuchos::null) IFPACK_CHK_ERRV(-5);
159 Indices_.resize(MaxNumEntries_);
160 Values_.resize(MaxNumEntries_);
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) {
191 IFPACK_CHK_ERRV(ExtractMyRowCopy(i, MaxNumEntries_, Nnz, &Values_[0], &Indices_[0]));
193 for (
int j = 0 ; j < Nnz; ++j) {
194 if (Indices_[j] == i) {
195 (*Diagonal_)[i] = Values_[j];
199 if (Nnz > ActualMaxNumEntries) {
200 ActualMaxNumEntries = Nnz;
203 NumMyNonzeros_ += Nnz;
204 NumEntries_[i] = Nnz;
207 SubComm_->SumAll(&NumMyNonzeros_, &NumGlobalNonzeros_, 1);
208 MaxNumEntries_ = ActualMaxNumEntries;
210 int gpid = Matrix->Comm().MyPID();
214 if (!(RowMatrixRowMap().SameAs(OperatorRangeMap()))) {
216 Exporter_ = rcp(
new Epetra_Export(RowMatrixRowMap(), OperatorRangeMap()));
219 printf(
"** * gpid %d: Ifpack_SubdomainFilter ctor: problem creating Exporter_ * **\n\n",gpid);
223 if (!(*colMap_).SameAs(*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]);
274 IFPACK_CHK_ERR(ierr);
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_;
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
322 std::cout <<
"Ifpack_SubdomainFilter::ApplyInverse not implemented.\n";
332 #endif //ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
int MyGlobalElements(int *MyGlobalElementList) const
int NumMyElements() const
Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.