00001 #ifndef ML_OPERATOR_H
00002 #define ML_OPERATOR_H
00003
00013
00014
00015
00016
00017 #include "ml_common.h"
00018 #include <iostream>
00019 #include "ml_operator.h"
00020 #include "ml_epetra_utils.h"
00021 #include "ml_RowMatrix.h"
00022 #include "Teuchos_RefCountPtr.hpp"
00023 #include "MLAPI_Error.h"
00024 #include "MLAPI_Space.h"
00025 #include "MLAPI_MultiVector.h"
00026 #include "MLAPI_BaseOperator.h"
00027 #include "MLAPI_CompObject.h"
00028 #include "MLAPI_TimeObject.h"
00029 #include "MLAPI_Workspace.h"
00030 #include "MLAPI_Operator_Box.h"
00031
00032 namespace MLAPI {
00033
00044 class Operator : public BaseOperator, public CompObject, public TimeObject {
00045
00046 public:
00048
00050 Operator()
00051 {
00052 RCPOperatorBox_ = Teuchos::null;
00053 }
00054
00056 Operator(const Space& DomainSpace, const Space& RangeSpace,
00057 ML_Operator* Op, bool Ownership = true,
00058 Teuchos::RefCountPtr<ML_Operator_Box> AuxOp = Teuchos::null)
00059 {
00060 Reshape(DomainSpace, RangeSpace, Op, Ownership, AuxOp);
00061 }
00062
00064 Operator(const Space& DomainSpace, const Space& RangeSpace,
00065 Epetra_RowMatrix* Matrix, bool Ownership = true,
00066 Teuchos::RefCountPtr<ML_Operator_Box> AuxOp = Teuchos::null)
00067 {
00068 Reshape(DomainSpace, RangeSpace, Matrix, Ownership, AuxOp);
00069 }
00070
00072 Operator(const Operator& RHS)
00073 {
00074 DomainSpace_ = RHS.GetDomainSpace();
00075 RangeSpace_ = RHS.GetRangeSpace();
00076 ColumnSpace_ = RHS.GetColumnSpace();
00077 RCPOperatorBox_ = RHS.GetRCPOperatorBox();
00078 RCPAuxOperatorBox_ = RHS.GetRCPAuxOperatorBox();
00079 RCPRowMatrix_ = RHS.GetRCPRowMatrix();
00080
00081 SetLabel(RHS.GetLabel());
00082 }
00083
00085 ~Operator()
00086 {
00087 Destroy();
00088 }
00089
00090
00091
00092
00094 void Reshape()
00095 {
00096 Destroy();
00097 }
00098
00100 void Reshape(const Space& DomainSpace, const Space& RangeSpace,
00101 ML_Operator* Op, bool Ownership = true,
00102 Teuchos::RefCountPtr<ML_Operator_Box> AuxOp = Teuchos::null)
00103 {
00104 StackPush();
00105
00106 RangeSpace_ = RangeSpace;
00107 DomainSpace_ = DomainSpace;
00108
00109 RCPOperatorBox_ = Teuchos::rcp(new ML_Operator_Box(Op,Ownership));
00110
00111 bool cheap;
00112 if (RangeSpace_ != DomainSpace_)
00113 cheap = true;
00114 else
00115 cheap = false;
00116 RCPRowMatrix_ = Teuchos::rcp(new ML_Epetra::RowMatrix(Op,&(GetEpetra_Comm()),
00117 cheap));
00118 RCPAuxOperatorBox_ = AuxOp;
00119
00120 StackPop();
00121 }
00122
00124 void Reshape(const Space& DomainSpace, const Space& RangeSpace,
00125 Epetra_RowMatrix* Matrix, bool Ownership = true,
00126 Teuchos::RefCountPtr<ML_Operator_Box> AuxOp = Teuchos::null)
00127 {
00128 StackPush();
00129
00130 RangeSpace_ = RangeSpace;
00131 DomainSpace_ = DomainSpace;
00132
00133 ML_Operator* Op = ML_Operator_Create(MLAPI::GetML_Comm());
00134 RCPOperatorBox_ = Teuchos::rcp(new ML_Operator_Box(Op,true));
00135 RCPAuxOperatorBox_ = AuxOp;
00136
00137 RCPRowMatrix_ = Teuchos::rcp(Matrix,Ownership);
00138 ML_Operator_WrapEpetraMatrix(RCPRowMatrix_.get(), GetML_Operator());
00139
00140 StackPop();
00141 }
00142
00143
00144
00145
00147 Operator& operator=(const Operator& RHS)
00148 {
00149 StackPush();
00150
00151 Destroy();
00152
00153 DomainSpace_ = RHS.GetDomainSpace();
00154 RangeSpace_ = RHS.GetRangeSpace();
00155 ColumnSpace_ = RHS.GetColumnSpace();
00156 RCPOperatorBox_ = RHS.GetRCPOperatorBox();
00157 RCPRowMatrix_ = RHS.GetRCPRowMatrix();
00158
00159 SetLabel(RHS.GetLabel());
00160
00161 StackPop();
00162
00163 return(*this);
00164 }
00165
00167 inline Operator& operator=(const string& Label)
00168 {
00169 SetLabel(Label);
00170 return(*this);
00171 }
00172
00173
00174
00175
00177 const Space GetOperatorDomainSpace() const {
00178 return(DomainSpace_);
00179 }
00180
00182 const Space GetOperatorRangeSpace() const {
00183 return(RangeSpace_);
00184 }
00185
00187 inline const Space GetDomainSpace() const {
00188 return(DomainSpace_);
00189 }
00190
00192 inline const Space GetRangeSpace() const {
00193 return(RangeSpace_);
00194 }
00195
00197 inline const Space GetColumnSpace() const
00198 {
00199 return(ColumnSpace_);
00200 }
00201
00203 inline int GetNumGlobalRows() const
00204 {
00205 return(GetRangeSpace().GetNumGlobalElements());
00206 }
00207
00209 inline int GetNumMyRows() const
00210 {
00211 return(GetRangeSpace().GetNumMyElements());
00212 }
00213
00215 inline int GetNumGlobalCols() const
00216 {
00217 return(GetRowMatrix()->NumGlobalCols());
00218 }
00219
00221 inline int GetNumMyCols() const
00222 {
00223 return(GetRowMatrix()->NumMyCols());
00224 }
00225
00227 inline int GetNumGlobalNonzeros() const
00228 {
00229 return(GetRowMatrix()->NumGlobalNonzeros());
00230 }
00231
00233 inline int GetNumMyNonzeros() const
00234 {
00235 return(GetRowMatrix()->NumMyNonzeros());
00236 }
00237
00239 inline const Epetra_RowMatrix* GetRowMatrix() const
00240 {
00241 return(RCPRowMatrix_.get());
00242 }
00243
00245 inline ML_Operator* GetML_Operator() const
00246 {
00247 return(GetRCPOperatorBox()->GetData());
00248 }
00249
00251 inline const Teuchos::RefCountPtr<ML_Operator_Box>& GetRCPOperatorBox() const
00252 {
00253 return(RCPOperatorBox_);
00254 }
00255
00257 inline const Teuchos::RefCountPtr<ML_Operator_Box>& GetRCPAuxOperatorBox() const
00258 {
00259 return(RCPAuxOperatorBox_);
00260 }
00262 inline const Teuchos::RefCountPtr<Epetra_RowMatrix>& GetRCPRowMatrix() const
00263 {
00264 return(RCPRowMatrix_);
00265 }
00266
00268 int GetGRID(const int LRID) const
00269 {
00270 #ifdef MLAPI_CHECK
00271 if (LRID < 0 || LRID >= GetNumMyRows())
00272 ML_THROW("LRID in invalid", -1);
00273 #endif
00274 return(GetRangeSpace()(LRID));
00275 }
00276
00278 int GetGCID(const int LCID) const
00279 {
00280 #ifdef MLAPI_CHECK
00281 if (LCID < 0 || LCID >= GetRowMatrix()->NumMyCols())
00282 ML_THROW("LRID in invalid", -1);
00283 #endif
00284 return(GetRowMatrix()->RowMatrixColMap().GID(LCID));
00285 }
00286
00287
00288
00289
00291 int Apply(const MultiVector& X, MultiVector& Y) const
00292 {
00293 ResetTimer();
00294 StackPush();
00295
00296 if (GetDomainSpace() != X.GetVectorSpace())
00297 ML_THROW("Domain spaces differ", -1);
00298 if (GetRangeSpace() != Y.GetVectorSpace())
00299 ML_THROW("Range spaces differ", -1);
00300 if (X.GetNumVectors() != Y.GetNumVectors())
00301 ML_THROW("Number of vectors differ", -1);
00302 if (GetML_Operator() == 0)
00303 ML_THROW("Operator not set", -1);
00304
00305 int (*func)(ML_Operator*,int,double*,int,double*) =
00306 GetML_Operator()->matvec->func_ptr;
00307
00308 for (int v = 0 ; v < X.GetNumVectors() ; ++v) {
00309 double* x_ptr = (double*)&X(0) + v * X.GetMyLength();
00310 double* y_ptr = (double*)&Y(0) + v * Y.GetMyLength();
00311 (*func)(GetML_Operator(),X.GetMyLength(),x_ptr,
00312 Y.GetMyLength(), y_ptr);
00313 }
00314
00315 StackPop();
00316
00317 UpdateFlops(2.0 * GetNumGlobalNonzeros());
00318 UpdateTime();
00319
00320 return(0);
00321 }
00322
00323
00324
00325
00327 ostream& Print(std::ostream& os, const bool verbose = true) const
00328 {
00329 if (GetRCPOperatorBox().get() == 0) {
00330 if (GetMyPID() == 0) {
00331 os << endl;
00332 os << "*** MLAPI::Operator ***" << endl;
00333 os << "Label = " << GetLabel() << endl;
00334 os << "Status = empty" << endl;
00335 os << endl;
00336 }
00337 return(os);
00338 }
00339
00340 StackPush();
00341
00342 int *bindx;
00343 double *val;
00344 int allocated, row_length;
00345 ML_Operator* matrix = GetML_Operator();
00346
00347 if (matrix->getrow == NULL)
00348 ML_THROW("getrow not set", -1);
00349
00350 if (GetMyPID() == 0) {
00351 os << endl;
00352 os << "*** MLAPI::Operator ***" << endl;
00353 os << "Label = " << GetLabel() << endl;
00354 os << "Number of rows = " << GetRangeSpace().GetNumGlobalElements() << endl;
00355 os << "Number of columns = " << GetDomainSpace().GetNumGlobalElements() << endl;
00356 os << "Flop count = " << GetFlops() << endl;
00357 os << "Cumulative time = " << GetTime() << endl;
00358 if (GetTime() != 0.0)
00359 os << "MFlops rate = " << 1.0e-6 * GetFlops() / GetTime() << endl;
00360 else
00361 os << "MFlops rate = 0.0" << endl;
00362 os << endl;
00363 }
00364
00365 if (!verbose)
00366 return(os);
00367
00368 allocated = 100;
00369 bindx = (int *) ML_allocate(allocated*sizeof(int ));
00370 val = (double *) ML_allocate(allocated*sizeof(double));
00371
00372 if (GetMyPID() == 0) {
00373 os.width(10);
00374 os << "ProcID";
00375 os.width(20);
00376 os << "Global Row";
00377 os.width(20);
00378 os << "Global Col";
00379 os.width(20);
00380 os << "Value" << endl;
00381 os << endl;
00382 }
00383
00384 for (int iproc = 0 ; iproc < GetNumProcs() ; ++iproc) {
00385
00386 if (GetMyPID() == iproc) {
00387
00388 for (int i = 0 ; i < matrix->getrow->Nrows; i++) {
00389 ML_get_matrix_row(matrix, 1, &i, &allocated, &bindx, &val,
00390 &row_length, 0);
00391 for (int j = 0; j < row_length; j++) {
00392 int GlobalRow = GetRangeSpace()(i);
00393
00394 int GlobalCol = GetRowMatrix()->RowMatrixColMap().GID(bindx[j]);
00395 os.width(10);
00396 os << iproc;
00397 os.width(20);
00398 os << GlobalRow;
00399 os.width(20);
00400 os << GlobalCol;
00401 os.width(20);
00402 os << val[j] << endl;
00403 }
00404 }
00405 }
00406 Barrier();
00407 }
00408
00409 if (GetMyPID() == 0)
00410 os << endl;
00411
00412 Barrier();
00413
00414 ML_free(val);
00415 ML_free(bindx);
00416
00417 StackPop();
00418
00419 return (os);
00420 }
00421
00423 void BuildColumnSpace()
00424 {
00425 StackPush();
00426
00427 if (GetNumProcs() == 1) {
00428 ColumnSpace_ = DomainSpace_;
00429 return;
00430 }
00431
00432 std::vector<double> dtemp;
00433 std::vector<int> GlobalElements;
00434
00435 int Nrows = GetML_Operator()->getrow->Nrows;
00436 int Nghosts;
00437 if (GetML_Operator()->getrow->pre_comm == NULL) Nghosts = 0;
00438 else {
00439 if (GetML_Operator()->getrow->pre_comm->total_rcv_length <= 0)
00440 ML_CommInfoOP_Compute_TotalRcvLength(GetML_Operator()->getrow->pre_comm);
00441 Nghosts = GetML_Operator()->getrow->pre_comm->total_rcv_length;
00442 }
00443
00444 dtemp.resize(Nrows + Nghosts);
00445
00446 for (int i = 0 ; i < Nrows ; ++i)
00447 dtemp[i] = 1.0 * GetDomainSpace()(i);
00448 for (int i = 0 ; i < Nghosts; ++i)
00449 dtemp[i + Nrows] = -1;
00450
00451 ML_exchange_bdry(&dtemp[0],GetML_Operator()->getrow->pre_comm,
00452 GetML_Operator()->outvec_leng,
00453 GetML_Comm(), ML_OVERWRITE,NULL);
00454
00455 GlobalElements.resize(Nrows + Nghosts);
00456
00457 for (int i = 0 ; i < Nrows + Nghosts ; ++i)
00458 GlobalElements[i] = (int)dtemp[i];
00459
00460 ColumnSpace_.Reshape(-1, Nrows + Nghosts, &GlobalElements[0]);
00461
00462 StackPop();
00463
00464 return;
00465 }
00466
00467
00468
00469 private:
00470
00472 void Destroy()
00473 {
00474 RangeSpace_.Reshape();
00475 DomainSpace_.Reshape();
00476 RCPOperatorBox_ = Teuchos::null;
00477 RCPRowMatrix_ = Teuchos::null;
00478 RCPAuxOperatorBox_ = Teuchos::null;
00479 }
00480
00482 Space DomainSpace_;
00484 Space RangeSpace_;
00486 Space ColumnSpace_;
00488 Teuchos::RefCountPtr<ML_Operator_Box> RCPOperatorBox_;
00490 Teuchos::RefCountPtr<ML_Operator_Box> RCPAuxOperatorBox_;
00492 Teuchos::RefCountPtr<Epetra_RowMatrix> RCPRowMatrix_;
00493
00494 };
00495
00496 }
00497
00498 #endif // ML_OPERATOR_H