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

MLAPI_Operator.h

Go to the documentation of this file.
00001 #ifndef ML_OPERATOR_H
00002 #define ML_OPERATOR_H
00003 
00013 /* ******************************************************************** */
00014 /* See the file COPYRIGHT for a complete copyright notice, contact      */
00015 /* person and disclaimer.                                               */        
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   // @{ \name Reshape methods
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     // NOTE: Should the code crash, change the following to `false'
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   // @{ \name Overloaded operators
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   // @{ \name Get and Set methods
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   // @{ \name Mathematical methods.
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   // @{ \name Miscellaneous methods
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             //int GlobalCol = GetColumnSpace()(bindx[j]);
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 }; // Operator
00495 
00496 } // namespace MLAPI
00497 
00498 #endif // ML_OPERATOR_H