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 // 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_LOCALFILTER_DEF_HPP
00044 #define IFPACK2_LOCALFILTER_DEF_HPP
00045 
00046 #include <Ifpack2_LocalFilter_decl.hpp>
00047 #include <Tpetra_Map.hpp>
00048 #include <Tpetra_MultiVector.hpp>
00049 #include <Tpetra_Vector.hpp>
00050 
00051 #ifdef HAVE_MPI
00052 #  include "Teuchos_DefaultMpiComm.hpp"
00053 #else
00054 #  include "Teuchos_DefaultSerialComm.hpp"
00055 #endif
00056 
00057 namespace Ifpack2 {
00058 
00059 
00060 template<class MatrixType>
00061 bool
00062 LocalFilter<MatrixType>::
00063 mapPairsAreFitted (const row_matrix_type& A)
00064 {
00065   const map_type& rangeMap = * (A.getRangeMap ());
00066   const map_type& rowMap = * (A.getRowMap ());
00067   const bool rangeAndRowFitted = mapPairIsFitted (rangeMap, rowMap);
00068 
00069   const map_type& domainMap = * (A.getDomainMap ());
00070   const map_type& columnMap = * (A.getColMap ());
00071   const bool domainAndColumnFitted = mapPairIsFitted (domainMap, columnMap);
00072 
00073   return rangeAndRowFitted && domainAndColumnFitted;
00074 }
00075 
00076 
00077 template<class MatrixType>
00078 bool
00079 LocalFilter<MatrixType>::
00080 mapPairIsFitted (const map_type& map1, const map_type& map2)
00081 {
00082   using Teuchos::ArrayView;
00083   using Teuchos::as;
00084   typedef global_ordinal_type GO; // a handy abbreviation
00085   typedef typename ArrayView<const GO>::size_type size_type;
00086 
00087   bool fitted = true;
00088   if (&map1 == &map2) {
00089     fitted = true;
00090   }
00091   else if (map1.isContiguous () && map2.isContiguous () &&
00092            map1.getMinGlobalIndex () == map2.getMinGlobalIndex () &&
00093            map1.getMaxGlobalIndex () <= map2.getMaxGlobalIndex ()) {
00094     // Special case where both Maps are contiguous.
00095     fitted = true;
00096   }
00097   else {
00098     ArrayView<const GO> inds_map2 = map2.getNodeElementList ();
00099     const size_type numInds_map1 = as<size_type> (map1.getNodeNumElements ());
00100 
00101     if (map1.isContiguous ()) {
00102       // Avoid calling getNodeElementList() on the always one-to-one
00103       // Map, if it is contiguous (a common case).  When called on a
00104       // contiguous Map, getNodeElementList() causes allocation of an
00105       // array that sticks around, even though the array isn't needed.
00106       // (The Map is contiguous, so you can compute the entries; you
00107       // don't need to store them.)
00108       if (numInds_map1 > inds_map2.size ()) {
00109         // There are fewer indices in map1 on this process than in
00110         // map2.  This case might be impossible.
00111         fitted = false;
00112       }
00113       else {
00114         // Do all the map1 indices match the initial map2 indices?
00115         const GO minInd_map1 = map1.getMinGlobalIndex ();
00116         for (size_type k = 0; k < numInds_map1; ++k) {
00117           const GO inds_map1_k = as<GO> (k) + minInd_map1;
00118           if (inds_map1_k != inds_map2[k]) {
00119             fitted = false;
00120             break;
00121           }
00122         }
00123       }
00124     }
00125     else { // map1 is not contiguous.
00126       // Get index lists from both Maps, and compare their indices.
00127       ArrayView<const GO> inds_map1 = map1.getNodeElementList ();
00128       if (numInds_map1 > inds_map2.size ()) {
00129         // There are fewer indices in map1 on this process than in
00130         // map2.  This case might be impossible.
00131         fitted = false;
00132       }
00133       else {
00134         // Do all the map1 indices match those in map2?
00135         for (size_type k = 0; k < numInds_map1; ++k) {
00136           if (inds_map1[k] != inds_map2[k]) {
00137             fitted = false;
00138             break;
00139           }
00140         }
00141       }
00142     }
00143   }
00144   return fitted;
00145 }
00146 
00147 
00148 template<class MatrixType>
00149 LocalFilter<MatrixType>::
00150 LocalFilter (const Teuchos::RCP<const row_matrix_type>& A) :
00151   A_ (A),
00152   NumNonzeros_ (0),
00153   MaxNumEntries_ (0),
00154   MaxNumEntriesA_ (0)
00155 {
00156   using Teuchos::as;
00157   using Teuchos::RCP;
00158   using Teuchos::rcp;
00159 
00160 #ifdef HAVE_IFPACK2_DEBUG
00161   TEUCHOS_TEST_FOR_EXCEPTION(
00162     ! mapPairsAreFitted (*A), std::invalid_argument, "Ifpack2::LocalFilter: "
00163     "A's Map pairs are not fitted to each other on Process "
00164     << A_->getRowMap ()->getComm ()->getRank () << " of the input matrix's "
00165     "communicator.  "
00166     "This means that LocalFilter does not currently know how to work with A.  "
00167     "This will change soon.  Please see discussion of Bug 5992.");
00168 #endif // HAVE_IFPACK2_DEBUG
00169 
00170   // Build the local communicator (containing this process only).
00171   RCP<const Teuchos::Comm<int> > localComm;
00172 #ifdef HAVE_MPI
00173   localComm = rcp (new Teuchos::MpiComm<int> (MPI_COMM_SELF));
00174 #else
00175   localComm = rcp (new Teuchos::SerialComm<int> ());
00176 #endif // HAVE_MPI
00177 
00178 
00179   // // FIXME (mfh 21 Nov 2013) Currently, the implementation implicitly
00180   // // assumes that the range Map is fitted to the row Map.  Otherwise,
00181   // // it probably won't work at all.
00182   // TEUCHOS_TEST_FOR_EXCEPTION(
00183   //   mapPairIsFitted (* (A_->getRangeMap ()), * (A_->getRowMap ())),
00184   //   std::logic_error, "Ifpack2::LocalFilter: Range Map of the input matrix "
00185   //   "is not fitted to its row Map.  LocalFilter does not currently work in "
00186   //   "this case.  See Bug 5992.");
00187 
00188   // Build the local row, domain, and range Maps.  They both use the
00189   // local communicator built above.  The global indices of each are
00190   // different than those of the corresponding original Map; they
00191   // actually are the same as the local indices of the original Map.
00192   //
00193   // It's even OK if signedness of local_ordinal_type and
00194   // global_ordinal_type are different.  (That would be a BAD IDEA but
00195   // some users seem to enjoy making trouble for themselves and us.)
00196   //
00197   // Both the local row and local range Maps must have the same number
00198   // of entries, namely, that of the local number of entries of A's
00199   // range Map.
00200 
00201   // FIXME (mfh 21 Nov 2013) For some reason, we have to use this,
00202   // even if it differs from the number of entries in the range Map.
00203   // Otherwise, AdditiveSchwarz Test1 fails, down in the local solve,
00204   // where the matrix has 8 columns but the local part of the vector
00205   // only has five rows.
00206   const size_t numRows = A_->getNodeNumRows ();
00207 
00208   // using std::cerr;
00209   // using std::endl;
00210   // cerr << "A_ has " << A_->getNodeNumRows () << " rows." << endl
00211   //      << "Range Map has " << A_->getRangeMap ()->getNodeNumElements () << " entries." << endl
00212   //      << "Row Map has " << A_->getRowMap ()->getNodeNumElements () << " entries." << endl;
00213 
00214   const global_ordinal_type indexBase = as<global_ordinal_type> (0);
00215 
00216   localRowMap_ =
00217     rcp (new map_type (numRows, indexBase, localComm,
00218                        Tpetra::GloballyDistributed, A_->getNode ()));
00219   // If the original matrix's range Map is not fitted to its row Map,
00220   // we'll have to do an Export when applying the matrix.
00221   localRangeMap_ = localRowMap_;
00222 
00223   // If the original matrix's domain Map is not fitted to its column
00224   // Map, we'll have to do an Import when applying the matrix.
00225   const size_t numCols = A_->getDomainMap ()->getNodeNumElements ();
00226   if (A_->getRangeMap ().getRawPtr () == A_->getDomainMap ().getRawPtr ()) {
00227     // The input matrix's domain and range Maps are identical, so the
00228     // locally filtered matrix's domain and range Maps can be also.
00229     localDomainMap_ = localRangeMap_;
00230   }
00231   else {
00232     localDomainMap_ =
00233       rcp (new map_type (numCols, indexBase, localComm,
00234                          Tpetra::GloballyDistributed, A_->getNode ()));
00235   }
00236 
00237   // NodeNumEntries_ will contain the actual number of nonzeros for
00238   // each localized row (that is, without external nodes, and always
00239   // with the diagonal entry)
00240   NumEntries_.resize (numRows);
00241 
00242   // tentative value for MaxNumEntries. This is the number of
00243   // nonzeros in the local matrix
00244   MaxNumEntries_  = A_->getNodeMaxNumRowEntries ();
00245   MaxNumEntriesA_ = A_->getNodeMaxNumRowEntries ();
00246 
00247   // Allocate temporary arrays for getLocalRowCopy().
00248   localIndices_.resize (MaxNumEntries_);
00249   Values_.resize (MaxNumEntries_);
00250 
00251   // now compute:
00252   // - the number of nonzero per row
00253   // - the total number of nonzeros
00254   // - the diagonal entries
00255 
00256   // compute nonzeros (total and per-row), and store the
00257   // diagonal entries (already modified)
00258   size_t ActualMaxNumEntries = 0;
00259 
00260   for (size_t i = 0; i < numRows; ++i) {
00261     NumEntries_[i] = 0;
00262     size_t Nnz, NewNnz = 0;
00263     A_->getLocalRowCopy (i, localIndices_, Values_, Nnz);
00264     for (size_t j = 0; j < Nnz; ++j) {
00265       // FIXME (mfh 03 Apr 2013) This assumes the following:
00266       //
00267       // 1. Row Map, range Map, and domain Map are all the same.
00268       //
00269       // 2. The column Map's list of GIDs on this process is the
00270       //    domain Map's list of GIDs, followed by remote GIDs.  Thus,
00271       //    for any GID in the domain Map on this process, its LID in
00272       //    the domain Map (and therefore in the row Map, by (1)) is
00273       //    the same as its LID in the column Map.  (Hence the
00274       //    less-than test, which if true, means that localIndices_[j]
00275       //    belongs to the row Map.)
00276       if (Teuchos::as<size_t> (localIndices_[j]) < numRows) {
00277         ++NewNnz;
00278       }
00279     }
00280 
00281     if (NewNnz > ActualMaxNumEntries) {
00282       ActualMaxNumEntries = NewNnz;
00283     }
00284 
00285     NumNonzeros_ += NewNnz;
00286     NumEntries_[i] = NewNnz;
00287   }
00288 
00289   MaxNumEntries_ = ActualMaxNumEntries;
00290 }
00291 
00292 
00293 template<class MatrixType>
00294 LocalFilter<MatrixType>::~LocalFilter()
00295 {}
00296 
00297 
00298 template<class MatrixType>
00299 Teuchos::RCP<const Teuchos::Comm<int> >
00300 LocalFilter<MatrixType>::getComm () const
00301 {
00302   return localRowMap_->getComm ();
00303 }
00304 
00305 
00306 template<class MatrixType>
00307 Teuchos::RCP<typename MatrixType::node_type>
00308 LocalFilter<MatrixType>::getNode () const
00309 {
00310   return A_->getNode ();
00311 }
00312 
00313 
00314 template<class MatrixType>
00315 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00316                                typename MatrixType::global_ordinal_type,
00317                                typename MatrixType::node_type> >
00318 LocalFilter<MatrixType>::getRowMap () const
00319 {
00320   return localRowMap_;
00321 }
00322 
00323 
00324 template<class MatrixType>
00325 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00326                                typename MatrixType::global_ordinal_type,
00327                                typename MatrixType::node_type> >
00328 LocalFilter<MatrixType>::getColMap() const
00329 {
00330   return localRowMap_; // FIXME (mfh 20 Nov 2013)
00331 }
00332 
00333 
00334 template<class MatrixType>
00335 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00336                                typename MatrixType::global_ordinal_type,
00337                                typename MatrixType::node_type> >
00338 LocalFilter<MatrixType>::getDomainMap() const
00339 {
00340   return localDomainMap_;
00341 }
00342 
00343 
00344 template<class MatrixType>
00345 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00346                                typename MatrixType::global_ordinal_type,
00347                                typename MatrixType::node_type> >
00348 LocalFilter<MatrixType>::getRangeMap() const
00349 {
00350   return localRangeMap_;
00351 }
00352 
00353 
00354 template<class MatrixType>
00355 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type,
00356                                      typename MatrixType::global_ordinal_type,
00357                                      typename MatrixType::node_type> >
00358 LocalFilter<MatrixType>::getGraph () const
00359 {
00360   // FIXME (mfh 20 Nov 2013) This is not what the documentation says
00361   // this method should do!  It should return the graph of the locally
00362   // filtered matrix, not the original matrix's graph.
00363   return A_->getGraph ();
00364 }
00365 
00366 
00367 template<class MatrixType>
00368 global_size_t LocalFilter<MatrixType>::getGlobalNumRows() const
00369 {
00370   return Teuchos::as<global_size_t> (localRangeMap_->getNodeNumElements ());
00371 }
00372 
00373 
00374 template<class MatrixType>
00375 global_size_t LocalFilter<MatrixType>::getGlobalNumCols() const
00376 {
00377   return Teuchos::as<global_size_t> (localDomainMap_->getNodeNumElements ());
00378 }
00379 
00380 
00381 template<class MatrixType>
00382 size_t LocalFilter<MatrixType>::getNodeNumRows() const
00383 {
00384   return Teuchos::as<size_t> (localRangeMap_->getNodeNumElements ());
00385 }
00386 
00387 
00388 template<class MatrixType>
00389 size_t LocalFilter<MatrixType>::getNodeNumCols() const
00390 {
00391   return Teuchos::as<size_t> (localDomainMap_->getNodeNumElements ());
00392 }
00393 
00394 
00395 template<class MatrixType>
00396 typename MatrixType::global_ordinal_type
00397 LocalFilter<MatrixType>::getIndexBase () const
00398 {
00399   return A_->getIndexBase ();
00400 }
00401 
00402 
00403 template<class MatrixType>
00404 global_size_t LocalFilter<MatrixType>::getGlobalNumEntries () const
00405 {
00406   return NumNonzeros_;
00407 }
00408 
00409 
00410 template<class MatrixType>
00411 size_t LocalFilter<MatrixType>::getNodeNumEntries () const
00412 {
00413   return NumNonzeros_;
00414 }
00415 
00416 
00417 template<class MatrixType>
00418 size_t
00419 LocalFilter<MatrixType>::
00420 getNumEntriesInGlobalRow (global_ordinal_type globalRow) const
00421 {
00422   const local_ordinal_type localRow = getRowMap ()->getLocalElement (globalRow);
00423   if (localRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
00424     // NOTE (mfh 26 Mar 2014) We return zero if globalRow is not in
00425     // the row Map on this process, since "get the number of entries
00426     // in the global row" refers only to what the calling process owns
00427     // in that row.  In this case, it owns no entries in that row,
00428     // since it doesn't own the row.
00429     return 0;
00430   } else {
00431     return NumEntries_[localRow];
00432   }
00433 }
00434 
00435 
00436 template<class MatrixType>
00437 size_t
00438 LocalFilter<MatrixType>::
00439 getNumEntriesInLocalRow (local_ordinal_type localRow) const
00440 {
00441   if (getRowMap ()->isNodeLocalElement (localRow)) {
00442     return NumEntries_[localRow];
00443   } else {
00444     // NOTE (mfh 26 Mar 2014) We return zero if localRow is not in the
00445     // row Map on this process, since "get the number of entries in
00446     // the local row" refers only to what the calling process owns in
00447     // that row.  In this case, it owns no entries in that row, since
00448     // it doesn't own the row.
00449     return 0;
00450   }
00451 }
00452 
00453 
00454 template<class MatrixType>
00455 global_size_t LocalFilter<MatrixType>::getGlobalNumDiags () const
00456 {
00457   return A_->getGlobalNumDiags ();
00458 }
00459 
00460 
00461 template<class MatrixType>
00462 size_t LocalFilter<MatrixType>::getNodeNumDiags () const
00463 {
00464   return A_->getNodeNumDiags ();
00465 }
00466 
00467 
00468 template<class MatrixType>
00469 size_t LocalFilter<MatrixType>::getGlobalMaxNumRowEntries () const
00470 {
00471   return MaxNumEntries_;
00472 }
00473 
00474 
00475 template<class MatrixType>
00476 size_t LocalFilter<MatrixType>::getNodeMaxNumRowEntries() const
00477 {
00478   return MaxNumEntries_;
00479 }
00480 
00481 
00482 template<class MatrixType>
00483 bool LocalFilter<MatrixType>::hasColMap () const
00484 {
00485   return true;
00486 }
00487 
00488 
00489 template<class MatrixType>
00490 bool LocalFilter<MatrixType>::isLowerTriangular () const
00491 {
00492   return A_->isLowerTriangular();
00493 }
00494 
00495 
00496 template<class MatrixType>
00497 bool LocalFilter<MatrixType>::isUpperTriangular () const
00498 {
00499   return A_->isUpperTriangular();
00500 }
00501 
00502 
00503 template<class MatrixType>
00504 bool LocalFilter<MatrixType>::isLocallyIndexed () const
00505 {
00506   return A_->isLocallyIndexed ();
00507 }
00508 
00509 
00510 template<class MatrixType>
00511 bool LocalFilter<MatrixType>::isGloballyIndexed () const
00512 {
00513   return A_->isGloballyIndexed();
00514 }
00515 
00516 
00517 template<class MatrixType>
00518 bool LocalFilter<MatrixType>::isFillComplete () const
00519 {
00520   return A_->isFillComplete ();
00521 }
00522 
00523 
00524 template<class MatrixType>
00525 void
00526 LocalFilter<MatrixType>::
00527 getGlobalRowCopy (global_ordinal_type globalRow,
00528                   const Teuchos::ArrayView<global_ordinal_type>& globalIndices,
00529                   const Teuchos::ArrayView<scalar_type>& values,
00530                   size_t& numEntries) const
00531 {
00532   const local_ordinal_type localRow = getRowMap ()->getLocalElement (globalRow);
00533   if (localRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
00534     // NOTE (mfh 26 Mar 2014) We return no entries if globalRow is not
00535     // in the row Map on this process, since "get a copy of the
00536     // entries in the global row" refers only to what the calling
00537     // process owns in that row.  In this case, it owns no entries in
00538     // that row, since it doesn't own the row.
00539     numEntries = 0;
00540   }
00541   else {
00542     // First get a copy of the current row using local indices.  Then,
00543     // convert to global indices using the input matrix's column Map.
00544     //
00545     numEntries = getNumEntriesInLocalRow (localRow);
00546     // FIXME (mfh 26 Mar 2014) If local_ordinal_type ==
00547     // global_ordinal_type, we could just alias the input array
00548     // instead of allocating a temporary array.
00549     Teuchos::Array<local_ordinal_type> localIndices (numEntries);
00550     getLocalRowCopy (localRow, localIndices (), values, numEntries);
00551 
00552     const map_type& colMap = * (getColMap ());
00553 
00554     typedef typename Teuchos::Array<local_ordinal_type>::size_type size_type;
00555     // Don't fill the output array beyond its size.
00556     const size_type numEnt =
00557       std::min (static_cast<size_type> (numEntries), globalIndices.size ());
00558     for (size_type k = 0; k < numEnt; ++k) {
00559       globalIndices[k] = colMap.getGlobalElement (localIndices[k]);
00560     }
00561   }
00562 }
00563 
00564 
00565 template<class MatrixType>
00566 void
00567 LocalFilter<MatrixType>::
00568 getLocalRowCopy (local_ordinal_type LocalRow,
00569                  const Teuchos::ArrayView<local_ordinal_type> &Indices,
00570                  const Teuchos::ArrayView<scalar_type> &Values,
00571                  size_t &NumEntries) const
00572 {
00573   const size_t numRows = localRowMap_->getNodeNumElements ();
00574   const int myRank = localRowMap_->getComm ()->getSize ();
00575 
00576   if (LocalRow < localRowMap_->getMinLocalIndex () ||
00577         static_cast<size_t> (LocalRow) >= numRows) {
00578     std::ostringstream err;
00579     err << "Ifpack2::LocalFilter::getLocalRowCopy: Invalid local row index "
00580         << LocalRow << ".  The valid range of row indices on this process "
00581         << myRank << " is ";
00582     if (numRows == 0) {
00583       err << "empty";
00584     }
00585     else {
00586       err << "[0, " << (numRows - static_cast<size_t> (1)) << "]";
00587     }
00588     err << ".";
00589     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, err.str ());
00590   }
00591   TEUCHOS_TEST_FOR_EXCEPTION(
00592     (size_t) Indices.size() <  NumEntries_[LocalRow], std::runtime_error,
00593     "Ifpack2::LocalFilter::getLocalRowCopy: Invalid output array length.  "
00594     "The output arrays must each have length at least " << NumEntries_[LocalRow]
00595     << " for local row " << LocalRow << " on process " << myRank << ".");
00596 
00597   size_t A_NumEntries=0;
00598   // Always extract using the object Values_ and localIndices_.  This is
00599   // because I may need more space than that given by the user.  The
00600   // users expects only the local (in the domain Map) column indices,
00601   // but I have to extract both local and remote (not in the domain
00602   // Map) column indices.
00603   A_->getLocalRowCopy (LocalRow, localIndices_ (), Values_ (), A_NumEntries);
00604 
00605   // populate the user's vectors
00606   NumEntries = 0;
00607   for (size_t j = 0 ; j < A_NumEntries; ++j) {
00608     // only local indices
00609     if ((size_t) localIndices_[j] < numRows) {
00610       Indices[NumEntries] = localIndices_[j];
00611       Values[NumEntries]  = Values_[j];
00612       NumEntries++;
00613     }
00614   }
00615 }
00616 
00617 
00618 template<class MatrixType>
00619 void
00620 LocalFilter<MatrixType>::
00621 getGlobalRowView (global_ordinal_type GlobalRow,
00622                   Teuchos::ArrayView<const global_ordinal_type> &indices,
00623                   Teuchos::ArrayView<const scalar_type> &values) const
00624 {
00625   TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00626     "Ifpack2::LocalFilter does not implement getGlobalRowView.");
00627 }
00628 
00629 
00630 template<class MatrixType>
00631 void
00632 LocalFilter<MatrixType>::
00633 getLocalRowView (local_ordinal_type LocalRow,
00634                  Teuchos::ArrayView<const local_ordinal_type> &indices,
00635                  Teuchos::ArrayView<const scalar_type> &values) const
00636 {
00637   TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00638     "Ifpack2::LocalFilter does not implement getLocalRowView.");
00639 }
00640 
00641 
00642 template<class MatrixType>
00643 void
00644 LocalFilter<MatrixType>::
00645 getLocalDiagCopy (Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& diag) const
00646 {
00647   using Teuchos::RCP;
00648   typedef Tpetra::Vector<scalar_type, local_ordinal_type,
00649                          global_ordinal_type, node_type> vector_type;
00650   // This is always correct, and doesn't require a collective check
00651   // for sameness of Maps.
00652   RCP<vector_type> diagView = diag.offsetViewNonConst (A_->getRowMap (), 0);
00653   A_->getLocalDiagCopy (*diagView);
00654 }
00655 
00656 
00657 template<class MatrixType>
00658 void
00659 LocalFilter<MatrixType>::
00660 leftScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
00661 {
00662   TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00663     "Ifpack2::LocalFilter does not implement leftScale.");
00664 }
00665 
00666 
00667 template<class MatrixType>
00668 void
00669 LocalFilter<MatrixType>::
00670 rightScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
00671 {
00672   TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00673     "Ifpack2::LocalFilter does not implement rightScale.");
00674 }
00675 
00676 
00677 template<class MatrixType>
00678 void
00679 LocalFilter<MatrixType>::
00680 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
00681        Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
00682        Teuchos::ETransp mode,
00683        scalar_type alpha,
00684        scalar_type beta) const
00685 {
00686   using Teuchos::as;
00687   using Teuchos::ArrayView;
00688   using Teuchos::ArrayRCP;
00689   typedef Teuchos::ScalarTraits<scalar_type> STS;
00690 
00691   TEUCHOS_TEST_FOR_EXCEPTION(
00692     X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00693     "Ifpack2::LocalFilter::apply: X and Y must have the same number of columns.  "
00694     "X has " << X.getNumVectors () << " columns, but Y has "
00695     << Y.getNumVectors () << " columns.");
00696 
00697   const scalar_type zero = STS::zero ();
00698   const scalar_type one = STS::one ();
00699 
00700   if (beta == zero) {
00701     Y.putScalar (zero);
00702   }
00703   else if (beta != one) {
00704     Y.scale (beta);
00705   }
00706 
00707   const size_t NumVectors = Y.getNumVectors ();
00708   const size_t numRows = localRowMap_->getNodeNumElements ();
00709 
00710   // FIXME (mfh 14 Apr 2014) At some point, we would like to
00711   // parallelize this using Kokkos.  This would require a
00712   // Kokkos-friendly version of getLocalRowCopy, perhaps with
00713   // thread-local storage.
00714 
00715   const bool constantStride = X.isConstantStride () && Y.isConstantStride ();
00716   if (constantStride) {
00717     // Since both X and Y have constant stride, we can work with "1-D"
00718     // views of their data.
00719     const size_t              x_stride = X.getStride();
00720     const size_t              y_stride = Y.getStride();
00721     ArrayRCP<scalar_type>        y_rcp = Y.get1dViewNonConst();
00722     ArrayRCP<const scalar_type>  x_rcp = X.get1dView();
00723     ArrayView<scalar_type>       y_ptr = y_rcp();
00724     ArrayView<const scalar_type> x_ptr = x_rcp();
00725     for (size_t i = 0; i < numRows; ++i) {
00726       size_t Nnz;
00727       // Use this class's getrow to make the below code simpler
00728       getLocalRowCopy (i, localIndices_ (), Values_ (), Nnz);
00729       if (mode == Teuchos::NO_TRANS) {
00730         for (size_t j = 0; j < Nnz; ++j) {
00731           const local_ordinal_type col = localIndices_[j];
00732           for (size_t k = 0; k < NumVectors; ++k) {
00733             y_ptr[i + y_stride*k] +=
00734               alpha * Values_[j] * x_ptr[col + x_stride*k];
00735           }
00736         }
00737       }
00738       else if (mode == Teuchos::TRANS) {
00739         for (size_t j = 0; j < Nnz; ++j) {
00740           const local_ordinal_type col = localIndices_[j];
00741           for (size_t k = 0; k < NumVectors; ++k) {
00742             y_ptr[col + y_stride*k] +=
00743               alpha * Values_[j] * x_ptr[i + x_stride*k];
00744           }
00745         }
00746       }
00747       else { //mode==Teuchos::CONJ_TRANS
00748         for (size_t j = 0; j < Nnz; ++j) {
00749           const local_ordinal_type col = localIndices_[j];
00750           for (size_t k = 0; k < NumVectors; ++k) {
00751             y_ptr[col + y_stride*k] +=
00752               alpha * STS::conjugate (Values_[j]) * x_ptr[i + x_stride*k];
00753           }
00754         }
00755       }
00756     }
00757   }
00758   else {
00759     // At least one of X or Y does not have constant stride.
00760     // This means we must work with "2-D" views of their data.
00761     ArrayRCP<ArrayRCP<const scalar_type> > x_ptr = X.get2dView();
00762     ArrayRCP<ArrayRCP<scalar_type> >       y_ptr = Y.get2dViewNonConst();
00763 
00764     for (size_t i = 0; i < numRows; ++i) {
00765       size_t Nnz;
00766       // Use this class's getrow to make the below code simpler
00767       getLocalRowCopy (i, localIndices_ (), Values_ (), Nnz);
00768       if (mode == Teuchos::NO_TRANS) {
00769         for (size_t k = 0; k < NumVectors; ++k) {
00770           ArrayView<const scalar_type> x_local = (x_ptr())[k]();
00771           ArrayView<scalar_type>       y_local = (y_ptr())[k]();
00772           for (size_t j = 0; j < Nnz; ++j) {
00773             y_local[i] += alpha * Values_[j] * x_local[localIndices_[j]];
00774           }
00775         }
00776       }
00777       else if (mode == Teuchos::TRANS) {
00778         for (size_t k = 0; k < NumVectors; ++k) {
00779           ArrayView<const scalar_type> x_local = (x_ptr())[k]();
00780           ArrayView<scalar_type>       y_local = (y_ptr())[k]();
00781           for (size_t j = 0; j < Nnz; ++j) {
00782             y_local[localIndices_[j]] += alpha * Values_[j] * x_local[i];
00783           }
00784         }
00785       }
00786       else { //mode==Teuchos::CONJ_TRANS
00787         for (size_t k = 0; k < NumVectors; ++k) {
00788           ArrayView<const scalar_type> x_local = (x_ptr())[k]();
00789           ArrayView<scalar_type>       y_local = (y_ptr())[k]();
00790           for (size_t j = 0; j < Nnz; ++j) {
00791             y_local[localIndices_[j]] +=
00792               alpha * STS::conjugate (Values_[j]) * x_local[i];
00793           }
00794         }
00795       }
00796     }
00797   }
00798 }
00799 
00800 template<class MatrixType>
00801 bool LocalFilter<MatrixType>::hasTransposeApply () const
00802 {
00803   return true;
00804 }
00805 
00806 
00807 template<class MatrixType>
00808 bool LocalFilter<MatrixType>::supportsRowViews () const
00809 {
00810   return false;
00811 }
00812 
00813 
00814 template<class MatrixType>
00815 typename
00816 Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00817 LocalFilter<MatrixType>::getFrobeniusNorm () const
00818 {
00819   using Teuchos::as;
00820   typedef Teuchos::ScalarTraits<scalar_type> STS;
00821   typedef Teuchos::ScalarTraits<magnitude_type> STM;
00822   typedef typename Teuchos::Array<scalar_type>::size_type size_type;
00823 
00824   const size_type maxNumRowEnt = getNodeMaxNumRowEntries ();
00825   Teuchos::Array<local_ordinal_type> ind (maxNumRowEnt);
00826   Teuchos::Array<scalar_type> val (maxNumRowEnt);
00827   const size_t numRows = as<size_t> (localRowMap_->getNodeNumElements ());
00828 
00829   // FIXME (mfh 03 Apr 2013) Scale during sum to avoid overflow.
00830   magnitude_type sumSquared = STM::zero ();
00831   if (STS::isComplex) {
00832     for (size_t i = 0; i < numRows; ++i) {
00833       size_t numEntries = 0;
00834       this->getLocalRowCopy (i, ind (), val (), numEntries);
00835       for (size_type k = 0; k < static_cast<size_type> (numEntries); ++k) {
00836         sumSquared += STS::real (val[k]) * STS::real (val[k]) +
00837           STS::imag (val[k]) * STS::imag (val[k]);
00838       }
00839     }
00840   }
00841   else {
00842     for (size_t i = 0; i < numRows; ++i) {
00843       size_t numEntries = 0;
00844       this->getLocalRowCopy (i, ind (), val (), numEntries);
00845       for (size_type k = 0; k < static_cast<size_type> (numEntries); ++k) {
00846         const magnitude_type v_k_abs = STS::magnitude (val[k]);
00847         sumSquared += v_k_abs * v_k_abs;
00848       }
00849     }
00850   }
00851   return STM::squareroot (sumSquared); // Different for each process; that's OK.
00852 }
00853 
00854 
00855 template<class MatrixType>
00856 TPETRA_DEPRECATED void
00857 LocalFilter<MatrixType>::
00858 getGlobalRowView (global_ordinal_type GlobalRow,
00859                   Teuchos::ArrayRCP<const global_ordinal_type> &indices,
00860                   Teuchos::ArrayRCP<const scalar_type> &values) const
00861 {
00862   TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00863     "Ifpack2::LocalFilter does not implement getGlobalRowView.");
00864 }
00865 
00866 
00867 template<class MatrixType>
00868 TPETRA_DEPRECATED
00869 void
00870 LocalFilter<MatrixType>::
00871 getLocalRowView (local_ordinal_type LocalRow,
00872                  Teuchos::ArrayRCP<const local_ordinal_type> &indices,
00873                  Teuchos::ArrayRCP<const scalar_type> &values) const
00874 {
00875   TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00876     "Ifpack2::LocalFilter does not implement getLocalRowView.");
00877 }
00878 
00879 
00880 template<class MatrixType>
00881 std::string
00882 LocalFilter<MatrixType>::description () const
00883 {
00884   using Teuchos::TypeNameTraits;
00885   std::ostringstream os;
00886 
00887   os << "Ifpack2::LocalFilter: {";
00888   os << "MatrixType: " << TypeNameTraits<MatrixType>::name ();
00889   if (this->getObjectLabel () != "") {
00890     os << ", Label: \"" << this->getObjectLabel () << "\"";
00891   }
00892   os << ", Number of rows: " << getGlobalNumRows ()
00893      << ", Number of columns: " << getGlobalNumCols ()
00894      << "}";
00895   return os.str ();
00896 }
00897 
00898 
00899 template<class MatrixType>
00900 void
00901 LocalFilter<MatrixType>::
00902 describe (Teuchos::FancyOStream &out,
00903           const Teuchos::EVerbosityLevel verbLevel) const
00904 {
00905   using Teuchos::OSTab;
00906   using Teuchos::TypeNameTraits;
00907   using std::endl;
00908 
00909   const Teuchos::EVerbosityLevel vl =
00910     (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
00911 
00912   if (vl > Teuchos::VERB_NONE) {
00913     // describe() starts with a tab, by convention.
00914     OSTab tab0 (out);
00915 
00916     out << "Ifpack2::LocalFilter:" << endl;
00917     OSTab tab1 (out);
00918     out << "MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
00919     if (this->getObjectLabel () != "") {
00920       out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
00921     }
00922     out << "Number of rows: " << getGlobalNumRows () << endl
00923         << "Number of columns: " << getGlobalNumCols () << endl
00924         << "Number of nonzeros: " << NumNonzeros_ << endl;
00925 
00926     if (vl > Teuchos::VERB_LOW) {
00927       out << "Row Map:" << endl;
00928       localRowMap_->describe (out, vl);
00929       out << "Domain Map:" << endl;
00930       localDomainMap_->describe (out, vl);
00931       out << "Range Map:" << endl;
00932       localRangeMap_->describe (out, vl);
00933     }
00934   }
00935 }
00936 
00937 
00938 } // namespace Ifpack2
00939 
00940 #define IFPACK2_LOCALFILTER_INSTANT(S,LO,GO,N)                            \
00941   template class Ifpack2::LocalFilter< Tpetra::CrsMatrix<S, LO, GO, N> >;
00942 
00943 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends