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