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