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