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