43 #include "Ifpack_ConfigDefs.h"
44 #include "Ifpack_SingletonFilter.h"
45 #include "Epetra_ConfigDefs.h"
46 #include "Epetra_RowMatrix.h"
47 #include "Epetra_Comm.h"
48 #include "Epetra_Map.h"
49 #include "Epetra_MultiVector.h"
50 #include "Epetra_Vector.h"
51 #include "Epetra_Import.h"
67 if (A_->Comm().NumProc() != 1) {
68 cerr <<
"Ifpack_SingletonFilter can be used with Comm().NumProc() == 1" << endl;
69 cerr <<
"only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
70 cerr <<
"and it is not meant to be used otherwise." << endl;
74 if ((A_->NumMyRows() != A_->NumGlobalRows64()) ||
75 (A_->NumMyRows() != A_->NumMyCols()))
78 NumRowsA_ = (A_->NumMyRows());
79 MaxNumEntriesA_ = A_->MaxNumEntries();
81 Indices_.resize(MaxNumEntriesA_);
82 Values_.resize(MaxNumEntriesA_);
83 Reorder_.resize(A_->NumMyRows());
85 for (
int i = 0 ; i < NumRowsA_ ; ++i)
89 for (
int i = 0 ; i < NumRowsA_ ; ++i) {
91 IFPACK_CHK_ERRV(A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
92 &Values_[0], &Indices_[0]));
94 Reorder_[i] = NumRows_++;
101 InvReorder_.resize(NumRows_);
102 for (
int i = 0 ; i < NumRowsA_ ; ++i) {
105 InvReorder_[Reorder_[i]] = i;
108 NumEntries_.resize(NumRows_);
109 SingletonIndex_.resize(NumSingletons_);
113 for (
int i = 0 ; i < A_->NumMyRows() ; ++i) {
116 IFPACK_CHK_ERRV(A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
117 &Values_[0], &Indices_[0]));
118 int ii = Reorder_[i];
122 NumEntries_[ii] = Nnz;
124 if (Nnz > MaxNumEntries_)
125 MaxNumEntries_ = Nnz;
128 SingletonIndex_[count] = i;
133 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
134 Map_ = Teuchos::rcp(
new Epetra_Map(NumRows_,0,Comm()) );
141 A_->ExtractDiagonalCopy(Diagonal);
142 for (
int i = 0 ; i < NumRows_ ; ++i) {
143 int ii = InvReorder_[i];
144 (*Diagonal_)[i] = Diagonal[ii];
150 int Ifpack_SingletonFilter::
151 ExtractMyRowCopy(
int MyRow,
int Length,
int & NumEntries,
152 double *Values,
int * Indices)
const
156 if (Length < NumEntries_[MyRow])
159 int Row = InvReorder_[MyRow];
160 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(Row,MaxNumEntriesA_,Nnz,
161 &Values_[0],&Indices_[0]));
163 for (
int i = 0 ; i < Nnz ; ++i) {
164 int ii = Reorder_[Indices_[i]];
166 Indices[NumEntries] = ii;
167 Values[NumEntries] = Values_[i];
175 int Ifpack_SingletonFilter::
178 Diagonal = *Diagonal_;
183 int Ifpack_SingletonFilter::
188 int NumVectors = X.NumVectors();
189 if (NumVectors != Y.NumVectors())
194 std::vector<int> Indices(MaxNumEntries_);
195 std::vector<double> Values(MaxNumEntries_);
198 for (
int i = 0 ; i < A_->NumMyRows() ; ++i) {
204 A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
205 &Values[0], &Indices[0]);
208 for (
int j = 0 ; j < NumVectors ; ++j) {
209 for (
int k = 0 ; k < Nnz ; ++k) {
210 if (Reorder_[i] >= 0)
211 Y[j][i] += Values[k] * X[j][Reorder_[Indices[k]]];
217 for (
int j = 0 ; j < NumVectors ; ++j) {
218 for (
int k = 0 ; k < Nnz ; ++k) {
219 if (Reorder_[i] >= 0)
220 Y[j][Reorder_[Indices[k]]] += Values[k] * X[j][i];
230 int Ifpack_SingletonFilter::
231 Solve(
bool ,
bool ,
bool ,
238 int Ifpack_SingletonFilter::
241 IFPACK_CHK_ERR(Multiply(
false,X,Y));
246 int Ifpack_SingletonFilter::
253 int Ifpack_SingletonFilter::
257 for (
int i = 0 ; i < NumSingletons_ ; ++i) {
258 int ii = SingletonIndex_[i];
261 A_->ExtractMyRowCopy(ii,MaxNumEntriesA_,Nnz,
262 &Values_[0], &Indices_[0]);
263 for (
int j = 0 ; j < Nnz ; ++j) {
264 if (Indices_[j] == ii) {
265 for (
int k = 0 ; k < LHS.NumVectors() ; ++k)
266 LHS[k][ii] = RHS[k][ii] / Values_[j];
275 int Ifpack_SingletonFilter::
280 int NumVectors = LHS.NumVectors();
282 for (
int i = 0 ; i < NumRows_ ; ++i)
283 for (
int k = 0 ; k < NumVectors ; ++k)
284 ReducedRHS[k][i] = RHS[k][InvReorder_[i]];
286 for (
int i = 0 ; i < NumRows_ ; ++i) {
287 int ii = InvReorder_[i];
289 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(ii,MaxNumEntriesA_,Nnz,
290 &Values_[0], &Indices_[0]));
292 for (
int j = 0 ; j < Nnz ; ++j) {
293 if (Reorder_[Indices_[j]] == -1) {
294 for (
int k = 0 ; k < NumVectors ; ++k)
295 ReducedRHS[k][i] -= Values_[j] * LHS[k][Indices_[j]];
303 int Ifpack_SingletonFilter::
307 for (
int i = 0 ; i < NumRows_ ; ++i)
308 for (
int k = 0 ; k < LHS.NumVectors() ; ++k)
309 LHS[k][InvReorder_[i]] = ReducedLHS[k][i];
Ifpack_SingletonFilter(const Teuchos::RefCountPtr< Epetra_RowMatrix > &Matrix)
Constructor.