Ifpack_LocalFilter.cpp

Go to the documentation of this file.
00001 #include "Ifpack_ConfigDefs.h"
00002 
00003 #include "Epetra_MultiVector.h"
00004 #include "Epetra_Vector.h"
00005 #include "Epetra_RowMatrix.h"
00006 #include "Epetra_Map.h"
00007 #include "Epetra_BlockMap.h"
00008 #include "Ifpack_LocalFilter.h"
00009 
00010 //==============================================================================
00011 Ifpack_LocalFilter::Ifpack_LocalFilter(const Epetra_RowMatrix* Matrix) :
00012   Matrix_(Matrix),
00013   SerialComm_(0),
00014   Map_(0),
00015   NumRows_(0),
00016   NumNonzeros_(0),
00017   MaxNumEntries_(0),
00018   MaxNumEntriesA_(0),
00019   Diagonal_(0)
00020 {
00021   sprintf(Label_,"%s","Ifpack_LocalFilter");
00022 
00023 #ifdef HAVE_MPI
00024   SerialComm_ = new Epetra_MpiComm(MPI_COMM_SELF);
00025 #else
00026   SerialComm_ = new Epetra_SerialComm;
00027 #endif
00028 
00029   // localized matrix has all the local rows of Matrix
00030   NumRows_ = Matrix->NumMyRows();
00031 
00032   // build a linear map, based on the serial communicator
00033   Map_ = new Epetra_Map(NumRows_,0,*SerialComm_);
00034 
00035   // NumEntries_ will contain the actual number of nonzeros
00036   // for each localized row (that is, without external nodes,
00037   // and always with the diagonal entry)
00038   NumEntries_.resize(NumRows_);
00039 
00040   // want to store the diagonal vector. FIXME: am I really useful?
00041   Diagonal_ = new Epetra_Vector(*Map_);
00042   if (Diagonal_ == 0) IFPACK_CHK_ERRV(-5);
00043   
00044   // store this for future access to ExtractMyRowCopy().
00045   // This is the # of nonzeros in the non-local matrix
00046   MaxNumEntriesA_ = Matrix->MaxNumEntries();
00047   // tentative value for MaxNumEntries. This is the number of
00048   // nonzeros in the local matrix
00049   MaxNumEntries_ = Matrix->MaxNumEntries();
00050 
00051   // ExtractMyRowCopy() will use these vectors
00052   Indices_.resize(MaxNumEntries_);
00053   Values_.resize(MaxNumEntries_);
00054 
00055   // now compute:
00056   // - the number of nonzero per row
00057   // - the total number of nonzeros
00058   // - the diagonal entries
00059 
00060   // ExtractMyRowCopy() uses Indices_ and Values_ for
00061   // its internal operators, here I need two other vectors
00062   vector<int> Ind(MaxNumEntries_);
00063   vector<double> Val(MaxNumEntries_);
00064 
00065   // compute nonzeros (total and per-row), and store the
00066   // diagonal entries (already modified)
00067   int ActualMaxNumEntries = 0;
00068 
00069   for (int i = 0 ; i < NumRows_ ; ++i) {
00070     
00071     NumEntries_[i] = 0;
00072     int Nnz;
00073     IFPACK_CHK_ERRV(ExtractMyRowCopy(i,MaxNumEntries_,Nnz,&Val[0],&Ind[0]));
00074 
00075     if (Nnz > ActualMaxNumEntries)
00076       ActualMaxNumEntries = Nnz;
00077 
00078     NumNonzeros_ += Nnz;
00079     NumEntries_[i] = Nnz;
00080 
00081     for (int j = 0 ; j < Nnz ; ++j) {
00082       if (Ind[j] == i)
00083   (*Diagonal_)[i] = Val[j];
00084     }
00085   }
00086  
00087   MaxNumEntries_ = ActualMaxNumEntries;
00088 }
00089 
00090 //==============================================================================
00091 Ifpack_LocalFilter::~Ifpack_LocalFilter()
00092 {
00093   if (Diagonal_)
00094     delete Diagonal_;
00095   if (Map_)
00096     delete Map_;
00097   if (SerialComm_)
00098     delete SerialComm_;
00099 }
00100 
00101 //==============================================================================
00102 int Ifpack_LocalFilter::
00103 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00104      double *Values, int * Indices) const
00105 {
00106   if ((MyRow < 0) || (MyRow >= NumRows_)) {
00107     IFPACK_CHK_ERR(-1); // range not valid
00108   }
00109 
00110   if (Length < NumEntries_[MyRow])
00111     return(-1);
00112 
00113   // always extract using the object Values_ and Indices_.
00114   // This is because I need more space than that given by
00115   // the user (for the external nodes)
00116   int Nnz;
00117   int ierr = Matrix_->ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz,
00118                &Values_[0],&Indices_[0]);
00119 
00120   IFPACK_CHK_ERR(ierr);
00121 
00122   // populate the user's vectors, add diagonal if not found
00123   NumEntries = 0;
00124 
00125   for (int j = 0 ; j < Nnz ; ++j) {
00126     // only local indices
00127     if (Indices_[j] < NumRows_ ) {
00128       Indices[NumEntries] = Indices_[j];
00129       Values[NumEntries] = Values_[j];
00130       ++NumEntries;
00131     }
00132   }
00133     
00134   return(0);
00135 
00136 }
00137 
00138 //==============================================================================
00139 int Ifpack_LocalFilter::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00140 {
00141   if (!Diagonal.Map().SameAs(*Map_))
00142     IFPACK_CHK_ERR(-1);
00143   Diagonal = *Diagonal_;
00144   return(0);
00145 }
00146 
00147 //==============================================================================
00148 int Ifpack_LocalFilter::Apply(const Epetra_MultiVector& X,
00149     Epetra_MultiVector& Y) const 
00150 {
00151 
00152   // skip expensive checks, I suppose input data are ok
00153 
00154   Y.PutScalar(0.0);
00155   int NumVectors = Y.NumVectors();
00156 
00157   double** X_ptr;
00158   double** Y_ptr;
00159   X.ExtractView(&X_ptr);
00160   Y.ExtractView(&Y_ptr);
00161 
00162   for (int i = 0 ; i < NumRows_ ; ++i) {
00163     
00164     int Nnz;
00165     int ierr = Matrix_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,&Values_[0],
00166                                          &Indices_[0]);
00167     IFPACK_CHK_ERR(ierr);
00168 
00169     for (int j = 0 ; j < Nnz ; ++j) {
00170       if (Indices_[j] < NumRows_ ) {
00171   for (int k = 0 ; k < NumVectors ; ++k)
00172     Y_ptr[k][i] += Values_[j] * X_ptr[k][Indices_[j]];
00173       }
00174     }
00175   }
00176 
00177   return(0);
00178 }
00179 
00180 //==============================================================================
00181 int Ifpack_LocalFilter::ApplyInverse(const Epetra_MultiVector& X,
00182      Epetra_MultiVector& Y) const
00183 {
00184   IFPACK_CHK_ERR(-1); // not implemented
00185 }
00186 
00187 //==============================================================================
00188 const Epetra_BlockMap& Ifpack_LocalFilter::Map() const
00189 {
00190   return(*Map_);
00191 }

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