Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_SingletonFilter_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_SINGLETONFILTER_DEF_HPP
00044 #define IFPACK2_SINGLETONFILTER_DEF_HPP
00045 #include "Ifpack2_SingletonFilter_decl.hpp"
00046 
00047 #include "Tpetra_ConfigDefs.hpp"
00048 #include "Tpetra_RowMatrix.hpp"
00049 #include "Tpetra_Map.hpp"
00050 #include "Tpetra_MultiVector.hpp"
00051 #include "Tpetra_Vector.hpp"
00052 
00053 namespace Ifpack2 {
00054 
00055 //==========================================================================
00056 template<class MatrixType>
00057 SingletonFilter<MatrixType>::SingletonFilter(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Matrix):
00058   A_(Matrix),
00059   NumSingletons_(0),
00060   NumRows_(0),
00061   NumNonzeros_(0),
00062   MaxNumEntries_(0),
00063   MaxNumEntriesA_(0)
00064 {
00065 
00066   // use this filter only on serial matrices
00067   if (A_->getComm()->getSize() != 1 || A_->getNodeNumRows() != A_->getGlobalNumRows()) {
00068     throw std::runtime_error("Ifpack2::SingeltonFilter 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.");
00069   }
00070 
00071   // Number of rows in A
00072   size_t NumRowsA_ = A_->getNodeNumRows();
00073 
00074   // tentative value for MaxNumEntries. This is the number of
00075   // nonzeros in the local matrix
00076   MaxNumEntriesA_ = A_->getNodeMaxNumRowEntries();
00077 
00078   // ExtractMyRowCopy() will use these vectors
00079   Indices_.resize(MaxNumEntriesA_);
00080   Values_.resize(MaxNumEntriesA_);
00081 
00082   // Initialize reordering vector to -1
00083   Reorder_.resize(NumRowsA_);
00084   Reorder_.assign(Reorder_.size(),-1);
00085 
00086   // first check how may singletons I do have
00087   NumRows_=0;
00088   for (size_t i = 0 ; i < NumRowsA_ ; ++i) {
00089     size_t Nnz;
00090     A_->getLocalRowCopy(i,Indices_,Values_,Nnz);
00091     if (Nnz != 1) {
00092       Reorder_[i] = NumRows_++;
00093     }
00094     else {
00095       NumSingletons_++;
00096     }
00097   }
00098 
00099   // build the inverse reordering
00100   InvReorder_.resize(NumRows_);
00101   for (size_t i = 0 ; i < NumRowsA_ ; ++i) {
00102     if (Reorder_[i] < 0)
00103       continue;
00104     InvReorder_[Reorder_[i]] = i;
00105   }
00106   NumEntries_.resize(NumRows_);
00107   SingletonIndex_.resize(NumSingletons_);
00108 
00109 
00110   // now compute the nonzeros per row
00111   size_t count = 0;
00112   for (size_t i = 0 ; i < NumRowsA_ ; ++i) {
00113     size_t Nnz;
00114     A_->getLocalRowCopy(i,Indices_,Values_,Nnz);
00115     LocalOrdinal ii = Reorder_[i];
00116     if (ii >= 0) {
00117       NumEntries_[ii] = Nnz;
00118       NumNonzeros_ += Nnz;
00119       if (Nnz > MaxNumEntries_)
00120         MaxNumEntries_ = Nnz;
00121     }
00122     else {
00123       SingletonIndex_[count] = i;
00124       count++;
00125     }
00126   }
00127 
00128   // Build the reduced map.  This map should be serial
00129   ReducedMap_ = Teuchos::rcp( new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>(NumRows_,0,A_->getComm()) );
00130 
00131   // and finish up with the diagonal entry
00132   Diagonal_ = Teuchos::rcp( new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(ReducedMap_) );
00133 
00134   Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> DiagonalA(A_->getRowMap());
00135   A_->getLocalDiagCopy(DiagonalA);
00136   const Teuchos::ArrayRCP<const Scalar> & DiagonalAview = DiagonalA.get1dView();
00137   for (size_t i = 0 ; i < NumRows_ ; ++i) {
00138     LocalOrdinal ii = InvReorder_[i];
00139     Diagonal_->replaceLocalValue((LocalOrdinal)i,DiagonalAview[ii]);
00140   }
00141 }
00142 
00143 //=========================================================================
00144 template<class MatrixType>
00145 SingletonFilter<MatrixType>::~SingletonFilter() { }
00146 
00147 //==========================================================================
00148 template<class MatrixType>
00149 Teuchos::RCP<const Teuchos::Comm<int> >
00150 SingletonFilter<MatrixType>::getComm() const
00151 {
00152   return A_->getComm();
00153 }
00154 
00155 //==========================================================================
00156 template<class MatrixType>
00157 Teuchos::RCP<typename MatrixType::node_type>
00158 SingletonFilter<MatrixType>::getNode() const
00159 {
00160   return A_->getNode();
00161 }
00162 
00163 //==========================================================================
00164 template<class MatrixType>
00165 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00166                                typename MatrixType::global_ordinal_type,
00167                                typename MatrixType::node_type> >
00168 SingletonFilter<MatrixType>::getRowMap() const
00169 {
00170   return ReducedMap_;
00171 }
00172 
00173 //==========================================================================
00174 template<class MatrixType>
00175 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00176                                typename MatrixType::global_ordinal_type,
00177                                typename MatrixType::node_type> >
00178 SingletonFilter<MatrixType>::getColMap() const
00179 {
00180   return ReducedMap_;
00181 }
00182 
00183 //==========================================================================
00184 template<class MatrixType>
00185 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00186                                typename MatrixType::global_ordinal_type,
00187                                typename MatrixType::node_type> >
00188 SingletonFilter<MatrixType>::getDomainMap() const
00189 {
00190   return ReducedMap_;
00191 }
00192 
00193 //==========================================================================
00194 template<class MatrixType>
00195 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00196                                typename MatrixType::global_ordinal_type,
00197                                typename MatrixType::node_type> >
00198 SingletonFilter<MatrixType>::getRangeMap() const
00199 {
00200   return ReducedMap_;
00201 }
00202 
00203 //==========================================================================
00204 template<class MatrixType>
00205 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type,
00206                                      typename MatrixType::global_ordinal_type,
00207                                      typename MatrixType::node_type> >
00208 SingletonFilter<MatrixType>::getGraph() const
00209 {
00210   throw std::runtime_error("Ifpack2::SingletonFilter: does not support getGraph.");
00211 }
00212 
00213 //==========================================================================
00214 template<class MatrixType>
00215 global_size_t SingletonFilter<MatrixType>::getGlobalNumRows() const
00216 {
00217   return NumRows_;
00218 }
00219 
00220 //==========================================================================
00221 template<class MatrixType>
00222 global_size_t SingletonFilter<MatrixType>::getGlobalNumCols() const
00223 {
00224   return NumRows_;
00225 }
00226 
00227 //==========================================================================
00228 template<class MatrixType>
00229 size_t SingletonFilter<MatrixType>::getNodeNumRows() const
00230 {
00231   return NumRows_;
00232 }
00233 
00234 //==========================================================================
00235 
00236 template<class MatrixType>
00237 size_t SingletonFilter<MatrixType>::getNodeNumCols() const
00238 {
00239   return NumRows_;
00240 }
00241 
00242 //==========================================================================
00243 template<class MatrixType>
00244 typename MatrixType::global_ordinal_type SingletonFilter<MatrixType>::getIndexBase() const
00245 {
00246   return A_->getIndexBase();
00247 }
00248 
00249 //==========================================================================
00250 template<class MatrixType>
00251 global_size_t SingletonFilter<MatrixType>::getGlobalNumEntries() const
00252 {
00253   return NumNonzeros_;
00254 }
00255 
00256 //==========================================================================
00257 template<class MatrixType>
00258 size_t SingletonFilter<MatrixType>::getNodeNumEntries() const
00259 {
00260   return NumNonzeros_;
00261 }
00262 
00263 //==========================================================================
00264 template<class MatrixType>
00265 size_t SingletonFilter<MatrixType>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
00266 {
00267   throw std::runtime_error("Ifpack2::SingletonFilter does not implement getNumEntriesInGlobalRow.");
00268 }
00269 
00270 //==========================================================================
00271 template<class MatrixType>
00272 size_t SingletonFilter<MatrixType>::getNumEntriesInLocalRow(LocalOrdinal localRow) const
00273 {
00274   return NumEntries_[localRow];
00275 }
00276 
00277 //==========================================================================
00278 template<class MatrixType>
00279 global_size_t SingletonFilter<MatrixType>::getGlobalNumDiags() const
00280 {
00281   return A_->getGlobalNumDiags();
00282 }
00283 
00284 //==========================================================================
00285 template<class MatrixType>
00286 size_t SingletonFilter<MatrixType>::getNodeNumDiags() const
00287 {
00288   return A_->getNodeNumDiags();
00289 }
00290 
00291 //==========================================================================
00292 template<class MatrixType>
00293 size_t SingletonFilter<MatrixType>::getGlobalMaxNumRowEntries() const
00294 {
00295   return MaxNumEntries_;
00296 }
00297 
00298 //==========================================================================
00299 template<class MatrixType>
00300 size_t SingletonFilter<MatrixType>::getNodeMaxNumRowEntries() const
00301 {
00302   return MaxNumEntries_;
00303 }
00304 
00305 //==========================================================================
00306 template<class MatrixType>
00307 bool SingletonFilter<MatrixType>::hasColMap() const
00308 {
00309   return true;
00310 }
00311 
00312 //==========================================================================
00313 template<class MatrixType>
00314 bool SingletonFilter<MatrixType>::isLowerTriangular() const
00315 {
00316   return A_->isLowerTriangular();
00317 }
00318 
00319 //==========================================================================
00320 template<class MatrixType>
00321 bool SingletonFilter<MatrixType>::isUpperTriangular() const
00322 {
00323   return A_->isUpperTriangular();
00324 }
00325 
00326 //==========================================================================
00327 template<class MatrixType>
00328 bool SingletonFilter<MatrixType>::isLocallyIndexed() const
00329 {
00330   return A_->isLocallyIndexed();
00331 }
00332 
00333 //==========================================================================
00334 template<class MatrixType>
00335 bool SingletonFilter<MatrixType>::isGloballyIndexed() const
00336 {
00337   return A_->isGloballyIndexed();
00338 }
00339 
00340 //==========================================================================
00341 template<class MatrixType>
00342 bool SingletonFilter<MatrixType>::isFillComplete() const
00343 {
00344   return A_->isFillComplete();
00345 }
00346 
00347 //==========================================================================
00348 template<class MatrixType>
00349 void SingletonFilter<MatrixType>::getGlobalRowCopy(GlobalOrdinal GlobalRow,
00350                                                   const Teuchos::ArrayView<GlobalOrdinal> &Indices,
00351                                                   const Teuchos::ArrayView<Scalar> &Values,
00352                                                   size_t &NumEntries) const
00353 {
00354   throw std::runtime_error("Ifpack2::SingletonFilter does not implement getGlobalRowCopy.");
00355 }
00356 
00357 //==========================================================================
00358 template<class MatrixType>
00359 void SingletonFilter<MatrixType>::getLocalRowCopy(LocalOrdinal LocalRow,
00360                                               const Teuchos::ArrayView<LocalOrdinal> &Indices,
00361                                               const Teuchos::ArrayView<Scalar> &Values,
00362                                               size_t &NumEntries) const
00363 {
00364   TEUCHOS_TEST_FOR_EXCEPTION((LocalRow < 0 || (size_t) LocalRow >=  NumRows_ || (size_t) Indices.size() <  NumEntries_[LocalRow]), std::runtime_error, "Ifpack2::SingletonFilter::getLocalRowCopy invalid row or array size.");
00365 
00366   size_t Nnz;
00367   LocalOrdinal ARow = InvReorder_[LocalRow];
00368   A_->getLocalRowCopy(ARow,Indices_(),Values_(),Nnz);
00369 
00370   // populate the user's vectors
00371   NumEntries = 0;
00372   for (size_t i = 0 ; i < Nnz ; ++i) {
00373     LocalOrdinal ii = Reorder_[Indices_[i]];
00374     if ( ii >= 0) {
00375       Indices[NumEntries] = ii;
00376       Values[NumEntries] = Values_[i];
00377       NumEntries++;
00378     }
00379   }
00380 
00381 }
00382 
00383 //==========================================================================
00384 template<class MatrixType>
00385 void SingletonFilter<MatrixType>::getGlobalRowView(GlobalOrdinal GlobalRow,
00386                                                   Teuchos::ArrayView<const GlobalOrdinal> &indices,
00387                                                   Teuchos::ArrayView<const Scalar> &values) const
00388 {
00389   throw std::runtime_error("Ifpack2::SingletonFilter: does not support getGlobalRowView.");
00390 }
00391 
00392 //==========================================================================
00393 template<class MatrixType>
00394 void SingletonFilter<MatrixType>::getLocalRowView(LocalOrdinal LocalRow,
00395                                                  Teuchos::ArrayView<const LocalOrdinal> &indices,
00396                                                  Teuchos::ArrayView<const Scalar> &values) const
00397 {
00398   throw std::runtime_error("Ifpack2::SingletonFilter: does not support getLocalRowView.");
00399 }
00400 
00401 //==========================================================================
00402 template<class MatrixType>
00403 void SingletonFilter<MatrixType>::getLocalDiagCopy(Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const
00404 {
00405   // This is somewhat dubious as to how the maps match.
00406   return A_->getLocalDiagCopy(diag);
00407 }
00408 
00409 //==========================================================================
00410 template<class MatrixType>
00411 void SingletonFilter<MatrixType>::leftScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x)
00412 {
00413   throw std::runtime_error("Ifpack2::SingletonFilter does not support leftScale.");
00414 }
00415 
00416 //==========================================================================
00417 template<class MatrixType>
00418 void SingletonFilter<MatrixType>::rightScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x)
00419 {
00420   throw std::runtime_error("Ifpack2::SingletonFilter does not support rightScale.");
00421 }
00422 
00423 //==========================================================================
00424 template<class MatrixType>
00425 void SingletonFilter<MatrixType>::apply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
00426                                        Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
00427                                        Teuchos::ETransp mode,
00428                                        Scalar alpha,
00429                                        Scalar beta) const
00430 {
00431   typedef Scalar DomainScalar;
00432   typedef Scalar RangeScalar;
00433 
00434   // Note: This isn't AztecOO compliant.  But neither was Ifpack's version.
00435 
00436   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00437                              "Ifpack2::SingletonFilter::apply ERROR: X.getNumVectors() != Y.getNumVectors().");
00438 
00439   RangeScalar zero = Teuchos::ScalarTraits<RangeScalar>::zero();
00440   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const DomainScalar> > x_ptr = X.get2dView();
00441   Teuchos::ArrayRCP<Teuchos::ArrayRCP<RangeScalar> >        y_ptr = Y.get2dViewNonConst();
00442 
00443   Y.putScalar(zero);
00444   size_t NumVectors = Y.getNumVectors();
00445 
00446 
00447   for (size_t i = 0 ; i < NumRows_ ; ++i) {
00448     size_t Nnz;
00449     // Use this class's getrow to make the below code simpler
00450     getLocalRowCopy(i,Indices_(),Values_(),Nnz);
00451     if (mode==Teuchos::NO_TRANS){
00452       for (size_t j = 0 ; j < Nnz ; ++j)
00453         for (size_t k = 0 ; k < NumVectors ; ++k)
00454           y_ptr[k][i] += (RangeScalar)Values_[j] * (RangeScalar)x_ptr[k][Indices_[j]];
00455     }
00456     else if (mode==Teuchos::TRANS){
00457       for (size_t j = 0 ; j < Nnz ; ++j)
00458         for (size_t k = 0 ; k < NumVectors ; ++k)
00459           y_ptr[k][Indices_[j]] += (RangeScalar)Values_[j] * (RangeScalar)x_ptr[k][i];
00460     }
00461     else { //mode==Teuchos::CONJ_TRANS
00462       for (size_t j = 0 ; j < Nnz ; ++j)
00463         for (size_t k = 0 ; k < NumVectors ; ++k)
00464           y_ptr[k][Indices_[j]] += Teuchos::ScalarTraits<RangeScalar>::conjugate((RangeScalar)Values_[j]) * (RangeScalar)x_ptr[k][i];
00465     }
00466   }
00467 }
00468 
00469 
00470 //==========================================================================
00471 template<class MatrixType>
00472 bool SingletonFilter<MatrixType>::hasTransposeApply() const
00473 {
00474   return true;
00475 }
00476 
00477 //==========================================================================
00478 template<class MatrixType>
00479 bool SingletonFilter<MatrixType>::supportsRowViews() const
00480 {
00481   return false;
00482 }
00483 
00484 //==============================================================================
00485 template<class MatrixType>
00486 void SingletonFilter<MatrixType>::SolveSingletons(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& RHS,
00487                                                   Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& LHS)
00488 {
00489   this->template SolveSingletonsTempl<Scalar,Scalar>(RHS, LHS);
00490 }
00491 
00492 //==============================================================================
00493 template<class MatrixType>
00494 template<class DomainScalar, class RangeScalar>
00495 void SingletonFilter<MatrixType>::SolveSingletonsTempl(const Tpetra::MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& RHS,
00496                                                        Tpetra::MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& LHS)
00497 {
00498   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const DomainScalar> > RHS_ptr = RHS.get2dView();
00499   Teuchos::ArrayRCP<Teuchos::ArrayRCP<RangeScalar> >        LHS_ptr = LHS.get2dViewNonConst();
00500 
00501   for (size_t i = 0 ; i < NumSingletons_ ; ++i) {
00502     LocalOrdinal ii = SingletonIndex_[i];
00503     // get the diagonal value for the singleton
00504     size_t Nnz;
00505     A_->getLocalRowCopy(ii,Indices_(),Values_(),Nnz);
00506     for (size_t j = 0 ; j < Nnz ; ++j) {
00507       if (Indices_[j] == ii) {
00508         for (size_t k = 0 ; k < LHS.getNumVectors() ; ++k)
00509           LHS_ptr[k][ii] = (RangeScalar)RHS_ptr[k][ii] / (RangeScalar)Values_[j];
00510       }
00511     }
00512   }
00513 }
00514 
00515 //==============================================================================
00516 template<class MatrixType>
00517 void SingletonFilter<MatrixType>::CreateReducedRHS(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& LHS,
00518                                                    const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& RHS,
00519                                                    Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& ReducedRHS)
00520 {
00521   this->template CreateReducedRHSTempl<Scalar,Scalar>(LHS, RHS, ReducedRHS);
00522 }
00523 
00524 //==============================================================================
00525 template<class MatrixType>
00526 template<class DomainScalar, class RangeScalar>
00527 void SingletonFilter<MatrixType>::CreateReducedRHSTempl(const Tpetra::MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& LHS,
00528                                                         const Tpetra::MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& RHS,
00529                                                         Tpetra::MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& ReducedRHS)
00530 {
00531   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const RangeScalar > > RHS_ptr = RHS.get2dView();
00532   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const DomainScalar> > LHS_ptr = LHS.get2dView();
00533   Teuchos::ArrayRCP<Teuchos::ArrayRCP<RangeScalar> >        ReducedRHS_ptr = ReducedRHS.get2dViewNonConst();
00534 
00535   size_t NumVectors = LHS.getNumVectors();
00536 
00537   for (size_t i = 0 ; i < NumRows_ ; ++i)
00538     for (size_t k = 0 ; k < NumVectors ; ++k)
00539       ReducedRHS_ptr[k][i] = RHS_ptr[k][InvReorder_[i]];
00540 
00541   for (size_t i = 0 ; i < NumRows_ ; ++i) {
00542     LocalOrdinal ii = InvReorder_[i];
00543     size_t Nnz;
00544     A_->getLocalRowCopy(ii,Indices_(),Values_(),Nnz);
00545 
00546     for (size_t j = 0 ; j < Nnz ; ++j) {
00547       if (Reorder_[Indices_[j]] == -1) {
00548         for (size_t k = 0 ; k < NumVectors ; ++k)
00549           ReducedRHS_ptr[k][i] -= (RangeScalar)Values_[j] * (RangeScalar)LHS_ptr[k][Indices_[j]];
00550       }
00551     }
00552   }
00553 }
00554 
00555 //==============================================================================
00556 template<class MatrixType>
00557 void SingletonFilter<MatrixType>::UpdateLHS(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& ReducedLHS,
00558                                             Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& LHS)
00559 {
00560   this->template UpdateLHSTempl<Scalar,Scalar>(ReducedLHS, LHS);
00561 }
00562 
00563 //==============================================================================
00564 template<class MatrixType>
00565 template<class DomainScalar, class RangeScalar>
00566 void SingletonFilter<MatrixType>::UpdateLHSTempl(const Tpetra::MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& ReducedLHS,
00567                                                  Tpetra::MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& LHS)
00568 {
00569 
00570   Teuchos::ArrayRCP<Teuchos::ArrayRCP<RangeScalar> >        LHS_ptr = LHS.get2dViewNonConst();
00571   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const DomainScalar> >  ReducedLHS_ptr = ReducedLHS.get2dView();
00572 
00573   for (size_t i = 0 ; i < NumRows_ ; ++i)
00574     for (size_t k = 0 ; k < LHS.getNumVectors() ; ++k)
00575       LHS_ptr[k][InvReorder_[i]] = (RangeScalar)ReducedLHS_ptr[k][i];
00576 }
00577 
00578 //==========================================================================
00579 template<class MatrixType>
00580 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType SingletonFilter<MatrixType>::getFrobeniusNorm() const
00581 {
00582   throw std::runtime_error("Ifpack2::SingletonFilter does not implement getFrobeniusNorm.");
00583 }
00584 
00585 //==========================================================================
00586 template<class MatrixType>
00587 TPETRA_DEPRECATED  void SingletonFilter<MatrixType>::getGlobalRowView(GlobalOrdinal GlobalRow,
00588                                                                      Teuchos::ArrayRCP<const GlobalOrdinal> &indices,
00589                                                                      Teuchos::ArrayRCP<const Scalar>        &values) const
00590 {
00591   throw std::runtime_error("Ifpack2::SingletonFilter does not implement getGlobalRowView.");
00592 }
00593 
00594 //==========================================================================
00595 template<class MatrixType>
00596 TPETRA_DEPRECATED  void SingletonFilter<MatrixType>::getLocalRowView(LocalOrdinal LocalRow,
00597                                                                     Teuchos::ArrayRCP<const LocalOrdinal> &indices,
00598                                                                     Teuchos::ArrayRCP<const Scalar>       &values) const
00599 {
00600   throw std::runtime_error("Ifpack2::SingletonFilter does not implement getLocalRowView.");
00601 }
00602 
00603 }// namespace Ifpack2
00604 
00605 #define IFPACK2_SINGLETONFILTER_INSTANT(S,LO,GO,N)                            \
00606   template class Ifpack2::SingletonFilter< Tpetra::CrsMatrix<S, LO, GO, N> >;
00607 
00608 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends