00001 #ifndef MLAPI_DISTRIBUTEDMATRIX_H
00002 #define MLAPI_DISTRIBUTEDMATRIX_H
00003
00013
00014
00015
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
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
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
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 };
00345
00346 }
00347
00348 #endif // ifndef MLAPI_DISTRIBUTEDMATRIX_H