Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_SparseContainer_def.hpp
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) 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 IFPACK2_SPARSECONTAINER_DEF_HPP
00031 #define IFPACK2_SPARSECONTAINER_DEF_HPP
00032 
00033 #include "Ifpack2_SparseContainer_decl.hpp"
00034 #ifdef HAVE_MPI
00035 #include <mpi.h>
00036 #include "Teuchos_DefaultMpiComm.hpp"
00037 #else
00038 #include "Teuchos_DefaultSerialComm.hpp"
00039 #endif
00040 
00041 namespace Ifpack2 {
00042 
00043 //==============================================================================
00044 template<class MatrixType, class InverseType>
00045 SparseContainer<MatrixType,InverseType>::SparseContainer(const size_t NumRows, const size_t NumVectors) :
00046   NumRows_(NumRows),
00047   NumVectors_(NumVectors),
00048   IsInitialized_(false),
00049   IsComputed_(false)
00050 {
00051 #ifdef HAVE_MPI
00052   LocalComm_ = Teuchos::rcp(new Teuchos::MpiComm<int>(Teuchos::opaqueWrapper((MPI_Comm)MPI_COMM_SELF)));
00053 #else
00054   LocalComm_ = Teuchos::rcp( new Teuchos::SerialComm<int>() );
00055 #endif
00056 }
00057 
00058 //==============================================================================
00059 template<class MatrixType, class InverseType>
00060 SparseContainer<MatrixType,InverseType>::SparseContainer(const SparseContainer<MatrixType,InverseType>& rhs)
00061 {
00062   // nobody should ever call this  
00063 }
00064 
00065 //==============================================================================
00066 template<class MatrixType, class InverseType>
00067 SparseContainer<MatrixType,InverseType>::~SparseContainer()
00068 {
00069   destroy();
00070 }
00071 
00072 //==============================================================================
00073 template<class MatrixType, class InverseType>
00074 size_t SparseContainer<MatrixType,InverseType>::getNumRows() const
00075 {
00076   if (isInitialized()) return NumRows_;
00077   else return 0;
00078 }
00079 
00080 //==============================================================================
00081 // Sets the number of vectors for X/Y
00082 template<class MatrixType, class InverseType>
00083 void SparseContainer<MatrixType,InverseType>::setNumVectors(const size_t NumVectors_in)
00084 {
00085   TEUCHOS_TEST_FOR_EXCEPTION(NumVectors_in<=0, std::runtime_error, "Ifpack2::SparseContainer::setNumVectors NumVectors cannot be negative.");  
00086 
00087   // Do nothing if the vectors are the right size.  Always do something if we're not initialized.
00088   if(!IsInitialized_  || NumVectors_!=NumVectors_in) {
00089     NumVectors_=NumVectors_in;
00090     X_ = Teuchos::rcp( new Tpetra::MultiVector<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode>(Map_,NumVectors_) );
00091     Y_ = Teuchos::rcp( new Tpetra::MultiVector<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode>(Map_,NumVectors_) );
00092   }
00093 
00094 }
00095 
00096 //==============================================================================
00097 // Returns the ID associated to local row i. 
00098 template<class MatrixType, class InverseType>
00099 typename MatrixType::local_ordinal_type & SparseContainer<MatrixType,InverseType>::ID(const size_t i)
00100 {
00101   return GID_[i];
00102 }
00103 
00104 //==============================================================================
00105 // Returns  true is the container has been successfully initialized.
00106 template<class MatrixType, class InverseType>
00107 bool SparseContainer<MatrixType,InverseType>::isInitialized() const
00108 {
00109   return IsInitialized_;
00110 }
00111 
00112 //==============================================================================
00113 // Returns  true is the container has been successfully computed.
00114 template<class MatrixType, class InverseType>
00115 bool SparseContainer<MatrixType,InverseType>::isComputed() const
00116 {
00117   return IsComputed_;
00118 }
00119 
00120 //==============================================================================
00121 // Sets all necessary parameters.
00122 template<class MatrixType, class InverseType>
00123 void SparseContainer<MatrixType,InverseType>::setParameters(const Teuchos::ParameterList& List)
00124 {
00125   List_ = List;
00126 }
00127  
00128 //==============================================================================
00129 template<class MatrixType, class InverseType>
00130 void SparseContainer<MatrixType,InverseType>::initialize()
00131 {
00132   if(IsInitialized_) destroy();
00133   IsInitialized_=false;
00134 
00135   // Allocations
00136   Map_ = Teuchos::rcp( new Tpetra::Map<InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode>(NumRows_,0,LocalComm_) );
00137   Matrix_ = Teuchos::rcp( new Tpetra::CrsMatrix<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode>(Map_,0) );
00138   GID_.resize(NumRows_);  
00139   setNumVectors(NumVectors_);
00140 
00141   // create the inverse
00142   Inverse_ = Teuchos::rcp( new InverseType(Matrix_) );
00143   TEUCHOS_TEST_FOR_EXCEPTION( Inverse_ == Teuchos::null, std::runtime_error, "Ifpack2::SparseContainer::initialize error in inverse constructor.");  
00144   Inverse_->setParameters(List_);
00145 
00146   // Note from IFPACK: Call Inverse_->Initialize() in Compute(). This saves
00147   // some time, because the diagonal blocks can be extracted faster,
00148   // and only once.
00149   IsInitialized_ = true;
00150 }
00151 
00152 //==============================================================================
00153 // Finalizes the linear system matrix and prepares for the application of the inverse.
00154 template<class MatrixType, class InverseType>
00155 void SparseContainer<MatrixType,InverseType>::compute(const Teuchos::RCP<const Tpetra::RowMatrix<MatrixScalar,MatrixLocalOrdinal,MatrixGlobalOrdinal,MatrixNode> >& Matrix)
00156 {
00157   IsComputed_=false;
00158   TEUCHOS_TEST_FOR_EXCEPTION( !IsInitialized_, std::runtime_error, "Ifpack2::SparseContainer::compute please call initialize first.");  
00159 
00160   // extract the submatrices
00161   extract(Matrix);
00162 
00163   // initialize & compute the inverse operator
00164   Inverse_->initialize();
00165   Inverse_->compute();
00166   IsComputed_=true;
00167 }
00168 
00169 //==============================================================================
00170 // Computes Y = alpha * M^{-1} X + beta*Y
00171 template<class MatrixType, class InverseType>
00172 void SparseContainer<MatrixType,InverseType>::apply(const Tpetra::MultiVector<MatrixScalar,MatrixLocalOrdinal,MatrixGlobalOrdinal,MatrixNode>& X,
00173                 Tpetra::MultiVector<MatrixScalar,MatrixLocalOrdinal,MatrixGlobalOrdinal,MatrixNode>& Y,
00174                 Teuchos::ETransp mode,
00175                 MatrixScalar alpha,
00176                 MatrixScalar beta)
00177 {
00178  TEUCHOS_TEST_FOR_EXCEPTION( !IsComputed_, std::runtime_error, "Ifpack2::SparseContainer::apply compute has not been called.");  
00179  TEUCHOS_TEST_FOR_EXCEPTION( X.getNumVectors() != Y.getNumVectors(), std::runtime_error, "Ifpack2::SparseContainer::apply X/Y have different numbers of vectors.");   
00180 
00181  // This is a noop if our internal vectors are already the right size.
00182  setNumVectors(X.getNumVectors());
00183 
00184  // Pull the Data Pointers
00185  Teuchos::ArrayRCP<Teuchos::ArrayRCP<const MatrixScalar> >  global_x_ptr = X.get2dView();
00186  Teuchos::ArrayRCP<Teuchos::ArrayRCP<InverseScalar> >       local_x_ptr  = X_->get2dViewNonConst();
00187  Teuchos::ArrayRCP<Teuchos::ArrayRCP<const InverseScalar> > local_y_ptr  = Y_->get2dView();
00188  Teuchos::ArrayRCP<Teuchos::ArrayRCP<MatrixScalar> >        global_y_ptr = Y.get2dViewNonConst();
00189 
00190  // Copy in
00191  for(size_t j=0; j < NumRows_; j++){
00192    MatrixLocalOrdinal LID = ID(j);
00193    for (size_t k = 0 ; k < NumVectors_ ; k++) 
00194      local_x_ptr[k][j] = global_x_ptr[k][LID];   
00195  }
00196 
00197  // Apply inverse
00198  Inverse_->apply(*X_, *Y_);
00199  
00200  // copy out into solution vector Y
00201  for (size_t j = 0 ; j < NumRows_ ; j++) {
00202    MatrixLocalOrdinal LID = ID(j);
00203    for (size_t k = 0 ; k < NumVectors_ ; k++) 
00204      global_y_ptr[k][LID] = alpha * local_y_ptr[k][j] + beta * global_y_ptr[k][LID];
00205  }
00206 
00207 }
00208 
00209 
00210 //==============================================================================
00211 // Computes Y = alpha * diag(D) * M^{-1} (diag(D) X) + beta*Y
00212 template<class MatrixType, class InverseType>
00213 void SparseContainer<MatrixType,InverseType>::weightedApply(const Tpetra::MultiVector<MatrixScalar,MatrixLocalOrdinal,MatrixGlobalOrdinal,MatrixNode>& X,
00214                   Tpetra::MultiVector<MatrixScalar,MatrixLocalOrdinal,MatrixGlobalOrdinal,MatrixNode>& Y,
00215                   const Tpetra::Vector<MatrixScalar,MatrixLocalOrdinal,MatrixGlobalOrdinal,MatrixNode>& D,
00216                   Teuchos::ETransp mode,
00217                   MatrixScalar alpha,
00218                   MatrixScalar beta)
00219 {
00220   TEUCHOS_TEST_FOR_EXCEPTION( !IsComputed_, std::runtime_error, "Ifpack2::SparseContainer::apply compute has not been called.");  
00221  TEUCHOS_TEST_FOR_EXCEPTION( X.getNumVectors() != Y.getNumVectors(), std::runtime_error, "Ifpack2::SparseContainer::apply X/Y have different numbers of vectors.");   
00222 
00223  // This is a noop if our internal vectors are already the right size.
00224  setNumVectors(X.getNumVectors());
00225 
00226  // Pull the data pointers
00227  Teuchos::ArrayRCP<const MatrixScalar >                     global_d_ptr = D.getData(0);
00228  Teuchos::ArrayRCP<Teuchos::ArrayRCP<const MatrixScalar> >  global_x_ptr = X.get2dView();
00229  Teuchos::ArrayRCP<Teuchos::ArrayRCP<InverseScalar> >       local_x_ptr  = X_->get2dViewNonConst();
00230  Teuchos::ArrayRCP<Teuchos::ArrayRCP<const InverseScalar> > local_y_ptr  = Y_->get2dView();
00231  Teuchos::ArrayRCP<Teuchos::ArrayRCP<MatrixScalar> >        global_y_ptr = Y.get2dViewNonConst();
00232 
00233  // Copy in
00234  for(size_t j=0; j < NumRows_; j++){
00235    MatrixLocalOrdinal LID = ID(j);
00236    for (size_t k = 0 ; k < NumVectors_ ; k++) 
00237      local_x_ptr[k][j] = global_d_ptr[LID] * global_x_ptr[k][LID];   
00238  }
00239 
00240  // Apply inverse
00241  Inverse_->apply(*X_, *Y_);
00242  
00243  // copy out into solution vector Y
00244  for (size_t j = 0 ; j < NumRows_ ; j++) {
00245    MatrixLocalOrdinal LID = ID(j);
00246    for (size_t k = 0 ; k < NumVectors_ ; k++) 
00247      global_y_ptr[k][LID] = alpha * global_d_ptr[LID] * local_y_ptr[k][j] + beta * global_y_ptr[k][LID];
00248  }
00249 }
00250 
00251 //==============================================================================
00252 // Destroys all data.
00253 template<class MatrixType, class InverseType>
00254 void SparseContainer<MatrixType,InverseType>::destroy()
00255 {
00256   IsInitialized_ = false;
00257   IsComputed_ = false;
00258 }
00259 
00260 //============================================================================== 
00261 // Prints basic information on iostream. This function is used by operator<<
00262 template<class MatrixType, class InverseType>
00263 std::ostream& SparseContainer<MatrixType,InverseType>::print(std::ostream& os) const
00264 {
00265   Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
00266   fos.setOutputToRootOnly(0);
00267   describe(fos);
00268   return(os);
00269 }
00270 
00271 
00272 //==============================================================================
00273 template<class MatrixType, class InverseType>
00274 std::string SparseContainer<MatrixType,InverseType>::description() const
00275 {
00276   std::ostringstream oss;
00277   oss << Teuchos::Describable::description();
00278   if (isInitialized()) {
00279     if (isComputed()) {
00280       oss << "{status = initialized, computed";
00281     }
00282     else {
00283       oss << "{status = initialized, not computed";
00284     }
00285   }
00286   else {
00287     oss << "{status = not initialized, not computed";
00288   }
00289 
00290   oss << "}";
00291   return oss.str();
00292 }
00293 
00294 //==============================================================================
00295 template<class MatrixType, class InverseType>
00296 void SparseContainer<MatrixType,InverseType>::describe(Teuchos::FancyOStream &os, const Teuchos::EVerbosityLevel verbLevel) const
00297 {
00298   using std::endl;
00299   if(verbLevel==Teuchos::VERB_NONE) return;
00300   os << "================================================================================" << endl;
00301   os << "Ifpack2_SparseContainer" << endl;
00302   os << "Number of rows          = " << NumRows_ << endl;
00303   os << "Number of vectors       = " << NumVectors_ << endl;
00304   os << "isInitialized()         = " << IsInitialized_ << endl;
00305   os << "isComputed()            = " << IsComputed_ << endl;
00306   os << "================================================================================" << endl;
00307   os << endl;
00308 }
00309 
00310 //============================================================================== 
00311 // Extract the submatrices identified by the ID set int ID().
00312 template<class MatrixType, class InverseType>
00313 void SparseContainer<MatrixType,InverseType>::extract(const Teuchos::RCP<const Tpetra::RowMatrix<MatrixScalar,MatrixLocalOrdinal,MatrixGlobalOrdinal,MatrixNode> >& Matrix_in) 
00314 {
00315   size_t MatrixInNumRows= Matrix_in->getNodeNumRows();
00316 
00317   // Sanity checking
00318   for (size_t j = 0 ; j < NumRows_ ; ++j) {
00319     TEUCHOS_TEST_FOR_EXCEPTION( GID_[j] < 0 || (size_t) GID_[j] >= MatrixInNumRows, std::runtime_error, "Ifpack2::SparseContainer::applyInverse compute has not been called.");  
00320   }  
00321 
00322   int Length = Matrix_in->getNodeMaxNumRowEntries();
00323   Teuchos::Array<MatrixScalar>         Values;
00324   Teuchos::Array<MatrixLocalOrdinal>   Indices;
00325   Teuchos::Array<InverseScalar>        Values_insert;
00326   Teuchos::Array<InverseGlobalOrdinal> Indices_insert;
00327 
00328   Values.resize(Length);
00329   Indices.resize(Length);
00330   Values_insert.resize(Length);
00331   Indices_insert.resize(Length);
00332 
00333   for (size_t j = 0 ; j < NumRows_ ; ++j) {
00334     MatrixLocalOrdinal LRID = ID(j);
00335     size_t NumEntries;
00336 
00337     Matrix_in->getLocalRowCopy(LRID,Indices(),Values(),NumEntries);
00338 
00339     size_t num_entries_found=0;
00340     for (size_t k = 0 ; k < NumEntries ; ++k) {
00341       MatrixLocalOrdinal LCID = Indices[k];
00342 
00343       // skip off-processor elements
00344       if ((size_t)LCID >= MatrixInNumRows) 
00345   continue;
00346 
00347       // for local column IDs, look for each ID in the list
00348       // of columns hosted by this object
00349       // FIXME: use STL
00350       InverseLocalOrdinal jj = -1;
00351       for (size_t kk = 0 ; kk < NumRows_ ; ++kk)
00352   if (ID(kk) == LCID)
00353     jj = kk;
00354 
00355       if (jj != -1) {
00356   Indices_insert[num_entries_found] = Map_->getGlobalElement(jj);
00357   Values_insert[num_entries_found]  = Values[k];
00358   num_entries_found++;
00359       }
00360 
00361     }
00362     Matrix_->insertGlobalValues(j,Indices_insert(0,num_entries_found),Values_insert(0,num_entries_found));
00363   }
00364 
00365   Matrix_->fillComplete();
00366 }
00367 
00368 
00369 } // namespace Ifpack2
00370 #endif // IFPACK2_SPARSECONTAINER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends