Ifpack_SparseContainer.h

Go to the documentation of this file.
00001 #ifndef IFPACK_SPARSECONTAINER_H
00002 #define IFPACK_SPARSECONTAINER_H
00003 
00004 #include "Ifpack_Container.h"
00005 #ifdef HAVE_IFPACK_TEUCHOS
00006 #include "Epetra_IntSerialDenseVector.h"
00007 #include "Epetra_MultiVector.h"
00008 #include "Epetra_Vector.h"
00009 #include "Epetra_Map.h"
00010 #include "Epetra_RowMatrix.h"
00011 #include "Epetra_CrsMatrix.h"
00012 #include "Epetra_LinearProblem.h"
00013 #include "Epetra_IntSerialDenseVector.h"
00014 #include "Teuchos_ParameterList.hpp"
00015 #ifdef HAVE_MPI
00016 #include "Epetra_MpiComm.h"
00017 #else
00018 #include "Epetra_SerialComm.h"
00019 #endif
00020 
00050 template<typename T>
00051 class Ifpack_SparseContainer : public Ifpack_Container {
00052 
00053 public:
00054 
00056 
00057   Ifpack_SparseContainer(const int NumRows, const int NumVectors = 1);
00058 
00060   Ifpack_SparseContainer(const Ifpack_SparseContainer<T>& rhs);
00061 
00063   virtual ~Ifpack_SparseContainer();
00065 
00067 
00069   Ifpack_SparseContainer& operator=(const Ifpack_SparseContainer<T>& rhs);
00071 
00073 
00074   virtual int NumRows() const;
00075 
00077   virtual int NumVectors() const
00078   {
00079     return(NumVectors_);
00080   }
00081 
00083   virtual int SetNumVectors(const int NumVectors)
00084   {
00085     if (NumVectors_ == NumVectors)
00086       return(0);
00087     IFPACK_CHK_ERR(-99); // STILL TO DO
00088   }
00089 
00091   virtual double& LHS(const int i, const int Vector = 0);
00092   
00094   virtual double& RHS(const int i, const int Vector = 0);
00095 
00097 
00106   virtual int& ID(const int i);
00107 
00109   virtual int SetMatrixElement(const int row, const int col,
00110              const double value);
00111 
00112 
00114   virtual bool IsInitialized() const
00115   {
00116     return(IsInitialized_);
00117   }
00118 
00120   virtual bool IsComputed() const
00121   {
00122     return(IsComputed_);
00123   }
00124 
00126   virtual int SetParameters(Teuchos::ParameterList& List);
00127 
00129   virtual const char* Label() const
00130   {
00131     return(Label_.c_str());
00132   }
00133   
00135   const Epetra_Map* Map() const
00136   {
00137     return(Map_);
00138   }
00139 
00141   const Epetra_MultiVector* LHS() const
00142   {
00143     return(LHS_);
00144   }
00145 
00147   const Epetra_MultiVector* RHS() const
00148   {
00149     return(RHS_);
00150   }
00151 
00153   const Epetra_CrsMatrix* Matrix() const
00154   {
00155     return(Matrix_);
00156   }
00157 
00159   const Epetra_IntSerialDenseVector* ID() const
00160   {
00161     return(GID_);
00162   }
00163 
00165   const T* Inverse() const
00166   {
00167     return(Inverse_);
00168   }
00170 
00172 
00179   virtual int Initialize();
00181   virtual int Compute(const Epetra_RowMatrix& Matrix);
00183   virtual int Apply();
00184 
00186   virtual int ApplyInverse();
00187 
00189 
00191 
00192   virtual int Destroy();
00194 
00196   virtual double InitializeFlops() const
00197   {
00198     if (Inverse_ == 0)
00199       return (0.0);
00200     else
00201       return(Inverse_->InitializeFlops());
00202   }
00203 
00205   virtual double ComputeFlops() const
00206   {
00207     if (Inverse_ == 0)
00208       return (0.0);
00209     else
00210       return(Inverse_->ComputeFlops());
00211   }
00212 
00214   virtual double ApplyFlops() const
00215   {
00216     return(ApplyFlops_);
00217   }
00218 
00220   virtual double ApplyInverseFlops() const
00221   {
00222     if (Inverse_ == 0)
00223       return (0.0);
00224     else
00225       return(Inverse_->ApplyInverseFlops());
00226   }
00227   
00229   virtual ostream& Print(std::ostream& os) const;
00230 
00231 private:
00232   
00234   virtual int Extract(const Epetra_RowMatrix& Matrix);
00235 
00237   int NumRows_; 
00239   int NumVectors_; 
00241   Epetra_Map* Map_;
00243   Epetra_CrsMatrix* Matrix_;
00245   Epetra_MultiVector* LHS_;
00247   Epetra_MultiVector* RHS_;
00249   Epetra_IntSerialDenseVector GID_;
00251   bool IsInitialized_;
00253   bool IsComputed_;
00255   Epetra_Comm* SerialComm_;
00257   T* Inverse_;
00259   string Label_;
00260   Teuchos::ParameterList List_;
00261   double ApplyFlops_;
00262 
00263 };
00264 
00265 //==============================================================================
00266 template<typename T>
00267 Ifpack_SparseContainer<T>::
00268 Ifpack_SparseContainer(const int NumRows, const int NumVectors) :
00269   NumRows_(NumRows),
00270   NumVectors_(NumVectors),
00271   Map_(0),
00272   Matrix_(0),
00273   LHS_(0),
00274   RHS_(0),
00275   IsInitialized_(false),
00276   IsComputed_(false),
00277   SerialComm_(0),
00278   Inverse_(0),
00279   ApplyFlops_(0.0)
00280 {
00281 
00282 #ifdef HAVE_MPI
00283   SerialComm_ = new Epetra_MpiComm(MPI_COMM_SELF);
00284 #else
00285   SerialComm_ = new Epetra_SerialComm;
00286 #endif
00287 
00288 }
00289 
00290 //==============================================================================
00291 template<typename T>
00292 Ifpack_SparseContainer<T>::
00293 Ifpack_SparseContainer(const Ifpack_SparseContainer<T>& rhs) :
00294   NumRows_(rhs.NumRows()),
00295   NumVectors_(rhs.NumVectors()),
00296   Map_(0),
00297   Matrix_(0),
00298   LHS_(0),
00299   RHS_(0),
00300   IsInitialized_(rhs.IsInitialized()),
00301   IsComputed_(rhs.IsComputed()),
00302   SerialComm_(0),
00303   Inverse_(0)
00304 {
00305 
00306 #ifdef HAVE_MPI
00307   SerialComm_ = new Epetra_MpiComm(MPI_COMM_SELF);
00308 #else
00309   SerialComm_ = new Epetra_SerialComm;
00310 #endif
00311 
00312   if (rhs.Map())
00313     Map_ = new Epetra_Map(*rhs.Map());
00314 
00315   if (rhs.Matrix())
00316     Matrix_ = new Epetra_CrsMatrix(*rhs.Matrix());
00317 
00318   if (rhs.LHS())
00319     LHS_ = new Epetra_MultiVector(*rhs.LHS());
00320 
00321   if (rhs.RHS())
00322     RHS_ = new Epetra_MultiVector(*rhs.RHS());
00323 
00324 }
00325 //==============================================================================
00326 template<typename T>
00327 Ifpack_SparseContainer<T>::~Ifpack_SparseContainer()
00328 {
00329   Destroy();
00330 }
00331 
00332 //==============================================================================
00333 template<typename T>
00334 int Ifpack_SparseContainer<T>::NumRows() const
00335 {
00336   if (IsInitialized() == false)
00337     return(0);
00338   else
00339     return(NumRows_);
00340 }
00341 
00342 //==============================================================================
00343 template<typename T>
00344 int Ifpack_SparseContainer<T>::Initialize()
00345 {
00346   
00347   if (IsInitialized_ == true)
00348     Destroy();
00349 
00350   IsInitialized_ = false;
00351 
00352   Map_ = new Epetra_Map(NumRows_,0,*SerialComm_);
00353 
00354   LHS_ = new Epetra_MultiVector(*Map_,NumVectors_);
00355   RHS_ = new Epetra_MultiVector(*Map_,NumVectors_);
00356   GID_.Reshape(NumRows_,1);
00357 
00358   Matrix_ = new Epetra_CrsMatrix(Copy,*Map_,0);
00359 
00360   // create the inverse
00361   Inverse_ = new T(Matrix_);
00362 
00363   if (Inverse_ == 0)
00364     IFPACK_CHK_ERR(-5);
00365 
00366   IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
00367 
00368   // Call Inverse_->Initialize() in Compute(). This saves
00369   // some time, because I can extract the diagonal blocks faster,
00370   // and only once.
00371 
00372   Label_ = "Ifpack_SparseContainer";
00373 
00374   IsInitialized_ = true;
00375   return(0);
00376   
00377 }
00378 
00379 //==============================================================================
00380 template<typename T>
00381 double& Ifpack_SparseContainer<T>::LHS(const int i, const int Vector)
00382 {
00383   return(((*LHS_)(Vector))->Values()[i]);
00384 }
00385   
00386 //==============================================================================
00387 template<typename T>
00388 double& Ifpack_SparseContainer<T>::RHS(const int i, const int Vector)
00389 {
00390   return(((*RHS_)(Vector))->Values()[i]);
00391 }
00392 
00393 //==============================================================================
00394 template<typename T>
00395 int Ifpack_SparseContainer<T>::
00396 SetMatrixElement(const int row, const int col, const double value)
00397 {
00398   if (!IsInitialized())
00399     IFPACK_CHK_ERR(-3); // problem not shaped yet
00400 
00401   if ((row < 0) || (row >= NumRows())) {
00402     IFPACK_CHK_ERR(-2); // not in range
00403   }
00404 
00405   if ((col < 0) || (col >= NumRows())) {
00406     IFPACK_CHK_ERR(-2); // not in range
00407   }
00408 
00409   int ierr = Matrix_->InsertGlobalValues((int)row,1,(double*)&value,(int*)&col);
00410   if (ierr < 0) {
00411     ierr = Matrix_->SumIntoGlobalValues((int)row,1,(double*)&value,(int*)&col);
00412     if (ierr < 0)
00413       IFPACK_CHK_ERR(-1);
00414   }
00415 
00416   return(0);
00417 
00418 }
00419 
00420 //==============================================================================
00421 template<typename T>
00422 int Ifpack_SparseContainer<T>::Compute(const Epetra_RowMatrix& Matrix)
00423 {
00424 
00425   IsComputed_ = false;
00426   if (!IsInitialized()) {
00427     IFPACK_CHK_ERR(Initialize()); 
00428   }
00429 
00430   // extract the submatrices
00431   IFPACK_CHK_ERR(Extract(Matrix));
00432 
00433   // initialize the inverse operator
00434   IFPACK_CHK_ERR(Inverse_->Initialize());
00435 
00436   // compute the inverse operator
00437   IFPACK_CHK_ERR(Inverse_->Compute());
00438 
00439   Label_ = "Ifpack_SparseContainer";
00440   
00441   IsComputed_ = true;
00442 
00443   return(0);
00444 }
00445 
00446 //==============================================================================
00447 template<typename T>
00448 int Ifpack_SparseContainer<T>::Apply()
00449 {
00450   if (IsComputed() == false) {
00451     IFPACK_CHK_ERR(-3); // not yet computed
00452   }
00453   
00454   IFPACK_CHK_ERR(Matrix_->Apply(*RHS_, *LHS_));
00455 
00456   ApplyFlops_ += 2 * Matrix_->NumGlobalNonzeros();
00457   return(0);
00458 }
00459 
00460 //==============================================================================
00461 template<typename T>
00462 int Ifpack_SparseContainer<T>::ApplyInverse()
00463 {
00464   if (!IsComputed())
00465     IFPACK_CHK_ERR(-1);
00466   
00467   IFPACK_CHK_ERR(Inverse_->ApplyInverse(*RHS_, *LHS_));
00468 
00469   return(0);
00470 }
00471  
00472 
00473 //==============================================================================
00474 template<typename T>
00475 int Ifpack_SparseContainer<T>::Destroy()
00476 {
00477   if (Map_)
00478     delete Map_;
00479 
00480   if (Matrix_)
00481     delete Matrix_;
00482 
00483   if (LHS_)
00484     delete LHS_;
00485 
00486   if (RHS_)
00487     delete RHS_;
00488 
00489   if (Inverse_)
00490     delete Inverse_;
00491 
00492   Map_ = 0;
00493   Matrix_ = 0;
00494   Inverse_ = 0;
00495   LHS_ = 0;
00496   RHS_ = 0;
00497 
00498   IsInitialized_ = false;
00499   IsComputed_ = false;
00500   return(0);
00501 }
00502 
00503 //==============================================================================
00504 template<typename T>
00505 int& Ifpack_SparseContainer<T>::ID(const int i)
00506 {
00507   return(GID_[i]);
00508 }
00509 
00510 //==============================================================================
00511 template<typename T>
00512 int Ifpack_SparseContainer<T>::
00513 SetParameters(Teuchos::ParameterList& List)
00514 {
00515   List_ = List;
00516   return(0);
00517 }
00518 
00519 //==============================================================================
00520 // FIXME: optimize performances of this guy...
00521 template<typename T>
00522 int Ifpack_SparseContainer<T>::Extract(const Epetra_RowMatrix& Matrix)
00523 {
00524 
00525   for (int j = 0 ; j < NumRows_ ; ++j) {
00526     // be sure that the user has set all the ID's
00527     if (ID(j) == -1)
00528       IFPACK_CHK_ERR(-1);
00529     // be sure that all are local indices
00530     if (ID(j) > Matrix.NumMyRows())
00531       IFPACK_CHK_ERR(-1);
00532   }
00533 
00534   int Length = Matrix.MaxNumEntries();
00535   vector<double> Values;
00536   Values.resize(Length);
00537   vector<int> Indices;
00538   Indices.resize(Length);
00539 
00540   for (int j = 0 ; j < NumRows_ ; ++j) {
00541 
00542     int LRID = ID(j);
00543 
00544     int NumEntries;
00545 
00546     int ierr = 
00547       Matrix.ExtractMyRowCopy(LRID, Length, NumEntries, 
00548              &Values[0], &Indices[0]);
00549     IFPACK_CHK_ERR(ierr);
00550 
00551     for (int k = 0 ; k < NumEntries ; ++k) {
00552 
00553       int LCID = Indices[k];
00554 
00555       // skip off-processor elements
00556       if (LCID >= Matrix.NumMyRows()) 
00557   continue;
00558 
00559       // for local column IDs, look for each ID in the list
00560       // of columns hosted by this object
00561       // FIXME: use STL
00562       int jj = -1;
00563       for (int kk = 0 ; kk < NumRows_ ; ++kk)
00564   if (ID(kk) == LCID)
00565     jj = kk;
00566 
00567       if (jj != -1)
00568   SetMatrixElement(j,jj,Values[k]);
00569 
00570     }
00571   }
00572 
00573   IFPACK_CHK_ERR(Matrix_->FillComplete());
00574 
00575   return(0);
00576 }
00577 
00578 //==============================================================================
00579 template<typename T>
00580 ostream& Ifpack_SparseContainer<T>::Print(ostream & os) const
00581 {
00582   os << "================================================================================" << endl;
00583   os << "Ifpack_SparseContainer" << endl;
00584   os << "Number of rows          = " << NumRows() << endl;
00585   os << "Number of vectors       = " << NumVectors() << endl;
00586   os << "IsInitialized()         = " << IsInitialized() << endl;
00587   os << "IsComputed()            = " << IsComputed() << endl;
00588   os << "Flops in Initialize()   = " << InitializeFlops() << endl; 
00589   os << "Flops in Compute()      = " << ComputeFlops() << endl; 
00590   os << "Flops in ApplyInverse() = " << ApplyInverseFlops() << endl; 
00591   os << "================================================================================" << endl;
00592   os << endl;
00593 
00594   return(os);
00595 }
00596 #endif // HAVE_IFPACK_TEUCHOS
00597 #endif // IFPACK_SPARSECONTAINER_H

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