Ifpack_SparsityFilter.cpp

00001 #include "Ifpack_ConfigDefs.h"
00002 #include "Ifpack_SparsityFilter.h"
00003 #include "Epetra_ConfigDefs.h"
00004 #include "Epetra_RowMatrix.h"
00005 #include "Epetra_Comm.h"
00006 #include "Epetra_Map.h"
00007 #include "Epetra_MultiVector.h"
00008 #include "Epetra_Vector.h"
00009 
00010 //==============================================================================
00011 Ifpack_SparsityFilter::Ifpack_SparsityFilter(Epetra_RowMatrix* Matrix,
00012                          int AllowedEntries, 
00013                          int AllowedBandwidth) :
00014   A_(*Matrix),
00015   MaxNumEntries_(0),
00016   MaxNumEntriesA_(0),
00017   AllowedBandwidth_(AllowedBandwidth),
00018   AllowedEntries_(AllowedEntries),
00019   NumNonzeros_(0),
00020   NumRows_(0)
00021 {
00022   // use this filter only on serial matrices
00023   if (A_.Comm().NumProc() != 1) {
00024     cerr << "Ifpack_DropFilter can be used with Comm().NumProc() == 1" << endl;
00025     cerr << "only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
00026     cerr << "and it is not meant to be used otherwise." << endl;
00027     exit(EXIT_FAILURE);
00028   }
00029 
00030   // only square serial matrices
00031   if ((A_.NumMyRows() != A_.NumMyCols()) ||
00032      (A_.NumMyRows() != A_.NumGlobalRows()))
00033     IFPACK_CHK_ERRV(-1);
00034 
00035   NumRows_ = A_.NumMyRows();
00036   MaxNumEntriesA_ = A_.MaxNumEntries();
00037   Indices_.resize(MaxNumEntriesA_);
00038   Values_.resize(MaxNumEntriesA_);
00039 
00040   // default value is to not consider bandwidth
00041   if (AllowedBandwidth_ == -1)
00042     AllowedBandwidth_ = NumRows_;
00043   
00044   // computes the number of nonzero elements per row in the 
00045   // dropped matrix. Stores this number in NumEntries_.
00046   // Also, computes the global number of nonzeros.
00047   vector<int>    Ind(MaxNumEntriesA_);
00048   vector<double> Val(MaxNumEntriesA_);
00049 
00050   NumEntries_.resize(NumRows_);
00051   for (int i = 0 ; i < NumRows_ ; ++i)
00052     NumEntries_[i] = MaxNumEntriesA_;
00053 
00054   for (int i = 0 ; i < A_.NumMyRows() ; ++i) {
00055     int Nnz;
00056     IFPACK_CHK_ERRV(ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00057                      &Val[0], &Ind[0]));
00058 
00059     NumEntries_[i] = Nnz;
00060     NumNonzeros_ += Nnz;
00061     if (Nnz > MaxNumEntries_)
00062       MaxNumEntries_ = Nnz;
00063 
00064   }
00065 }
00066 
00067 //==============================================================================
00068 Ifpack_SparsityFilter::~Ifpack_SparsityFilter()
00069 {
00070 }
00071 
00072 //==============================================================================
00073 int Ifpack_SparsityFilter::
00074 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00075          double *Values, int * Indices) const
00076 {
00077   if (Length < NumEntries_[MyRow])
00078     IFPACK_CHK_ERR(-1);
00079 
00080   int Nnz;
00081   IFPACK_CHK_ERR(A_.ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz,
00082                      &Values_[0],&Indices_[0]));
00083 
00084   double Threshold = 0.0;
00085     
00086   // this `if' is to define the cut-off value
00087   if (Nnz > AllowedEntries_) {
00088  
00089     vector<double> Values2(Nnz);
00090     int count = 0;
00091     for (int i = 0 ; i < Nnz ; ++i) {
00092       // skip diagonal entry (which is always inserted)
00093       if (Indices_[i] == MyRow)
00094     continue;
00095       // put absolute value
00096       Values2[count] = IFPACK_ABS(Values_[i]);
00097       count++;
00098     }
00099 
00100     for (int i = count ; i < Nnz ; ++i)
00101       Values2[i] = 0.0;
00102 
00103     // sort in descending order
00104     sort(Values2.rbegin(),Values2.rend());
00105     // get the cut-off value
00106     Threshold = Values2[AllowedEntries_ - 1];
00107 
00108   }
00109 
00110   // loop over all nonzero elements of row MyRow,
00111   // and drop elements below specified threshold.
00112   // Also, add value to zero diagonal
00113   NumEntries = 0;
00114 
00115   for (int i = 0 ; i < Nnz ; ++i) {
00116 
00117     if (IFPACK_ABS(Indices_[i] - MyRow) > AllowedBandwidth_)
00118       continue;
00119 
00120     if ((Indices_[i] != MyRow) && (IFPACK_ABS(Values_[i]) < Threshold))
00121       continue;
00122 
00123     Values[NumEntries] = Values_[i];
00124     Indices[NumEntries] = Indices_[i];
00125 
00126     NumEntries++;
00127     if (NumEntries > AllowedEntries_)
00128       break;
00129   }
00130 
00131   return(0);
00132 }
00133 
00134 //==============================================================================
00135 int Ifpack_SparsityFilter::
00136 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00137 {
00138   IFPACK_RETURN(A_.ExtractDiagonalCopy(Diagonal));
00139 }
00140 
00141 //==============================================================================
00142 int Ifpack_SparsityFilter::
00143 Multiply(bool TransA, const Epetra_MultiVector& X, 
00144      Epetra_MultiVector& Y) const
00145 {
00146 
00147   int NumVectors = X.NumVectors();
00148   if (NumVectors != Y.NumVectors())
00149     IFPACK_CHK_ERR(-1);
00150 
00151   Y.PutScalar(0.0);
00152 
00153   vector<int> Indices(MaxNumEntries_);
00154   vector<double> Values(MaxNumEntries_);
00155 
00156   for (int i = 0 ; i < A_.NumMyRows() ; ++i) {
00157 
00158     int Nnz;
00159     ExtractMyRowCopy(i,MaxNumEntries_,Nnz,
00160              &Values[0], &Indices[0]);
00161     if (!TransA) {
00162       // no transpose first
00163       for (int j = 0 ; j < NumVectors ; ++j) {
00164     for (int k = 0 ; k < Nnz ; ++k) {
00165       Y[j][i] += Values[k] * X[j][Indices[k]];
00166     }
00167       }
00168     }
00169     else {
00170       // transpose here
00171       for (int j = 0 ; j < NumVectors ; ++j) {
00172     for (int k = 0 ; k < Nnz ; ++k) {
00173       Y[j][Indices[k]] += Values[k] * X[j][i];
00174     }
00175       }
00176     }
00177   }
00178 
00179   return(0);
00180 }
00181 
00182 //==============================================================================
00183 int Ifpack_SparsityFilter::
00184 Solve(bool Upper, bool Trans, bool UnitDiagonal, 
00185       const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00186 {
00187   IFPACK_CHK_ERR(-98);
00188 }
00189 
00190 //==============================================================================
00191 int Ifpack_SparsityFilter::
00192 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00193 {
00194   IFPACK_RETURN(Multiply(UseTranspose(),X,Y));
00195 }
00196 
00197 //==============================================================================
00198 int Ifpack_SparsityFilter::
00199 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00200 {
00201   IFPACK_CHK_ERR(-98); 
00202 }

Generated on Thu Sep 18 12:37:08 2008 for IFPACK by doxygen 1.3.9.1