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