43 #include "Ifpack_ConfigDefs.h"
44 #include "Ifpack_SparsityFilter.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"
53 Ifpack_SparsityFilter::Ifpack_SparsityFilter(
const Teuchos::RefCountPtr<Epetra_RowMatrix>& Matrix,
55 int AllowedBandwidth) :
59 AllowedBandwidth_(AllowedBandwidth),
60 AllowedEntries_(AllowedEntries),
68 if (A_->Comm().NumProc() != 1) {
69 cerr <<
"Ifpack_SparsityFilter can be used with Comm().NumProc() == 1" << endl;
70 cerr <<
"only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
71 cerr <<
"and it is not meant to be used otherwise." << endl;
76 if ((A_->NumMyRows() != A_->NumMyCols()) ||
77 (A_->NumMyRows() != A_->NumGlobalRows64()))
80 NumRows_ = A_->NumMyRows();
81 MaxNumEntriesA_ = A_->MaxNumEntries();
82 Indices_.resize(MaxNumEntriesA_);
83 Values_.resize(MaxNumEntriesA_);
86 if (AllowedBandwidth_ == -1)
87 AllowedBandwidth_ = NumRows_;
92 std::vector<int> Ind(MaxNumEntriesA_);
93 std::vector<double> Val(MaxNumEntriesA_);
95 NumEntries_.resize(NumRows_);
96 for (
int i = 0 ; i < NumRows_ ; ++i)
97 NumEntries_[i] = MaxNumEntriesA_;
99 for (
int i = 0 ; i < A_->NumMyRows() ; ++i) {
101 IFPACK_CHK_ERRV(ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
104 NumEntries_[i] = Nnz;
106 if (Nnz > MaxNumEntries_)
107 MaxNumEntries_ = Nnz;
113 int Ifpack_SparsityFilter::
114 ExtractMyRowCopy(
int MyRow,
int Length,
int & NumEntries,
115 double *Values,
int * Indices)
const
117 if (Length < NumEntries_[MyRow])
121 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz,
122 &Values_[0],&Indices_[0]));
124 double Threshold = 0.0;
127 if (Nnz > AllowedEntries_) {
129 std::vector<double> Values2(Nnz);
131 for (
int i = 0 ; i < Nnz ; ++i) {
133 if (Indices_[i] == MyRow)
136 Values2[count] = IFPACK_ABS(Values_[i]);
140 for (
int i = count ; i < Nnz ; ++i)
144 std::sort(Values2.rbegin(),Values2.rend());
146 Threshold = Values2[AllowedEntries_ - 1];
155 for (
int i = 0 ; i < Nnz ; ++i) {
157 if (IFPACK_ABS(Indices_[i] - MyRow) > AllowedBandwidth_)
160 if ((Indices_[i] != MyRow) && (IFPACK_ABS(Values_[i]) < Threshold))
163 Values[NumEntries] = Values_[i];
164 Indices[NumEntries] = Indices_[i];
167 if (NumEntries > AllowedEntries_)
175 int Ifpack_SparsityFilter::
178 int ierr = A_->ExtractDiagonalCopy(Diagonal);
183 int Ifpack_SparsityFilter::
188 int NumVectors = X.NumVectors();
189 if (NumVectors != Y.NumVectors())
194 std::vector<int> Indices(MaxNumEntries_);
195 std::vector<double> Values(MaxNumEntries_);
197 for (
int i = 0 ; i < A_->NumMyRows() ; ++i) {
200 ExtractMyRowCopy(i,MaxNumEntries_,Nnz,
201 &Values[0], &Indices[0]);
204 for (
int j = 0 ; j < NumVectors ; ++j) {
205 for (
int k = 0 ; k < Nnz ; ++k) {
206 Y[j][i] += Values[k] * X[j][Indices[k]];
212 for (
int j = 0 ; j < NumVectors ; ++j) {
213 for (
int k = 0 ; k < Nnz ; ++k) {
214 Y[j][Indices[k]] += Values[k] * X[j][i];
224 int Ifpack_SparsityFilter::
225 Solve(
bool ,
bool ,
bool ,
232 int Ifpack_SparsityFilter::
235 int ierr = Multiply(UseTranspose(),X,Y);
240 int Ifpack_SparsityFilter::