Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_DiagonalFilter_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_DIAGONALFILTER_DEF_HPP
00044 #define IFPACK2_DIAGONALFILTER_DEF_HPP
00045 #include "Ifpack2_DiagonalFilter_decl.hpp"
00046 #include <vector>
00047 
00048 #include "Teuchos_Comm.hpp"
00049 #include "Tpetra_ConfigDefs.hpp"
00050 #include "Tpetra_RowMatrix.hpp"
00051 #include "Tpetra_Map.hpp"
00052 #include "Tpetra_MultiVector.hpp"
00053 #include "Tpetra_Vector.hpp"
00054 
00055 
00056 namespace Ifpack2 {
00057 
00058 //==========================================================================
00059 template<class MatrixType>
00060 DiagonalFilter<MatrixType>::DiagonalFilter(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Matrix,
00061              typename Teuchos::ScalarTraits<Scalar>::magnitudeType AbsoluteThreshold,
00062              typename Teuchos::ScalarTraits<Scalar>::magnitudeType RelativeThreshold):
00063   A_(Matrix),  
00064   AbsoluteThreshold_(AbsoluteThreshold),
00065   RelativeThreshold_(RelativeThreshold)
00066 {
00067   pos_.resize(getNodeNumRows());
00068   val_=Teuchos::rcp(new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A_->getRowMap()));
00069   
00070   std::vector<LocalOrdinal> Indices(getNodeMaxNumRowEntries());
00071   std::vector<Scalar> Values(getNodeMaxNumRowEntries());
00072   size_t NumEntries;
00073   magnitudeType mysign;
00074  
00075 
00076   for (size_t MyRow = 0 ; MyRow < getNodeNumRows() ; ++MyRow) {    
00077     pos_[MyRow] = -1;
00078     A_->getLocalRowCopy(MyRow,Indices,Values,NumEntries);
00079 
00080     for (size_t i = 0 ; i < NumEntries ; ++i) {
00081       if ((size_t)Indices[i] == MyRow) {
00082   pos_[MyRow] = i;
00083   if(Teuchos::ScalarTraits<Scalar>::real(Values[i]) < 0)
00084     mysign=-Teuchos::ScalarTraits<magnitudeType>::one();
00085   else
00086     mysign=Teuchos::ScalarTraits<magnitudeType>::one();
00087 
00088 
00089   val_->replaceLocalValue(MyRow, Values[i] * (RelativeThreshold_ - 1) +
00090         AbsoluteThreshold_ * mysign);
00091   break;
00092       }
00093     }
00094   }
00095 }
00096 
00097 //=========================================================================
00098 template<class MatrixType>
00099 DiagonalFilter<MatrixType>::~DiagonalFilter() { }
00100 
00101 //==========================================================================
00102 template<class MatrixType>
00103 Teuchos::RCP<const Teuchos::Comm<int> > DiagonalFilter<MatrixType>::getComm() const
00104 {
00105   return A_->getComm();
00106 }
00107 
00108 //==========================================================================
00109 template<class MatrixType>
00110 Teuchos::RCP <typename MatrixType::node_type> DiagonalFilter<MatrixType>::getNode() const
00111 {
00112   return A_->getNode();
00113 }
00114 
00115 //==========================================================================
00116 template<class MatrixType>
00117 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00118              typename MatrixType::global_ordinal_type,
00119              typename MatrixType::node_type> >
00120 DiagonalFilter<MatrixType>::getRowMap() const
00121 {
00122   return A_->getRowMap();
00123 }
00124 
00125 //==========================================================================
00126 template<class MatrixType>
00127 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00128              typename MatrixType::global_ordinal_type,
00129              typename MatrixType::node_type> >
00130 DiagonalFilter<MatrixType>::getColMap() const
00131 {
00132   return A_->getColMap();
00133 }
00134 
00135 //==========================================================================
00136 template<class MatrixType>
00137 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00138              typename MatrixType::global_ordinal_type,
00139              typename MatrixType::node_type> >
00140 DiagonalFilter<MatrixType>::getDomainMap() const
00141 {
00142   return A_->getDomainMap();
00143 }
00144 
00145 //==========================================================================
00146 template<class MatrixType>
00147 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00148              typename MatrixType::global_ordinal_type,
00149              typename MatrixType::node_type> >
00150 DiagonalFilter<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 DiagonalFilter<MatrixType>::getGraph() const
00161 {
00162   return A_->getGraph();
00163 }
00164 
00165 //==========================================================================
00166 template<class MatrixType>
00167 global_size_t DiagonalFilter<MatrixType>::getGlobalNumRows() const
00168 {
00169   return A_->getGlobalNumRows();
00170 }
00171 
00172 //==========================================================================
00173 template<class MatrixType>
00174 global_size_t DiagonalFilter<MatrixType>::getGlobalNumCols() const
00175 {
00176   return A_->getGlobalNumCols();
00177 }
00178 
00179 //==========================================================================
00180 template<class MatrixType>
00181 size_t DiagonalFilter<MatrixType>::getNodeNumRows() const
00182 {
00183   return A_->getNodeNumRows();
00184 }
00185 
00186 //==========================================================================
00187  
00188 template<class MatrixType>
00189 size_t DiagonalFilter<MatrixType>::getNodeNumCols() const
00190 {
00191   return A_->getNodeNumCols();
00192 }
00193 
00194 //==========================================================================  
00195 template<class MatrixType>
00196 typename MatrixType::global_ordinal_type DiagonalFilter<MatrixType>::getIndexBase() const
00197 {
00198   return A_->getIndexBase();
00199 }
00200 
00201 //========================================================================== 
00202 template<class MatrixType>
00203 global_size_t DiagonalFilter<MatrixType>::getGlobalNumEntries() const
00204 {
00205   return A_->getGlobalNumEntries();
00206 }
00207 
00208 //========================================================================== 
00209 template<class MatrixType>
00210 size_t DiagonalFilter<MatrixType>::getNodeNumEntries() const
00211 {
00212   return A_->getNodeNumEntries();
00213 }
00214 
00215 //==========================================================================
00216 template<class MatrixType> 
00217 size_t DiagonalFilter<MatrixType>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
00218 {
00219   return A_->getNumEntriesInGlobalRow(globalRow);
00220 }
00221 
00222 //==========================================================================  
00223 template<class MatrixType> 
00224 size_t DiagonalFilter<MatrixType>::getNumEntriesInLocalRow(LocalOrdinal localRow) const
00225 {
00226   return A_->getNumEntriesInLocalRow(localRow);
00227 }
00228 
00229 //==========================================================================  
00230 template<class MatrixType>   
00231 global_size_t DiagonalFilter<MatrixType>::getGlobalNumDiags() const
00232 {
00233   return A_->getGlobalNumDiags();
00234 }
00235 
00236 //==========================================================================  
00237 template<class MatrixType>   
00238 size_t DiagonalFilter<MatrixType>::getNodeNumDiags() const
00239 {
00240   return A_->getNodeNumDiags();
00241 }
00242 
00243 //==========================================================================  
00244 template<class MatrixType>   
00245 size_t DiagonalFilter<MatrixType>::getGlobalMaxNumRowEntries() const
00246 {
00247   return A_->getGlobalMaxNumRowEntries();
00248 }
00249 
00250 //==========================================================================  
00251 template<class MatrixType>   
00252 size_t DiagonalFilter<MatrixType>::getNodeMaxNumRowEntries() const
00253 {
00254   return A_->getNodeMaxNumRowEntries();
00255 }
00256 
00257 //==========================================================================  
00258 template<class MatrixType>   
00259 bool DiagonalFilter<MatrixType>::hasColMap() const
00260 {
00261   return A_->hasColMap();
00262 }
00263 
00264 //==========================================================================  
00265 template<class MatrixType>   
00266 bool DiagonalFilter<MatrixType>::isLowerTriangular() const
00267 {
00268   return A_->isLowerTriangular();
00269 }
00270 
00271 //==========================================================================  
00272 template<class MatrixType>   
00273 bool DiagonalFilter<MatrixType>::isUpperTriangular() const
00274 {
00275   return A_->isUpperTriangular();
00276 }
00277 
00278 //==========================================================================  
00279 template<class MatrixType>   
00280 bool DiagonalFilter<MatrixType>::isLocallyIndexed() const
00281 {
00282   return A_->isLocallyIndexed();
00283 }
00284 
00285 //==========================================================================  
00286 template<class MatrixType>   
00287 bool DiagonalFilter<MatrixType>::isGloballyIndexed() const
00288 {
00289   return A_->isGloballyIndexed();
00290 }
00291 
00292 //==========================================================================  
00293 template<class MatrixType>   
00294 bool DiagonalFilter<MatrixType>::isFillComplete() const
00295 {
00296   return A_->isFillComplete();
00297 }
00298   
00299 //==========================================================================
00300 template<class MatrixType> 
00301 void DiagonalFilter<MatrixType>::getGlobalRowCopy(GlobalOrdinal GlobalRow,
00302               const Teuchos::ArrayView<GlobalOrdinal> &Indices,
00303               const Teuchos::ArrayView<Scalar> &Values,
00304               size_t &NumEntries) const
00305 {
00306   Teuchos::ArrayRCP< const Scalar > myvals=val_->get1dView();
00307   LocalOrdinal LocalRow=getRowMap()->getLocalElement(GlobalRow);
00308 
00309   A_->getGlobalRowCopy(GlobalRow, Indices,Values,NumEntries);
00310 
00311   if (pos_[LocalRow] != -1)
00312     Values[pos_[LocalRow]] += myvals[LocalRow];
00313 }
00314 
00315 //==========================================================================  
00316 template<class MatrixType> 
00317 void DiagonalFilter<MatrixType>::getLocalRowCopy(LocalOrdinal LocalRow, 
00318              const Teuchos::ArrayView<LocalOrdinal> &Indices, 
00319              const Teuchos::ArrayView<Scalar> &Values,
00320              size_t &NumEntries) const 
00321 { 
00322   Teuchos::ArrayRCP< const Scalar > myvals=val_->get1dView();
00323 
00324   A_->getLocalRowCopy(LocalRow, Indices,Values,NumEntries);
00325 
00326   if (pos_[LocalRow] != -1)
00327     Values[pos_[LocalRow]] += myvals[LocalRow];
00328 }
00329 
00330 //==========================================================================  
00331 template<class MatrixType> 
00332 void DiagonalFilter<MatrixType>::getGlobalRowView(GlobalOrdinal GlobalRow, 
00333               Teuchos::ArrayView<const GlobalOrdinal> &indices, 
00334               Teuchos::ArrayView<const Scalar> &values) const
00335 {
00336   throw std::runtime_error("Ifpack2::DiagonalFilter: does not support getGlobalRowView.");
00337 }
00338 
00339 //==========================================================================  
00340 template<class MatrixType> 
00341 void DiagonalFilter<MatrixType>::getLocalRowView(LocalOrdinal LocalRow, 
00342              Teuchos::ArrayView<const LocalOrdinal> &indices, 
00343              Teuchos::ArrayView<const Scalar> &values) const
00344 {
00345   throw std::runtime_error("Ifpack2::DiagonalFilter: does not support getLocalRowView.");
00346 }
00347 
00348 //==========================================================================  
00349 template<class MatrixType> 
00350 void DiagonalFilter<MatrixType>::getLocalDiagCopy(Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const
00351 {
00352   // Returns the matrix's actual diagonal, rather than the "filtered" diagonal.
00353   // This is dubious, but it duplicates the functionality of Old Ifpack.
00354   return A_->getLocalDiagCopy(diag);
00355 }
00356 
00357 //========================================================================== 
00358 template<class MatrixType> 
00359 void DiagonalFilter<MatrixType>::leftScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x) 
00360 {
00361   throw std::runtime_error("Ifpack2::DiagonalFilter does not support leftScale.");
00362 }
00363 
00364 //==========================================================================  
00365 template<class MatrixType> 
00366 void DiagonalFilter<MatrixType>::rightScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x) 
00367 {
00368   throw std::runtime_error("Ifpack2::DiagonalFilter does not support rightScale.");
00369 }
00370 
00371 //==========================================================================  
00372 template<class MatrixType> 
00373 void DiagonalFilter<MatrixType>::apply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X, 
00374                Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 
00375                Teuchos::ETransp mode, 
00376                Scalar alpha,
00377                Scalar beta) const
00378 {
00379   Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00380   A_->apply(X,Y,mode,alpha,beta);
00381   Y.elementWiseMultiply(one,*val_,X,one);
00382 }
00383   
00384 
00385 //==========================================================================  
00386 template<class MatrixType> 
00387 bool DiagonalFilter<MatrixType>::hasTransposeApply() const
00388 {
00389   return A_->hasTransposeApply();
00390 }
00391 
00392 //==========================================================================
00393 template<class MatrixType> 
00394 bool DiagonalFilter<MatrixType>::supportsRowViews() const
00395 {
00396   return false;
00397 }
00398 
00399 //==========================================================================  
00400 template<class MatrixType> 
00401 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType DiagonalFilter<MatrixType>::getFrobeniusNorm() const
00402 {
00403   throw std::runtime_error("Ifpack2::DiagonalFilter does not implement getFrobeniusNorm.");
00404 }
00405 
00406 //==========================================================================  
00407 template<class MatrixType> 
00408 TPETRA_DEPRECATED  void DiagonalFilter<MatrixType>::getGlobalRowView(GlobalOrdinal GlobalRow, 
00409                      Teuchos::ArrayRCP<const GlobalOrdinal> &indices,
00410                      Teuchos::ArrayRCP<const Scalar>        &values) const
00411 {
00412   throw std::runtime_error("Ifpack2::DiagonalFilter does not implement getGlobalRowView.");
00413 }
00414 
00415 //==========================================================================  
00416 template<class MatrixType> 
00417 TPETRA_DEPRECATED  void DiagonalFilter<MatrixType>::getLocalRowView(LocalOrdinal LocalRow,
00418                     Teuchos::ArrayRCP<const LocalOrdinal> &indices,
00419                     Teuchos::ArrayRCP<const Scalar>       &values) const
00420 {
00421   throw std::runtime_error("Ifpack2::DiagonalFilter does not implement getLocalRowView.");
00422 }
00423 
00424 }// namespace Ifpack2
00425 
00426 #define IFPACK2_DIAGONALFILTER_INSTANT(S,LO,GO,N)                            \
00427   template class Ifpack2::DiagonalFilter< Tpetra::CrsMatrix<S, LO, GO, N> >;
00428 
00429 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends