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