Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_SparsityFilter_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_SPARSITYFILTER_DEF_HPP
00031 #define IFPACK2_SPARSITYFILTER_DEF_HPP
00032 #include "Ifpack2_SparsityFilter_decl.hpp"
00033 #include <vector>
00034 
00035 #include "Tpetra_ConfigDefs.hpp"
00036 #include "Tpetra_RowMatrix.hpp"
00037 #include "Tpetra_Map.hpp"
00038 #include "Tpetra_MultiVector.hpp"
00039 #include "Tpetra_Vector.hpp"
00040 
00041 #include <vector>
00042 #include <algorithm>
00043 
00044 namespace Ifpack2 {
00045 
00046 //==========================================================================
00047 template<class MatrixType>
00048 SparsityFilter<MatrixType>::SparsityFilter(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Matrix, 
00049              size_t AllowedNumEntries,
00050              LocalOrdinal AllowedBandwidth):
00051   A_(Matrix),
00052   AllowedNumEntries_(AllowedNumEntries),
00053   AllowedBandwidth_(AllowedBandwidth),
00054   NumRows_(0),
00055   NumNonzeros_(0),
00056   MaxNumEntries_(0),
00057   MaxNumEntriesA_(0)
00058 {
00059 
00060  // use this filter only on serial matrices
00061   if (A_->getComm()->getSize() != 1 || A_->getNodeNumRows() != A_->getGlobalNumRows()) {
00062     throw std::runtime_error("Ifpack2::SparsityFilter can be used with Comm().getSize() == 1 only. This class is a tool for Ifpack2_AdditiveSchwarz, and it is not meant to be used otherwise.");
00063   }
00064 
00065 
00066   // localized matrix has all the local rows of Matrix
00067   NumRows_ = A_->getNodeNumRows();
00068 
00069   // NodeNumEntries_ will contain the actual number of nonzeros
00070   // for each localized row (that is, without external nodes,
00071   // and always with the diagonal entry)
00072   NumEntries_.resize(NumRows_);
00073 
00074   // tentative value for MaxNumEntries. This is the number of
00075   // nonzeros in the local matrix
00076   MaxNumEntries_  = A_->getNodeMaxNumRowEntries();
00077   MaxNumEntriesA_ = A_->getNodeMaxNumRowEntries();
00078 
00079   // ExtractMyRowCopy() will use these vectors
00080   Indices_.resize(MaxNumEntries_);
00081   Values_.resize(MaxNumEntries_);
00082 
00083   size_t ActualMaxNumEntries = 0;
00084   for (size_t i = 0 ; i < NumRows_ ; ++i) {
00085     NumEntries_[i] = MaxNumEntriesA_;
00086     size_t Nnz;
00087     A_->getLocalRowCopy(i,Indices_,Values_,Nnz);
00088 
00089     NumNonzeros_ += Nnz;
00090     NumEntries_[i] = Nnz;
00091     if (Nnz > ActualMaxNumEntries)
00092       ActualMaxNumEntries = Nnz;
00093   }
00094   
00095   MaxNumEntries_ = ActualMaxNumEntries;
00096 }
00097 
00098 //=========================================================================
00099 template<class MatrixType>
00100 SparsityFilter<MatrixType>::~SparsityFilter() { }
00101 
00102 //==========================================================================
00103 template<class MatrixType>
00104 const Teuchos::RCP<const Teuchos::Comm<int> > & SparsityFilter<MatrixType>::getComm() const
00105 {
00106   return A_->getComm();
00107 }
00108 
00109 //==========================================================================
00110 template<class MatrixType>
00111 Teuchos::RCP <typename MatrixType::node_type> SparsityFilter<MatrixType>::getNode() const
00112 {
00113   return A_->getNode();
00114 }
00115 
00116 //==========================================================================
00117 template<class MatrixType>
00118 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00119                                      typename MatrixType::global_ordinal_type,
00120                                      typename MatrixType::node_type> >&
00121 SparsityFilter<MatrixType>::getRowMap() const
00122 {
00123   return A_->getRowMap();
00124 }
00125 
00126 //==========================================================================
00127 template<class MatrixType>
00128 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00129                                      typename MatrixType::global_ordinal_type,
00130                                      typename MatrixType::node_type> >&
00131 SparsityFilter<MatrixType>::getColMap() const
00132 {
00133   return A_->getColMap();
00134 }
00135 
00136 //==========================================================================
00137 template<class MatrixType>
00138 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00139                                      typename MatrixType::global_ordinal_type,
00140                                      typename MatrixType::node_type> >&
00141 SparsityFilter<MatrixType>::getDomainMap() const
00142 {
00143   return A_->getDomainMap();
00144 }
00145 
00146 //==========================================================================
00147 template<class MatrixType>
00148 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00149                                      typename MatrixType::global_ordinal_type,
00150                                      typename MatrixType::node_type> >&
00151 SparsityFilter<MatrixType>::getRangeMap() const
00152 {
00153   return A_->getRangeMap();
00154 }
00155 
00156 //==========================================================================
00157 template<class MatrixType>
00158 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type,
00159                                      typename MatrixType::global_ordinal_type,
00160                                      typename MatrixType::node_type> >
00161 SparsityFilter<MatrixType>::getGraph() const
00162 {
00163   throw std::runtime_error("Ifpack2::SparsityFilter: does not support getGraph.");
00164 }
00165 
00166 //==========================================================================
00167 template<class MatrixType>
00168 global_size_t SparsityFilter<MatrixType>::getGlobalNumRows() const
00169 {
00170   return NumRows_;
00171 }
00172 
00173 //==========================================================================
00174 template<class MatrixType>
00175 global_size_t SparsityFilter<MatrixType>::getGlobalNumCols() const
00176 {
00177   return NumRows_;
00178 }
00179 
00180 //==========================================================================
00181 template<class MatrixType>
00182 size_t SparsityFilter<MatrixType>::getNodeNumRows() const
00183 {
00184   return NumRows_;
00185 }
00186 
00187 //==========================================================================
00188  
00189 template<class MatrixType>
00190 size_t SparsityFilter<MatrixType>::getNodeNumCols() const
00191 {
00192   return NumRows_;
00193 }
00194 
00195 //==========================================================================  
00196 template<class MatrixType>
00197 typename MatrixType::global_ordinal_type SparsityFilter<MatrixType>::getIndexBase() const
00198 {
00199   return A_->getIndexBase();
00200 }
00201 
00202 //========================================================================== 
00203 template<class MatrixType>
00204 global_size_t SparsityFilter<MatrixType>::getGlobalNumEntries() const
00205 {
00206   return NumNonzeros_;
00207 }
00208 
00209 //========================================================================== 
00210 template<class MatrixType>
00211 size_t SparsityFilter<MatrixType>::getNodeNumEntries() const
00212 {
00213   return NumNonzeros_;
00214 }
00215 
00216 //==========================================================================
00217 template<class MatrixType> 
00218 size_t SparsityFilter<MatrixType>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
00219 {
00220   throw std::runtime_error("Ifpack2::SparsityFilter does not implement getNumEntriesInGlobalRow.");
00221 }
00222 
00223 //==========================================================================  
00224 template<class MatrixType> 
00225 size_t SparsityFilter<MatrixType>::getNumEntriesInLocalRow(LocalOrdinal localRow) const
00226 {
00227   return NumEntries_[localRow];
00228 }
00229 
00230 //==========================================================================  
00231 template<class MatrixType>   
00232 global_size_t SparsityFilter<MatrixType>::getGlobalNumDiags() const
00233 {
00234   return A_->getGlobalNumDiags();
00235 }
00236 
00237 //==========================================================================  
00238 template<class MatrixType>   
00239 size_t SparsityFilter<MatrixType>::getNodeNumDiags() const
00240 {
00241   return A_->getNodeNumDiags();
00242 }
00243 
00244 //==========================================================================  
00245 template<class MatrixType>   
00246 size_t SparsityFilter<MatrixType>::getGlobalMaxNumRowEntries() const
00247 {
00248   return MaxNumEntries_;
00249 }
00250 
00251 //==========================================================================  
00252 template<class MatrixType>   
00253 size_t SparsityFilter<MatrixType>::getNodeMaxNumRowEntries() const
00254 {
00255   return MaxNumEntries_;
00256 }
00257 
00258 //==========================================================================  
00259 template<class MatrixType>   
00260 bool SparsityFilter<MatrixType>::hasColMap() const
00261 {
00262   return true;
00263 }
00264 
00265 //==========================================================================  
00266 template<class MatrixType>   
00267 bool SparsityFilter<MatrixType>::isLowerTriangular() const
00268 {
00269   return A_->isLowerTriangular();
00270 }
00271 
00272 //==========================================================================  
00273 template<class MatrixType>   
00274 bool SparsityFilter<MatrixType>::isUpperTriangular() const
00275 {
00276   return A_->isUpperTriangular();
00277 }
00278 
00279 //==========================================================================  
00280 template<class MatrixType>   
00281 bool SparsityFilter<MatrixType>::isLocallyIndexed() const
00282 {
00283   return A_->isLocallyIndexed();
00284 }
00285 
00286 //==========================================================================  
00287 template<class MatrixType>   
00288 bool SparsityFilter<MatrixType>::isGloballyIndexed() const
00289 {
00290   return A_->isGloballyIndexed();
00291 }
00292 
00293 //==========================================================================  
00294 template<class MatrixType>   
00295 bool SparsityFilter<MatrixType>::isFillComplete() const
00296 {
00297   return A_->isFillComplete();
00298 }
00299   
00300 //==========================================================================
00301 template<class MatrixType> 
00302 void SparsityFilter<MatrixType>::getGlobalRowCopy(GlobalOrdinal GlobalRow,
00303               const Teuchos::ArrayView<GlobalOrdinal> &Indices,
00304               const Teuchos::ArrayView<Scalar> &Values,
00305               size_t &NumEntries) const
00306 {
00307   throw std::runtime_error("Ifpack2::SparsityFilter does not implement getGlobalRowCopy.");
00308 }
00309 
00310 //==========================================================================  
00311 template<class MatrixType> 
00312 void SparsityFilter<MatrixType>::getLocalRowCopy(LocalOrdinal LocalRow, 
00313                 const Teuchos::ArrayView<LocalOrdinal> &Indices, 
00314                 const Teuchos::ArrayView<Scalar> &Values,
00315                 size_t &NumEntries) const 
00316 { 
00317   TEUCHOS_TEST_FOR_EXCEPTION((LocalRow < 0 || (size_t) LocalRow >=  NumRows_ || (size_t) Indices.size() <  NumEntries_[LocalRow]), std::runtime_error, "Ifpack2::SparsityFilter::getLocalRowCopy invalid row or array size.");
00318 
00319   // Note: This function will work correctly if called by apply, say, with Indices_ and Values_ as
00320   // parameters.  The structure of the loop below should make that obvious.
00321 
00322   // always extract using the object Values_ and Indices_.
00323   // This is because I need more space than that given by
00324   // the user (for the external nodes)
00325   size_t A_NumEntries=0;
00326   A_->getLocalRowCopy(LocalRow,Indices_(),Values_(),A_NumEntries);
00327   magnitudeType Threshold = Teuchos::ScalarTraits<magnitudeType>::zero();
00328   std::vector<magnitudeType> Values2(A_NumEntries,Teuchos::ScalarTraits<magnitudeType>::zero());
00329 
00330   if (A_NumEntries > AllowedNumEntries_) {
00331     size_t count = 0;
00332     for (size_t i = 0 ; i < A_NumEntries ; ++i) {
00333       // skip diagonal entry (which is always inserted)
00334       if (Indices_[i] == LocalRow)
00335   continue;
00336       // put magnitude
00337       Values2[count] = Teuchos::ScalarTraits<Scalar>::magnitude(Values_[i]);
00338       count++;
00339     }
00340 
00341     // sort in descending order
00342     std::sort(Values2.rbegin(),Values2.rend());
00343     // get the cut-off value
00344     Threshold = Values2[AllowedNumEntries_ - 1];
00345   }
00346   
00347 
00348   // loop over all nonzero elements of row MyRow,
00349   // and drop elements below specified threshold.
00350   // Also, add value to zero diagonal
00351   NumEntries = 0;
00352   for (size_t i = 0 ; i < A_NumEntries ; ++i) {
00353     if (std::abs(Indices_[i] - LocalRow) > AllowedBandwidth_)
00354       continue;
00355     
00356     if ((Indices_[i] != LocalRow) && (Teuchos::ScalarTraits<Scalar>::magnitude(Values_[i]) < Threshold))
00357       continue;
00358 
00359     Values[NumEntries] = Values_[i];
00360     Indices[NumEntries] = Indices_[i];
00361 
00362     NumEntries++;
00363     if (NumEntries > AllowedNumEntries_)
00364       break;
00365   }
00366 
00367 }
00368 
00369 //==========================================================================  
00370 template<class MatrixType> 
00371 void SparsityFilter<MatrixType>::getGlobalRowView(GlobalOrdinal GlobalRow, 
00372               Teuchos::ArrayView<const GlobalOrdinal> &indices, 
00373               Teuchos::ArrayView<const Scalar> &values) const
00374 {
00375   throw std::runtime_error("Ifpack2::SparsityFilter: does not support getGlobalRowView.");
00376 }
00377 
00378 //==========================================================================  
00379 template<class MatrixType> 
00380 void SparsityFilter<MatrixType>::getLocalRowView(LocalOrdinal LocalRow, 
00381              Teuchos::ArrayView<const LocalOrdinal> &indices, 
00382              Teuchos::ArrayView<const Scalar> &values) const
00383 {
00384   throw std::runtime_error("Ifpack2::SparsityFilter: does not support getLocalRowView.");
00385 }
00386 
00387 //==========================================================================  
00388 template<class MatrixType> 
00389 void SparsityFilter<MatrixType>::getLocalDiagCopy(Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const
00390 {
00391   // This is somewhat dubious as to how the maps match.
00392   return A_->getLocalDiagCopy(diag);
00393 }
00394 
00395 //========================================================================== 
00396 template<class MatrixType> 
00397 void SparsityFilter<MatrixType>::leftScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x) 
00398 {
00399   throw std::runtime_error("Ifpack2::SparsityFilter does not support leftScale.");
00400 }
00401 
00402 //==========================================================================  
00403 template<class MatrixType> 
00404 void SparsityFilter<MatrixType>::rightScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x) 
00405 {
00406   throw std::runtime_error("Ifpack2::SparsityFilter does not support rightScale.");
00407 }
00408 
00409 //==========================================================================  
00410 template<class MatrixType> 
00411 void SparsityFilter<MatrixType>::apply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X, 
00412                Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 
00413                Teuchos::ETransp mode, 
00414                Scalar alpha,
00415                Scalar beta) const
00416 {  
00417   // Note: This isn't AztecOO compliant.  But neither was Ifpack's version.
00418   // Note: The localized maps mean the matvec is trivial (and has no import)
00419   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00420      "Ifpack2::SparsityFilter::apply ERROR: X.getNumVectors() != Y.getNumVectors().");
00421  
00422   Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00423   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > x_ptr = X.get2dView();
00424   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >       y_ptr = Y.get2dViewNonConst();
00425 
00426   Y.putScalar(zero);
00427   size_t NumVectors = Y.getNumVectors();
00428 
00429   for (size_t i = 0 ; i < NumRows_ ; ++i) {
00430     size_t Nnz;
00431     // Use this class's getrow to make the below code simpler
00432     getLocalRowCopy(i,Indices_(),Values_(),Nnz);
00433     if (mode==Teuchos::NO_TRANS){
00434       for (size_t j = 0 ; j < Nnz ; ++j) 
00435   for (size_t k = 0 ; k < NumVectors ; ++k)
00436     y_ptr[k][i] += Values_[j] * x_ptr[k][Indices_[j]];      
00437     }
00438     else if (mode==Teuchos::TRANS){
00439       for (size_t j = 0 ; j < Nnz ; ++j) 
00440   for (size_t k = 0 ; k < NumVectors ; ++k)
00441     y_ptr[k][Indices_[j]] += Values_[j] * x_ptr[k][i];
00442     }
00443     else { //mode==Teuchos::CONJ_TRANS
00444       for (size_t j = 0 ; j < Nnz ; ++j) 
00445   for (size_t k = 0 ; k < NumVectors ; ++k)
00446     y_ptr[k][Indices_[j]] += Teuchos::ScalarTraits<Scalar>::conjugate(Values_[j]) * x_ptr[k][i];
00447     }
00448   }
00449 }
00450   
00451 
00452 //==========================================================================  
00453 template<class MatrixType> 
00454 bool SparsityFilter<MatrixType>::hasTransposeApply() const
00455 {
00456   return true;
00457 }
00458 
00459 //==========================================================================
00460 template<class MatrixType> 
00461 bool SparsityFilter<MatrixType>::supportsRowViews() const
00462 {
00463   return false;
00464 }
00465 
00466 //==========================================================================  
00467 template<class MatrixType> 
00468 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType SparsityFilter<MatrixType>::getFrobeniusNorm() const
00469 {
00470   throw std::runtime_error("Ifpack2::SparsityFilter does not implement getFrobeniusNorm.");
00471 }
00472 
00473 //==========================================================================  
00474 template<class MatrixType> 
00475 TPETRA_DEPRECATED  void SparsityFilter<MatrixType>::getGlobalRowView(GlobalOrdinal GlobalRow, 
00476                      Teuchos::ArrayRCP<const GlobalOrdinal> &indices,
00477                      Teuchos::ArrayRCP<const Scalar>        &values) const
00478 {
00479   throw std::runtime_error("Ifpack2::SparsityFilter does not implement getGlobalRowView.");
00480 }
00481 
00482 //==========================================================================  
00483 template<class MatrixType> 
00484 TPETRA_DEPRECATED  void SparsityFilter<MatrixType>::getLocalRowView(LocalOrdinal LocalRow,
00485                     Teuchos::ArrayRCP<const LocalOrdinal> &indices,
00486                     Teuchos::ArrayRCP<const Scalar>       &values) const
00487 {
00488   throw std::runtime_error("Ifpack2::SparsityFilter does not implement getLocalRowView.");
00489 }
00490 
00491 }// namespace Ifpack2
00492 
00493 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends