Ifpack_SingletonFilter.cpp

Go to the documentation of this file.
00001 #include "Ifpack_ConfigDefs.h"
00002 #include "Ifpack_SingletonFilter.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 #include "Epetra_Import.h"
00010 
00011 //==============================================================================
00012 Ifpack_SingletonFilter::Ifpack_SingletonFilter(Epetra_RowMatrix* Matrix) :
00013   A_(*Matrix),
00014   NumSingletons_(0),
00015   NumRows_(0),
00016   NumRowsA_(0),
00017   MaxNumEntries_(0),
00018   MaxNumEntriesA_(0),
00019   NumNonzeros_(0),
00020   Map_(0),
00021   Diagonal_(0)
00022 {
00023   // use this filter only on serial matrices
00024   if (A_.Comm().NumProc() != 1) {
00025     cerr << "Ifpack_DropFilter can be used with Comm().NumProc() == 1" << endl;
00026     cerr << "only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
00027     cerr << "and it is not meant to be used otherwise." << endl;
00028     exit(EXIT_FAILURE);
00029   }
00030   
00031   if ((A_.NumMyRows() != A_.NumGlobalRows()) ||
00032      (A_.NumMyRows() != A_.NumMyCols()))
00033     IFPACK_CHK_ERRV(-1);
00034   
00035   NumRowsA_ = (A_.NumMyRows());
00036   MaxNumEntriesA_ = A_.MaxNumEntries();
00037 
00038   Indices_.resize(MaxNumEntriesA_);
00039   Values_.resize(MaxNumEntriesA_);
00040   Reorder_.resize(A_.NumMyRows());
00041 
00042   for (int i = 0 ; i < NumRowsA_ ; ++i)
00043     Reorder_[i] = -1;
00044 
00045   // first check how may singletons I do have
00046   for (int i = 0 ; i < NumRowsA_ ; ++i) {
00047     int Nnz;
00048     IFPACK_CHK_ERRV(A_.ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00049           &Values_[0], &Indices_[0]));
00050     if (Nnz != 1) {
00051       Reorder_[i] = NumRows_++;
00052     }
00053     else {
00054       NumSingletons_++;
00055     }
00056   }
00057 
00058   InvReorder_.resize(NumRows_);
00059   for (int i = 0 ; i < NumRowsA_ ; ++i) {
00060     if (Reorder_[i] < 0)
00061       continue;
00062     InvReorder_[Reorder_[i]] = i;
00063   }
00064  
00065   NumEntries_.resize(NumRows_);
00066   SingletonIndex_.resize(NumSingletons_);
00067 
00068   // now compute the nonzeros per row
00069   int count = 0;
00070   for (int i = 0 ; i < A_.NumMyRows() ; ++i) {
00071 
00072     int Nnz;
00073     IFPACK_CHK_ERRV(A_.ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00074             &Values_[0], &Indices_[0]));
00075     int ii = Reorder_[i];
00076     if (ii >= 0) {
00077       assert (Nnz != 1);
00078 
00079       NumEntries_[ii] = Nnz;
00080       NumNonzeros_ += Nnz;
00081       if (Nnz > MaxNumEntries_)
00082   MaxNumEntries_ = Nnz;
00083     }
00084     else {
00085       SingletonIndex_[count] = i;
00086       count++;
00087     }
00088   }
00089 
00090   Map_ = new Epetra_Map(NumRows_,0,Comm());
00091 
00092   // and finish up with the diagonal entry
00093   Diagonal_ = new Epetra_Vector(*Map_);
00094 
00095   Epetra_Vector Diagonal(A_.Map());
00096   A_.ExtractDiagonalCopy(Diagonal);
00097   for (int i = 0 ; i < NumRows_ ; ++i) {
00098     int ii = InvReorder_[i];
00099     (*Diagonal_)[i] = Diagonal[ii];
00100   }
00101   
00102 }
00103 
00104 //==============================================================================
00105 Ifpack_SingletonFilter::~Ifpack_SingletonFilter()
00106 {
00107   if (Diagonal_)
00108     delete Diagonal_;
00109 
00110   if (Map_)
00111     delete Map_;
00112 }
00113 
00114 //==============================================================================
00115 int Ifpack_SingletonFilter::
00116 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00117      double *Values, int * Indices) const
00118 {
00119   int Nnz;
00120 
00121   if (Length < NumEntries_[MyRow])
00122     IFPACK_CHK_ERR(-1);
00123 
00124   int Row = InvReorder_[MyRow];
00125   IFPACK_CHK_ERR(A_.ExtractMyRowCopy(Row,MaxNumEntriesA_,Nnz,
00126              &Values_[0],&Indices_[0]));
00127   NumEntries = 0;
00128   for (int i = 0 ; i < Nnz ; ++i) {
00129     int ii = Reorder_[Indices_[i]];
00130     if ( ii >= 0) {
00131       Indices[NumEntries] = ii;
00132       Values[NumEntries] = Values_[i];
00133       NumEntries++;
00134     }
00135   }
00136   return(0);   
00137 }
00138 
00139 //==============================================================================
00140 int Ifpack_SingletonFilter::
00141 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00142 {
00143   Diagonal = *Diagonal_;
00144   return(0);
00145 }
00146 
00147 //==============================================================================
00148 int Ifpack_SingletonFilter::
00149 Multiply(bool TransA, const Epetra_MultiVector& X, 
00150    Epetra_MultiVector& Y) const
00151 {
00152 
00153   int NumVectors = X.NumVectors();
00154   if (NumVectors != Y.NumVectors())
00155     IFPACK_CHK_ERR(-1);
00156 
00157   Y.PutScalar(0.0);
00158 
00159   vector<int> Indices(MaxNumEntries_);
00160   vector<double> Values(MaxNumEntries_);
00161 
00162   // cycle over all rows of the original matrix
00163   for (int i = 0 ; i < A_.NumMyRows() ; ++i) {
00164 
00165     if (Reorder_[i] < 0)
00166       continue; // skip singleton rows
00167     
00168     int Nnz;
00169     A_.ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00170          &Values[0], &Indices[0]);
00171     if (!TransA) {
00172       // no transpose first
00173       for (int j = 0 ; j < NumVectors ; ++j) {
00174   for (int k = 0 ; k < Nnz ; ++k) {
00175     if (Reorder_[i] >= 0)
00176       Y[j][i] += Values[k] * X[j][Reorder_[Indices[k]]];
00177   }
00178       }
00179     }
00180     else {
00181       // transpose here
00182       for (int j = 0 ; j < NumVectors ; ++j) {
00183   for (int k = 0 ; k < Nnz ; ++k) {
00184     if (Reorder_[i] >= 0)
00185       Y[j][Reorder_[Indices[k]]] += Values[k] * X[j][i];
00186   }
00187       }
00188     }
00189   }
00190 
00191   return(0);
00192 }
00193 
00194 //==============================================================================
00195 int Ifpack_SingletonFilter::
00196 Solve(bool Upper, bool Trans, bool UnitDiagonal, 
00197       const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00198 {
00199   return(-1);
00200 }
00201 
00202 //==============================================================================
00203 int Ifpack_SingletonFilter::
00204 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00205 {
00206   IFPACK_CHK_ERR(Multiply(false,X,Y));
00207   return(0);
00208 }
00209 
00210 //==============================================================================
00211 int Ifpack_SingletonFilter::
00212 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00213 {
00214   return(-1); // NOT IMPLEMENTED AT THIS STAGE
00215 }
00216 
00217 //==============================================================================
00218 int Ifpack_SingletonFilter::
00219 SolveSingletons(const Epetra_MultiVector& RHS, 
00220     Epetra_MultiVector& LHS)
00221 {
00222   for (int i = 0 ; i < NumSingletons_ ; ++i) {
00223     int ii = SingletonIndex_[i];
00224     // get the diagonal value for the singleton
00225     int Nnz;
00226     A_.ExtractMyRowCopy(ii,MaxNumEntriesA_,Nnz,
00227       &Values_[0], &Indices_[0]);
00228     for (int j = 0 ; j < Nnz ; ++j) {
00229       if (Indices_[j] == ii) {
00230   for (int k = 0 ; k < LHS.NumVectors() ; ++k)
00231     LHS[k][ii] = RHS[k][ii] / Values_[j];
00232       }
00233     }
00234   }
00235 
00236   return(0);
00237 }
00238 
00239 //==============================================================================
00240 int Ifpack_SingletonFilter::
00241 CreateReducedRHS(const Epetra_MultiVector& LHS,
00242      const Epetra_MultiVector& RHS, 
00243      Epetra_MultiVector& ReducedRHS)
00244 {
00245   int NumVectors = LHS.NumVectors();
00246 
00247   for (int i = 0 ; i < NumRows_ ; ++i)
00248     for (int k = 0 ; k < NumVectors ; ++k)
00249       ReducedRHS[k][i] = RHS[k][InvReorder_[i]];
00250 
00251   for (int i = 0 ; i < NumRows_ ; ++i) {
00252     int ii = InvReorder_[i];
00253     int Nnz;
00254     IFPACK_CHK_ERR(A_.ExtractMyRowCopy(ii,MaxNumEntriesA_,Nnz,
00255           &Values_[0], &Indices_[0]));
00256 
00257     for (int j = 0 ; j < Nnz ; ++j) {
00258       if (Reorder_[Indices_[j]] == -1) {
00259   for (int k = 0 ; k < NumVectors ; ++k)
00260     ReducedRHS[k][i] -= Values_[j] * LHS[k][Indices_[j]];
00261       }
00262     }
00263   }
00264   return(0);
00265 }
00266 
00267 //==============================================================================
00268 int Ifpack_SingletonFilter::
00269 UpdateLHS(const Epetra_MultiVector& ReducedLHS,
00270     Epetra_MultiVector& LHS)
00271 {
00272   for (int i = 0 ; i < NumRows_ ; ++i)
00273     for (int k = 0 ; k < LHS.NumVectors() ; ++k)
00274       LHS[k][InvReorder_[i]] = ReducedLHS[k][i];
00275 
00276   return(0);
00277 }

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