Ifpack_DropFilter.cpp

Go to the documentation of this file.
00001 #include "Ifpack_ConfigDefs.h"
00002 #include "Ifpack_DropFilter.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_DropFilter::Ifpack_DropFilter(Epetra_RowMatrix* Matrix,
00012              double DropTol) :
00013   A_(*Matrix),
00014   DropTol_(DropTol),
00015   MaxNumEntries_(0),
00016   MaxNumEntriesA_(0),
00017   NumNonzeros_(0)
00018 {
00019   // use this filter only on serial matrices
00020   if (A_.Comm().NumProc() != 1) {
00021     cerr << "Ifpack_DropFilter can be used with Comm().NumProc() == 1" << endl;
00022     cerr << "only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
00023     cerr << "and it is not meant to be used otherwise." << endl;
00024     exit(EXIT_FAILURE);
00025   }
00026   
00027   if ((A_.NumMyRows() != A_.NumGlobalRows()) ||
00028       (A_.NumMyRows() != A_.NumMyCols()))
00029     IFPACK_CHK_ERRV(-2);
00030 
00031   NumRows_ = A_.NumMyRows();
00032   MaxNumEntriesA_ = A_.MaxNumEntries();
00033 
00034   NumEntries_.resize(NumRows_);
00035   Indices_.resize(MaxNumEntriesA_);
00036   Values_.resize(MaxNumEntriesA_);
00037 
00038   vector<int>    Ind(MaxNumEntriesA_);
00039   vector<double> Val(MaxNumEntriesA_);
00040 
00041   for (int i = 0 ; i < NumRows_ ; ++i) {
00042     NumEntries_[i] = MaxNumEntriesA_;
00043     int Nnz;
00044     IFPACK_CHK_ERRV(ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00045              &Val[0], &Ind[0]));
00046  
00047     NumEntries_[i] = Nnz;
00048     NumNonzeros_ += Nnz;
00049     if (Nnz > MaxNumEntries_)
00050       MaxNumEntries_ = Nnz;
00051   }
00052 
00053 }
00054 
00055 //==============================================================================
00056 Ifpack_DropFilter::~Ifpack_DropFilter()
00057 {
00058 }
00059 
00060 //==============================================================================
00061 int Ifpack_DropFilter::
00062 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00063      double *Values, int * Indices) const
00064 {
00065   if (Length < NumEntries_[MyRow])
00066     IFPACK_CHK_ERR(-1);
00067 
00068   int Nnz;
00069 
00070   IFPACK_CHK_ERR(A_.ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz,
00071              &Values_[0],&Indices_[0]));
00072 
00073   // loop over all nonzero elements of row MyRow,
00074   // and drop elements below specified threshold.
00075   // Also, add value to zero diagonal
00076   int count = 0;
00077   for (int i = 0 ; i < Nnz ; ++i) {
00078 
00079     // if element is above specified tol, add to the
00080     // user's defined arrays. Check that we are not
00081     // exceeding allocated space. Do not drop any diagonal entry.
00082     if ((Indices_[i] == MyRow) || (IFPACK_ABS(Values_[i]) >= DropTol_)) {
00083       if (count == Length)
00084   IFPACK_CHK_ERR(-1);
00085       Values[count] = Values_[i];
00086       Indices[count] = Indices_[i];
00087       count++;
00088     }
00089   }
00090 
00091   NumEntries = count;
00092   return(0);
00093 }
00094 
00095 //==============================================================================
00096 int Ifpack_DropFilter::
00097 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00098 {
00099   IFPACK_CHK_ERR(A_.ExtractDiagonalCopy(Diagonal));
00100   return(0);
00101 }
00102 
00103 //==============================================================================
00104 int Ifpack_DropFilter::
00105 Multiply(bool TransA, const Epetra_MultiVector& X, 
00106    Epetra_MultiVector& Y) const
00107 {
00108   // NOTE: I suppose that the matrix has been localized,
00109   // hence all maps are trivial.
00110   int NumVectors = X.NumVectors();
00111   if (NumVectors != Y.NumVectors())
00112     IFPACK_CHK_ERR(-1);
00113 
00114   Y.PutScalar(0.0);
00115 
00116   vector<int> Indices(MaxNumEntries_);
00117   vector<double> Values(MaxNumEntries_);
00118 
00119   for (int i = 0 ; i < NumRows_ ; ++i) {
00120 
00121     int Nnz;
00122     ExtractMyRowCopy(i,MaxNumEntries_,Nnz,
00123          &Values[0], &Indices[0]);
00124     if (!TransA) {
00125       // no transpose first
00126       for (int j = 0 ; j < NumVectors ; ++j) {
00127   for (int k = 0 ; k < Nnz ; ++k) {
00128     Y[j][i] += Values[k] * X[j][Indices[k]];
00129   }
00130       }
00131     }
00132     else {
00133       // transpose here
00134       for (int j = 0 ; j < NumVectors ; ++j) {
00135   for (int k = 0 ; k < Nnz ; ++k) {
00136     Y[j][Indices[k]] += Values[k] * X[j][i];
00137   }
00138       }
00139     }
00140   }
00141 
00142   return(0);
00143 }
00144 
00145 //==============================================================================
00146 int Ifpack_DropFilter::
00147 Solve(bool Upper, bool Trans, bool UnitDiagonal, 
00148       const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00149 {
00150   IFPACK_CHK_ERR(-99);
00151 }
00152 
00153 //==============================================================================
00154 int Ifpack_DropFilter::
00155 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00156 {
00157   IFPACK_RETURN(Multiply(UseTranspose(),X,Y));
00158 }
00159 
00160 //==============================================================================
00161 int Ifpack_DropFilter::
00162 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00163 {
00164   IFPACK_CHK_ERR(-99);
00165 }
00166 
00167 //==============================================================================
00168 int Ifpack_DropFilter::InvRowSums(Epetra_Vector& x) const
00169 {
00170   IFPACK_CHK_ERR(-1);
00171 }
00172 
00173 //==============================================================================
00174 int Ifpack_DropFilter::InvColSums(Epetra_Vector& x) const
00175 {
00176   IFPACK_CHK_ERR(-1);
00177 }

Generated on Thu Sep 18 12:37:21 2008 for Ifpack Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1