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"
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();
150 NumGlobalCols_ = NumGlobalRows_;
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_)) {
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
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
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
char Label_[80]
Label for this object.
std::vector< int > Indices_
Used in ExtractMyRowCopy, to avoid allocation each time.
bool SameAs(const Epetra_BlockMap &Map) const
virtual const Epetra_Map & RowMatrixRowMap() const
int MyGlobalElements(int *MyGlobalElementList) const
int MaxNumEntries_
Maximum entries in each row.
const Epetra_Map & OperatorRangeMap() const
std::vector< double > Values_
Used in ExtractMyRowCopy, to avoid allocation each time.
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.
std::vector< int > NumEntries_
#define IFPACK_CHK_ERR(ifpack_err)
static bool find(Parser_dh p, char *option, OptionsNode **ptr)