Ifpack_SparseContainer.h

Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #ifndef IFPACK_SPARSECONTAINER_H
00031 #define IFPACK_SPARSECONTAINER_H
00032 
00033 #include "Ifpack_Container.h"
00034 #include "Epetra_IntSerialDenseVector.h"
00035 #include "Epetra_MultiVector.h"
00036 #include "Epetra_Vector.h"
00037 #include "Epetra_Map.h"
00038 #include "Epetra_RowMatrix.h"
00039 #include "Epetra_CrsMatrix.h"
00040 #include "Epetra_LinearProblem.h"
00041 #include "Epetra_IntSerialDenseVector.h"
00042 #include "Teuchos_ParameterList.hpp"
00043 #include "Teuchos_RefCountPtr.hpp"
00044 #ifdef HAVE_MPI
00045 #include "Epetra_MpiComm.h"
00046 #else
00047 #include "Epetra_SerialComm.h"
00048 #endif
00049 
00079 template<typename T>
00080 class Ifpack_SparseContainer : public Ifpack_Container {
00081 
00082 public:
00083 
00085 
00086   Ifpack_SparseContainer(const int NumRows, const int NumVectors = 1);
00087 
00089   Ifpack_SparseContainer(const Ifpack_SparseContainer<T>& rhs);
00090 
00092   virtual ~Ifpack_SparseContainer();
00094 
00096 
00098   Ifpack_SparseContainer& operator=(const Ifpack_SparseContainer<T>& rhs);
00100 
00102 
00103   virtual int NumRows() const;
00104 
00106   virtual int NumVectors() const
00107   {
00108     return(NumVectors_);
00109   }
00110 
00112   virtual int SetNumVectors(const int NumVectors_in)
00113   {
00114     if (NumVectors_ == NumVectors_in)
00115       return(0);
00116     IFPACK_CHK_ERR(-99); // STILL TO DO
00117   }
00118 
00120   virtual double& LHS(const int i, const int Vector = 0);
00121   
00123   virtual double& RHS(const int i, const int Vector = 0);
00124 
00126 
00135   virtual int& ID(const int i);
00136 
00138   virtual int SetMatrixElement(const int row, const int col,
00139              const double value);
00140 
00141 
00143   virtual bool IsInitialized() const
00144   {
00145     return(IsInitialized_);
00146   }
00147 
00149   virtual bool IsComputed() const
00150   {
00151     return(IsComputed_);
00152   }
00153 
00155   virtual int SetParameters(Teuchos::ParameterList& List);
00156 
00158   virtual const char* Label() const
00159   {
00160     return(Label_.c_str());
00161   }
00162   
00164   const Epetra_Map* Map() const
00165   {
00166     return(Map_);
00167   }
00168 
00170   const Epetra_MultiVector* LHS() const
00171   {
00172     return(LHS_);
00173   }
00174 
00176   const Epetra_MultiVector* RHS() const
00177   {
00178     return(RHS_);
00179   }
00180 
00182   const Epetra_CrsMatrix* Matrix() const
00183   {
00184     return(Matrix_);
00185   }
00186 
00188   const Epetra_IntSerialDenseVector* ID() const
00189   {
00190     return(GID_);
00191   }
00192 
00194   const T* Inverse() const
00195   {
00196     return(Inverse_);
00197   }
00199 
00201 
00208   virtual int Initialize();
00210   virtual int Compute(const Epetra_RowMatrix& Matrix_in);
00212   virtual int Apply();
00213 
00215   virtual int ApplyInverse();
00216 
00218 
00220 
00221   virtual int Destroy();
00223 
00225   virtual double InitializeFlops() const
00226   {
00227     if (Inverse_ == Teuchos::null)
00228       return (0.0);
00229     else
00230       return(Inverse_->InitializeFlops());
00231   }
00232 
00234   virtual double ComputeFlops() const
00235   {
00236     if (Inverse_ == Teuchos::null)
00237       return (0.0);
00238     else
00239       return(Inverse_->ComputeFlops());
00240   }
00241 
00243   virtual double ApplyFlops() const
00244   {
00245     return(ApplyFlops_);
00246   }
00247 
00249   virtual double ApplyInverseFlops() const
00250   {
00251     if (Inverse_ == Teuchos::null)
00252       return (0.0);
00253     else
00254       return(Inverse_->ApplyInverseFlops());
00255   }
00256   
00258   virtual ostream& Print(std::ostream& os) const;
00259 
00260 private:
00261   
00263   virtual int Extract(const Epetra_RowMatrix& Matrix_in);
00264 
00266   int NumRows_; 
00268   int NumVectors_; 
00270   Teuchos::RefCountPtr<Epetra_Map> Map_;
00272   Teuchos::RefCountPtr<Epetra_CrsMatrix> Matrix_;
00274   Teuchos::RefCountPtr<Epetra_MultiVector> LHS_;
00276   Teuchos::RefCountPtr<Epetra_MultiVector> RHS_;
00278   Epetra_IntSerialDenseVector GID_;
00280   bool IsInitialized_;
00282   bool IsComputed_;
00284   Teuchos::RefCountPtr<Epetra_Comm> SerialComm_;
00286   Teuchos::RefCountPtr<T> Inverse_;
00288   string Label_;
00289   Teuchos::ParameterList List_;
00290   double ApplyFlops_;
00291 
00292 };
00293 
00294 //==============================================================================
00295 template<typename T>
00296 Ifpack_SparseContainer<T>::
00297 Ifpack_SparseContainer(const int NumRows_in, const int NumVectors_in) :
00298   NumRows_(NumRows_in),
00299   NumVectors_(NumVectors_in),
00300   IsInitialized_(false),
00301   IsComputed_(false),
00302   ApplyFlops_(0.0)
00303 {
00304 
00305 #ifdef HAVE_MPI
00306   SerialComm_ = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_SELF) );
00307 #else
00308   SerialComm_ = Teuchos::rcp( new Epetra_SerialComm );
00309 #endif
00310 
00311 }
00312 
00313 //==============================================================================
00314 template<typename T>
00315 Ifpack_SparseContainer<T>::
00316 Ifpack_SparseContainer(const Ifpack_SparseContainer<T>& rhs) :
00317   NumRows_(rhs.NumRows()),
00318   NumVectors_(rhs.NumVectors()),
00319   IsInitialized_(rhs.IsInitialized()),
00320   IsComputed_(rhs.IsComputed())
00321 {
00322 
00323 #ifdef HAVE_MPI
00324   SerialComm_ = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_SELF) );
00325 #else
00326   SerialComm_ = Teuchos::rcp( new Epetra_SerialComm );
00327 #endif
00328 
00329   if (rhs.Map())
00330     Map_ = Teuchos::rcp( new Epetra_Map(*rhs.Map()) );
00331 
00332   if (rhs.Matrix())
00333     Matrix_ = Teuchos::rcp( new Epetra_CrsMatrix(*rhs.Matrix()) );
00334 
00335   if (rhs.LHS())
00336     LHS_ = Teuchos::rcp( new Epetra_MultiVector(*rhs.LHS()) );
00337 
00338   if (rhs.RHS())
00339     RHS_ = Teuchos::rcp( new Epetra_MultiVector(*rhs.RHS()) );
00340 
00341 }
00342 //==============================================================================
00343 template<typename T>
00344 Ifpack_SparseContainer<T>::~Ifpack_SparseContainer()
00345 {
00346   Destroy();
00347 }
00348 
00349 //==============================================================================
00350 template<typename T>
00351 int Ifpack_SparseContainer<T>::NumRows() const
00352 {
00353   if (IsInitialized() == false)
00354     return(0);
00355   else
00356     return(NumRows_);
00357 }
00358 
00359 //==============================================================================
00360 template<typename T>
00361 int Ifpack_SparseContainer<T>::Initialize()
00362 {
00363   
00364   if (IsInitialized_ == true)
00365     Destroy();
00366 
00367   IsInitialized_ = false;
00368 
00369   Map_ = Teuchos::rcp( new Epetra_Map(NumRows_,0,*SerialComm_) );
00370 
00371   LHS_ = Teuchos::rcp( new Epetra_MultiVector(*Map_,NumVectors_) );
00372   RHS_ = Teuchos::rcp( new Epetra_MultiVector(*Map_,NumVectors_) );
00373   GID_.Reshape(NumRows_,1);
00374 
00375   Matrix_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy,*Map_,0) );
00376 
00377   // create the inverse
00378   Inverse_ = Teuchos::rcp( new T(Matrix_.get()) );
00379 
00380   if (Inverse_ == Teuchos::null)
00381     IFPACK_CHK_ERR(-5);
00382 
00383   IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
00384 
00385   // Call Inverse_->Initialize() in Compute(). This saves
00386   // some time, because I can extract the diagonal blocks faster,
00387   // and only once.
00388 
00389   Label_ = "Ifpack_SparseContainer";
00390 
00391   IsInitialized_ = true;
00392   return(0);
00393   
00394 }
00395 
00396 //==============================================================================
00397 template<typename T>
00398 double& Ifpack_SparseContainer<T>::LHS(const int i, const int Vector)
00399 {
00400   return(((*LHS_)(Vector))->Values()[i]);
00401 }
00402   
00403 //==============================================================================
00404 template<typename T>
00405 double& Ifpack_SparseContainer<T>::RHS(const int i, const int Vector)
00406 {
00407   return(((*RHS_)(Vector))->Values()[i]);
00408 }
00409 
00410 //==============================================================================
00411 template<typename T>
00412 int Ifpack_SparseContainer<T>::
00413 SetMatrixElement(const int row, const int col, const double value)
00414 {
00415   if (!IsInitialized())
00416     IFPACK_CHK_ERR(-3); // problem not shaped yet
00417 
00418   if ((row < 0) || (row >= NumRows())) {
00419     IFPACK_CHK_ERR(-2); // not in range
00420   }
00421 
00422   if ((col < 0) || (col >= NumRows())) {
00423     IFPACK_CHK_ERR(-2); // not in range
00424   }
00425 
00426   int ierr = Matrix_->InsertGlobalValues((int)row,1,(double*)&value,(int*)&col);
00427   if (ierr < 0) {
00428     ierr = Matrix_->SumIntoGlobalValues((int)row,1,(double*)&value,(int*)&col);
00429     if (ierr < 0)
00430       IFPACK_CHK_ERR(-1);
00431   }
00432 
00433   return(0);
00434 
00435 }
00436 
00437 //==============================================================================
00438 template<typename T>
00439 int Ifpack_SparseContainer<T>::Compute(const Epetra_RowMatrix& Matrix_in)
00440 {
00441 
00442   IsComputed_ = false;
00443   if (!IsInitialized()) {
00444     IFPACK_CHK_ERR(Initialize()); 
00445   }
00446 
00447   // extract the submatrices
00448   IFPACK_CHK_ERR(Extract(Matrix_in));
00449 
00450   // initialize the inverse operator
00451   IFPACK_CHK_ERR(Inverse_->Initialize());
00452 
00453   // compute the inverse operator
00454   IFPACK_CHK_ERR(Inverse_->Compute());
00455 
00456   Label_ = "Ifpack_SparseContainer";
00457   
00458   IsComputed_ = true;
00459 
00460   return(0);
00461 }
00462 
00463 //==============================================================================
00464 template<typename T>
00465 int Ifpack_SparseContainer<T>::Apply()
00466 {
00467   if (IsComputed() == false) {
00468     IFPACK_CHK_ERR(-3); // not yet computed
00469   }
00470   
00471   IFPACK_CHK_ERR(Matrix_->Apply(*RHS_, *LHS_));
00472 
00473   ApplyFlops_ += 2 * Matrix_->NumGlobalNonzeros();
00474   return(0);
00475 }
00476 
00477 //==============================================================================
00478 template<typename T>
00479 int Ifpack_SparseContainer<T>::ApplyInverse()
00480 {
00481   if (!IsComputed())
00482     IFPACK_CHK_ERR(-1);
00483   
00484   IFPACK_CHK_ERR(Inverse_->ApplyInverse(*RHS_, *LHS_));
00485 
00486   return(0);
00487 }
00488  
00489 
00490 //==============================================================================
00491 template<typename T>
00492 int Ifpack_SparseContainer<T>::Destroy()
00493 {
00494   IsInitialized_ = false;
00495   IsComputed_ = false;
00496   return(0);
00497 }
00498 
00499 //==============================================================================
00500 template<typename T>
00501 int& Ifpack_SparseContainer<T>::ID(const int i)
00502 {
00503   return(GID_[i]);
00504 }
00505 
00506 //==============================================================================
00507 template<typename T>
00508 int Ifpack_SparseContainer<T>::
00509 SetParameters(Teuchos::ParameterList& List)
00510 {
00511   List_ = List;
00512   return(0);
00513 }
00514 
00515 //==============================================================================
00516 // FIXME: optimize performances of this guy...
00517 template<typename T>
00518 int Ifpack_SparseContainer<T>::Extract(const Epetra_RowMatrix& Matrix_in)
00519 {
00520 
00521   for (int j = 0 ; j < NumRows_ ; ++j) {
00522     // be sure that the user has set all the ID's
00523     if (ID(j) == -1)
00524       IFPACK_CHK_ERR(-1);
00525     // be sure that all are local indices
00526     if (ID(j) > Matrix_in.NumMyRows())
00527       IFPACK_CHK_ERR(-1);
00528   }
00529 
00530   int Length = Matrix_in.MaxNumEntries();
00531   std::vector<double> Values;
00532   Values.resize(Length);
00533   std::vector<int> Indices;
00534   Indices.resize(Length);
00535 
00536   for (int j = 0 ; j < NumRows_ ; ++j) {
00537 
00538     int LRID = ID(j);
00539 
00540     int NumEntries;
00541 
00542     int ierr = 
00543       Matrix_in.ExtractMyRowCopy(LRID, Length, NumEntries, 
00544              &Values[0], &Indices[0]);
00545     IFPACK_CHK_ERR(ierr);
00546 
00547     for (int k = 0 ; k < NumEntries ; ++k) {
00548 
00549       int LCID = Indices[k];
00550 
00551       // skip off-processor elements
00552       if (LCID >= Matrix_in.NumMyRows()) 
00553   continue;
00554 
00555       // for local column IDs, look for each ID in the list
00556       // of columns hosted by this object
00557       // FIXME: use STL
00558       int jj = -1;
00559       for (int kk = 0 ; kk < NumRows_ ; ++kk)
00560   if (ID(kk) == LCID)
00561     jj = kk;
00562 
00563       if (jj != -1)
00564   SetMatrixElement(j,jj,Values[k]);
00565 
00566     }
00567   }
00568 
00569   IFPACK_CHK_ERR(Matrix_->FillComplete());
00570 
00571   return(0);
00572 }
00573 
00574 //==============================================================================
00575 template<typename T>
00576 ostream& Ifpack_SparseContainer<T>::Print(ostream & os) const
00577 {
00578   os << "================================================================================" << endl;
00579   os << "Ifpack_SparseContainer" << endl;
00580   os << "Number of rows          = " << NumRows() << endl;
00581   os << "Number of vectors       = " << NumVectors() << endl;
00582   os << "IsInitialized()         = " << IsInitialized() << endl;
00583   os << "IsComputed()            = " << IsComputed() << endl;
00584   os << "Flops in Initialize()   = " << InitializeFlops() << endl; 
00585   os << "Flops in Compute()      = " << ComputeFlops() << endl; 
00586   os << "Flops in ApplyInverse() = " << ApplyInverseFlops() << endl; 
00587   os << "================================================================================" << endl;
00588   os << endl;
00589 
00590   return(os);
00591 }
00592 #endif // IFPACK_SPARSECONTAINER_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:35 2011 for Ifpack Package Browser (Single Doxygen Collection) by  doxygen 1.6.3