IFPACK Development
Ifpack_SparseContainer.h
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        {
00116        NumVectors_=NumVectors_in;
00117        LHS_=Teuchos::rcp(new Epetra_MultiVector(*Map_,NumVectors_));
00118        RHS_=Teuchos::rcp(new Epetra_MultiVector(*Map_,NumVectors_));
00119        }
00120    return(0);
00121    }
00122 
00124   virtual double& LHS(const int i, const int Vector = 0);
00125   
00127   virtual double& RHS(const int i, const int Vector = 0);
00128 
00130 
00139   virtual int& ID(const int i);
00140 
00142   virtual int SetMatrixElement(const int row, const int col,
00143                    const double value);
00144 
00145 
00147   virtual bool IsInitialized() const
00148   {
00149     return(IsInitialized_);
00150   }
00151 
00153   virtual bool IsComputed() const
00154   {
00155     return(IsComputed_);
00156   }
00157 
00159   virtual int SetParameters(Teuchos::ParameterList& List);
00160 
00162   virtual const char* Label() const
00163   {
00164     return(Label_.c_str());
00165   }
00166   
00168   Teuchos::RCP<const Epetra_Map> Map() const
00169   {
00170     return(Map_);
00171   }
00172 
00174   Teuchos::RCP<const Epetra_MultiVector> LHS() const
00175   {
00176     return(LHS_);
00177   }
00178 
00180   Teuchos::RCP<const Epetra_MultiVector> RHS() const
00181   {
00182     return(RHS_);
00183   }
00184 
00186   Teuchos::RCP<const Epetra_CrsMatrix> Matrix() const
00187   {
00188     return(Matrix_);
00189   }
00190 
00192   const Epetra_IntSerialDenseVector& ID() const
00193   {
00194     return(GID_);
00195   }
00196 
00198   Teuchos::RCP<const T> Inverse() const
00199   {
00200     return(Inverse_);
00201   }
00203 
00205 
00212   virtual int Initialize();
00214   virtual int Compute(const Epetra_RowMatrix& Matrix_in);
00216   virtual int Apply();
00217 
00219   virtual int ApplyInverse();
00220 
00222 
00224 
00225   virtual int Destroy();
00227 
00229   virtual double InitializeFlops() const
00230   {
00231     if (Inverse_ == Teuchos::null)
00232       return (0.0);
00233     else
00234       return(Inverse_->InitializeFlops());
00235   }
00236 
00238   virtual double ComputeFlops() const
00239   {
00240     if (Inverse_ == Teuchos::null)
00241       return (0.0);
00242     else
00243       return(Inverse_->ComputeFlops());
00244   }
00245 
00247   virtual double ApplyFlops() const
00248   {
00249     return(ApplyFlops_);
00250   }
00251 
00253   virtual double ApplyInverseFlops() const
00254   {
00255     if (Inverse_ == Teuchos::null)
00256       return (0.0);
00257     else
00258       return(Inverse_->ApplyInverseFlops());
00259   }
00260   
00262   virtual ostream& Print(std::ostream& os) const;
00263 
00264 private:
00265   
00267   virtual int Extract(const Epetra_RowMatrix& Matrix_in);
00268 
00270   int NumRows_; 
00272   int NumVectors_; 
00274   Teuchos::RefCountPtr<Epetra_Map> Map_;
00276   Teuchos::RefCountPtr<Epetra_CrsMatrix> Matrix_;
00278   Teuchos::RefCountPtr<Epetra_MultiVector> LHS_;
00280   Teuchos::RefCountPtr<Epetra_MultiVector> RHS_;
00282   Epetra_IntSerialDenseVector GID_;
00284   bool IsInitialized_;
00286   bool IsComputed_;
00288   Teuchos::RefCountPtr<Epetra_Comm> SerialComm_;
00290   Teuchos::RefCountPtr<T> Inverse_;
00292   string Label_;
00293   Teuchos::ParameterList List_;
00294   double ApplyFlops_;
00295 
00296 };
00297 
00298 //==============================================================================
00299 template<typename T>
00300 Ifpack_SparseContainer<T>::
00301 Ifpack_SparseContainer(const int NumRows_in, const int NumVectors_in) :
00302   NumRows_(NumRows_in),
00303   NumVectors_(NumVectors_in),
00304   IsInitialized_(false),
00305   IsComputed_(false),
00306   ApplyFlops_(0.0)
00307 {
00308 
00309 #ifdef HAVE_MPI
00310   SerialComm_ = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_SELF) );
00311 #else
00312   SerialComm_ = Teuchos::rcp( new Epetra_SerialComm );
00313 #endif
00314 
00315 }
00316 
00317 //==============================================================================
00318 template<typename T>
00319 Ifpack_SparseContainer<T>::
00320 Ifpack_SparseContainer(const Ifpack_SparseContainer<T>& rhs) :
00321   NumRows_(rhs.NumRows()),
00322   NumVectors_(rhs.NumVectors()),
00323   IsInitialized_(rhs.IsInitialized()),
00324   IsComputed_(rhs.IsComputed())
00325 {
00326 
00327 #ifdef HAVE_MPI
00328   SerialComm_ = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_SELF) );
00329 #else
00330   SerialComm_ = Teuchos::rcp( new Epetra_SerialComm );
00331 #endif
00332 
00333   if (rhs.Map())
00334     Map_ = Teuchos::rcp( new Epetra_Map(*rhs.Map()) );
00335 
00336   if (rhs.Matrix())
00337     Matrix_ = Teuchos::rcp( new Epetra_CrsMatrix(*rhs.Matrix()) );
00338 
00339   if (rhs.LHS())
00340     LHS_ = Teuchos::rcp( new Epetra_MultiVector(*rhs.LHS()) );
00341 
00342   if (rhs.RHS())
00343     RHS_ = Teuchos::rcp( new Epetra_MultiVector(*rhs.RHS()) );
00344 
00345 }
00346 //==============================================================================
00347 template<typename T>
00348 Ifpack_SparseContainer<T>::~Ifpack_SparseContainer()
00349 {
00350   Destroy();
00351 }
00352 
00353 //==============================================================================
00354 template<typename T>
00355 int Ifpack_SparseContainer<T>::NumRows() const
00356 {
00357   if (IsInitialized() == false)
00358     return(0);
00359   else
00360     return(NumRows_);
00361 }
00362 
00363 //==============================================================================
00364 template<typename T>
00365 int Ifpack_SparseContainer<T>::Initialize()
00366 {
00367   
00368   if (IsInitialized_ == true)
00369     Destroy();
00370 
00371   IsInitialized_ = false;
00372 
00373   Map_ = Teuchos::rcp( new Epetra_Map(NumRows_,0,*SerialComm_) );
00374 
00375   LHS_ = Teuchos::rcp( new Epetra_MultiVector(*Map_,NumVectors_) );
00376   RHS_ = Teuchos::rcp( new Epetra_MultiVector(*Map_,NumVectors_) );
00377   GID_.Reshape(NumRows_,1);
00378 
00379   Matrix_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy,*Map_,0) );
00380 
00381   // create the inverse
00382   Inverse_ = Teuchos::rcp( new T(Matrix_.get()) );
00383 
00384   if (Inverse_ == Teuchos::null)
00385     IFPACK_CHK_ERR(-5);
00386 
00387   IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
00388 
00389   // Call Inverse_->Initialize() in Compute(). This saves
00390   // some time, because I can extract the diagonal blocks faster,
00391   // and only once.
00392 
00393   Label_ = "Ifpack_SparseContainer";
00394 
00395   IsInitialized_ = true;
00396   return(0);
00397   
00398 }
00399 
00400 //==============================================================================
00401 template<typename T>
00402 double& Ifpack_SparseContainer<T>::LHS(const int i, const int Vector)
00403 {
00404   return(((*LHS_)(Vector))->Values()[i]);
00405 }
00406   
00407 //==============================================================================
00408 template<typename T>
00409 double& Ifpack_SparseContainer<T>::RHS(const int i, const int Vector)
00410 {
00411   return(((*RHS_)(Vector))->Values()[i]);
00412 }
00413 
00414 //==============================================================================
00415 template<typename T>
00416 int Ifpack_SparseContainer<T>::
00417 SetMatrixElement(const int row, const int col, const double value)
00418 {
00419   if (!IsInitialized())
00420     IFPACK_CHK_ERR(-3); // problem not shaped yet
00421 
00422   if ((row < 0) || (row >= NumRows())) {
00423     IFPACK_CHK_ERR(-2); // not in range
00424   }
00425 
00426   if ((col < 0) || (col >= NumRows())) {
00427     IFPACK_CHK_ERR(-2); // not in range
00428   }
00429 
00430   int ierr = Matrix_->InsertGlobalValues((int)row,1,(double*)&value,(int*)&col);
00431   if (ierr < 0) {
00432     ierr = Matrix_->SumIntoGlobalValues((int)row,1,(double*)&value,(int*)&col);
00433     if (ierr < 0)
00434       IFPACK_CHK_ERR(-1);
00435   }
00436 
00437   return(0);
00438 
00439 }
00440 
00441 //==============================================================================
00442 template<typename T>
00443 int Ifpack_SparseContainer<T>::Compute(const Epetra_RowMatrix& Matrix_in)
00444 {
00445 
00446   IsComputed_ = false;
00447   if (!IsInitialized()) {
00448     IFPACK_CHK_ERR(Initialize()); 
00449   }
00450 
00451   // extract the submatrices
00452   IFPACK_CHK_ERR(Extract(Matrix_in));
00453 
00454   // initialize the inverse operator
00455   IFPACK_CHK_ERR(Inverse_->Initialize());
00456 
00457   // compute the inverse operator
00458   IFPACK_CHK_ERR(Inverse_->Compute());
00459 
00460   Label_ = "Ifpack_SparseContainer";
00461   
00462   IsComputed_ = true;
00463 
00464   return(0);
00465 }
00466 
00467 //==============================================================================
00468 template<typename T>
00469 int Ifpack_SparseContainer<T>::Apply()
00470 {
00471   if (IsComputed() == false) {
00472     IFPACK_CHK_ERR(-3); // not yet computed
00473   }
00474   
00475   IFPACK_CHK_ERR(Matrix_->Apply(*RHS_, *LHS_));
00476 
00477   ApplyFlops_ += 2 * Matrix_->NumGlobalNonzeros();
00478   return(0);
00479 }
00480 
00481 //==============================================================================
00482 template<typename T>
00483 int Ifpack_SparseContainer<T>::ApplyInverse()
00484 {
00485   if (!IsComputed())
00486     IFPACK_CHK_ERR(-1);
00487   
00488   IFPACK_CHK_ERR(Inverse_->ApplyInverse(*RHS_, *LHS_));
00489 
00490   return(0);
00491 }
00492  
00493 
00494 //==============================================================================
00495 template<typename T>
00496 int Ifpack_SparseContainer<T>::Destroy()
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_in)
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_in.NumMyRows())
00531       IFPACK_CHK_ERR(-1);
00532   }
00533 
00534   int Length = Matrix_in.MaxNumEntries();
00535   std::vector<double> Values;
00536   Values.resize(Length);
00537   std::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_in.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_in.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 // IFPACK_SPARSECONTAINER_H
 All Classes Files Functions Variables Enumerations Friends