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