• Main Page
  • Related Pages
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

MLAPI_DistributedMatrix.h

Go to the documentation of this file.
00001 #ifndef MLAPI_DISTRIBUTEDMATRIX_H
00002 #define MLAPI_DISTRIBUTEDMATRIX_H
00003 
00013 /* ******************************************************************** */
00014 /* See the file COPYRIGHT for a complete copyright notice, contact      */
00015 /* person and disclaimer.                                               */        
00016 /* ******************************************************************** */
00017 
00018 #include "ml_common.h"
00019 
00020 #include "MLAPI_Error.h"
00021 #include "MLAPI_Operator.h"
00022 #include "MLAPI_Space.h"
00023 #include "MLAPI_BaseLinearCombination.h"
00024 #include "Epetra_Map.h"
00025 #include "Epetra_FECrsMatrix.h"
00026 
00027 namespace MLAPI {
00028 
00029 class DistributedMatrix : public Epetra_RowMatrix, public Operator {
00030 
00031 public:
00032 
00033   DistributedMatrix(const Space& RowSpace, const Space& ColSpace)
00034   {
00035     FillCompleted_ = false;
00036 
00037     RowSpace_ = RowSpace;
00038     ColSpace_ = ColSpace;
00039 
00040     int NumMyRows = RowSpace_.GetNumMyElements();
00041     int NumMyCols = ColSpace_.GetNumMyElements();
00042     
00043     // FIXME: add MyGlobalElements()
00044     RangeMap_ = new Epetra_Map(-1, NumMyRows, 0, GetEpetra_Comm());
00045     DomainMap_ = new Epetra_Map(-1, NumMyCols, 0, GetEpetra_Comm());
00046 
00047     Matrix_ = new Epetra_FECrsMatrix(Copy, *RangeMap_, 0);
00048   }
00049 
00050   virtual int NumMyRowEntries(int MyRow, int& NumEntries) const
00051   {
00052     ML_RETURN(Matrix_->NumMyRowEntries(MyRow, NumEntries));
00053   }
00054 
00055   virtual int MaxNumEntries() const
00056   {
00057     return(Matrix_->MaxNumEntries());
00058   }
00059 
00060   virtual int ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00061                                double* Values, int* Indices) const
00062   {
00063     ML_RETURN(Matrix_->ExtractMyRowCopy(MyRow, Length, NumEntries,
00064                                         Values, Indices));
00065   }
00066 
00067   virtual int ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00068   {
00069     ML_RETURN(Matrix_->ExtractDiagonalCopy(Diagonal));
00070   }
00071 
00072   virtual int Multiply(bool TransA, const Epetra_MultiVector& X, 
00073                        Epetra_MultiVector& Y) const
00074   {
00075     ML_RETURN(Matrix_->Multiply(TransA, X, Y));
00076   }
00077 
00078   virtual int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, 
00079                     Epetra_MultiVector& Y) const
00080   {
00081     ML_RETURN(Matrix_->Solve(Upper, Trans, UnitDiagonal, X, Y));
00082   }
00083 
00084   virtual int InvRowSums(Epetra_Vector& x) const
00085   {
00086     ML_RETURN(Matrix_->InvRowSums(x));
00087   }
00088 
00089   virtual int LeftScale(const Epetra_Vector& x)
00090   {
00091     ML_RETURN(Matrix_->LeftScale(x));
00092   }
00093 
00094   virtual int InvColSums(Epetra_Vector& x) const
00095   {
00096     ML_RETURN(Matrix_->InvColSums(x));
00097   }
00098 
00099   virtual int RightScale(const Epetra_Vector& x)
00100   {
00101     ML_RETURN(Matrix_->RightScale(x));
00102   }
00103 
00104   virtual bool Filled() const
00105   {
00106     return(Matrix_->Filled());
00107   }
00108 
00109   virtual double NormInf() const 
00110   {
00111     return(Matrix_->NormInf());
00112   }
00113 
00114   virtual double NormOne() const
00115   {
00116     return(Matrix_->NormOne());
00117   }
00118 
00119   virtual int NumGlobalNonzeros() const
00120   {
00121     return(Matrix_->NumGlobalNonzeros());
00122   }
00123 
00124   virtual int NumGlobalRows() const
00125   {
00126     return(Matrix_->NumGlobalRows());
00127   }
00128 
00129   virtual int NumGlobalCols() const
00130   {
00131     return(Matrix_->NumGlobalCols());
00132   }
00133 
00134   virtual int NumGlobalDiagonals() const
00135   {
00136     return(Matrix_->NumGlobalDiagonals());
00137   }
00138 
00139   virtual int NumMyNonzeros() const
00140   {
00141     return(Matrix_->NumMyNonzeros());
00142   }
00143 
00144   virtual int NumMyRows() const
00145   {
00146     return(Matrix_->NumMyRows());
00147   }
00148 
00149   virtual int NumMyCols() const
00150   {
00151     return(Matrix_->NumMyCols());
00152   }
00153 
00154   virtual int NumMyDiagonals() const
00155   {
00156     return(Matrix_->NumMyDiagonals());
00157   }
00158 
00159   virtual bool LowerTriangular() const
00160   {
00161     return(Matrix_->LowerTriangular());
00162   }
00163 
00164   virtual bool UpperTriangular() const 
00165   {
00166     return(Matrix_->UpperTriangular());
00167   }
00168 
00169   virtual const Epetra_Map & RowMatrixRowMap() const
00170   {
00171     return(Matrix_->RowMatrixRowMap());
00172   }
00173 
00174   virtual const Epetra_Map & RowMatrixColMap() const
00175   {
00176     return(Matrix_->RowMatrixColMap());
00177   }
00178 
00179   virtual const Epetra_Import * RowMatrixImporter() const
00180   {
00181     return(Matrix_->RowMatrixImporter());
00182   }
00183 
00184   virtual const Epetra_Map& OperatorDomainMap() const
00185   {
00186     return(Matrix_->OperatorDomainMap());
00187   }
00188 
00189   virtual const Epetra_Map& OperatorRangeMap() const
00190   {
00191     return(Matrix_->OperatorRangeMap());
00192   }
00193 
00194   virtual const Epetra_Map& Map() const
00195   {
00196     return(Matrix_->RowMatrixRowMap());
00197   }
00198     
00200 
00201   virtual int SetUseTranspose(bool what)
00202   {
00203     return(Matrix_->SetUseTranspose(what));
00204   }
00205   
00206   int Apply(const MultiVector& X, MultiVector& Y) const
00207   {
00208     return(Operator::Apply(X, Y));
00209   }
00210 
00211   virtual int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00212   {
00213     if (!IsFillCompleted())
00214       ML_THROW("Matrix not yet FillComplete()'d", -1);
00215 
00216     return(Matrix_->Apply(X, Y));
00217   }
00218 
00219   virtual int ApplyInverse(const Epetra_MultiVector& X,
00220                            Epetra_MultiVector& Y) const
00221   {
00222     if (!IsFillCompleted())
00223       ML_THROW("Matrix not yet FillComplete()'d", -1);
00224 
00225     return(Matrix_->ApplyInverse(X, Y));
00226   }
00227 
00228   virtual const char* Label() const
00229   {
00230     return(Matrix_->Label());
00231   }
00232 
00233   virtual bool UseTranspose() const
00234   {
00235     return(Matrix_->UseTranspose());
00236   }
00237 
00238   virtual bool HasNormInf() const
00239   {
00240     return(Matrix_->HasNormInf());
00241   }
00242 
00243   virtual const Epetra_Comm& Comm() const
00244   {
00245     return(Matrix_->Comm()); 
00246   }
00247 
00248   std::ostream& Print(std::ostream& os, const bool verbose = true) const;
00249 
00250   Space GetDomainSpace() const
00251   {
00252     return(ColSpace_);
00253   }
00254 
00255   Space GetRangeSpace() const
00256   {
00257     return(RowSpace_);
00258   }
00259 
00260   inline double& operator()(const int GRID, const int GCID)
00261   {
00262     if (IsFillCompleted())
00263     {
00264       ML_THROW("Matrix already FillCompleted()'d", -1);
00265     }
00266     else
00267     {
00268       rows_.push_back(GRID);
00269       cols_.push_back(GCID);
00270       vals_.push_back(0.0);
00271       return(vals_[vals_.size() - 1]);
00272     }
00273   }
00274   
00275   inline void ReplaceElement(const int GRID, const int GCID, 
00276                              const double value)
00277   {
00278     if (!IsFillCompleted())
00279       ML_THROW("Matrix not FillCompleted()'d yet", -1);
00280 
00281     int LRID = RangeMap_->LID(GRID);
00282     int LCID = Matrix_->ColMap().LID(GCID);
00283     if (Matrix_->ReplaceMyValues((int)LRID, 1, (double*)&value, 
00284                                  (int*)&LCID) < 0)
00285       ML_THROW("Can only replace locally owned elements", -1);
00286   }
00287 
00288   void FillComplete()
00289   {
00290     // populate the matrix here
00291     for (unsigned int i = 0 ; i < vals_.size() ; ++i)
00292     {
00293       int    GRID  = rows_[i];
00294       int    GCID  = cols_[i];
00295       double value = vals_[i];
00296       if (Matrix_->ReplaceGlobalValues(1, &GRID, 1, &GCID, &value) > 0)
00297         Matrix_->InsertGlobalValues(1, &GRID, 1, &GCID, &value);
00298     }
00299     rows_.resize(0);
00300     cols_.resize(0);
00301     vals_.resize(0);
00302 
00303     // freeze the matrix
00304     if (Matrix_->GlobalAssemble())
00305       ML_THROW("Error in GlobalAssemble()", -1);
00306 
00307     if (Matrix_->FillComplete(*DomainMap_, *RangeMap_))
00308       ML_THROW("Error in FillComplete()", -1);
00309 
00310     FillCompleted_ = true;
00311 
00312     Reshape(ColSpace_, RowSpace_, Matrix_, true);
00313   }
00314 
00315   bool IsFillCompleted() const 
00316   {
00317     return(FillCompleted_);
00318   }
00319 
00320 private:
00321 
00322   DistributedMatrix(const DistributedMatrix& rhs)
00323   {
00324   }
00325 
00326   DistributedMatrix& operator=(const DistributedMatrix& rhs)
00327   {
00328     return(*this);
00329   }
00330 
00331   mutable std::vector<int>    rows_;
00332   mutable std::vector<int>    cols_;
00333   mutable std::vector<double> vals_;
00334 
00335   Epetra_FECrsMatrix* Matrix_;
00336   Space ColSpace_;
00337   Space RowSpace_;
00338 
00339   Epetra_Map* DomainMap_;
00340   Epetra_Map* RangeMap_;
00341 
00342   bool FillCompleted_;
00343 
00344 }; // class DistributedMatrix
00345 
00346 } // namespace MLAPI
00347 
00348 #endif // ifndef MLAPI_DISTRIBUTEDMATRIX_H