Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_LocalFilter_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #ifndef IFPACK2_LOCALFILTER_DEF_HPP
44 #define IFPACK2_LOCALFILTER_DEF_HPP
45 
46 #include <Ifpack2_LocalFilter_decl.hpp>
47 #include <Tpetra_Map.hpp>
48 #include <Tpetra_MultiVector.hpp>
49 #include <Tpetra_Vector.hpp>
50 
51 #ifdef HAVE_MPI
52 # include "Teuchos_DefaultMpiComm.hpp"
53 #else
54 # include "Teuchos_DefaultSerialComm.hpp"
55 #endif
56 
57 namespace Ifpack2 {
58 
59 
60 template<class MatrixType>
61 bool
62 LocalFilter<MatrixType>::
63 mapPairsAreFitted (const row_matrix_type& A)
64 {
65  const map_type& rangeMap = * (A.getRangeMap ());
66  const map_type& rowMap = * (A.getRowMap ());
67  const bool rangeAndRowFitted = mapPairIsFitted (rangeMap, rowMap);
68 
69  const map_type& domainMap = * (A.getDomainMap ());
70  const map_type& columnMap = * (A.getColMap ());
71  const bool domainAndColumnFitted = mapPairIsFitted (domainMap, columnMap);
72 
73  return rangeAndRowFitted && domainAndColumnFitted;
74 }
75 
76 
77 template<class MatrixType>
78 bool
79 LocalFilter<MatrixType>::
80 mapPairIsFitted (const map_type& map1, const map_type& map2)
81 {
82  return map1.isLocallyFitted (map2);
83 }
84 
85 
86 template<class MatrixType>
88 LocalFilter (const Teuchos::RCP<const row_matrix_type>& A) :
89  A_ (A),
90  NumNonzeros_ (0),
91  MaxNumEntries_ (0),
92  MaxNumEntriesA_ (0)
93 {
94  using Teuchos::RCP;
95  using Teuchos::rcp;
96 
97 #ifdef HAVE_IFPACK2_DEBUG
98  TEUCHOS_TEST_FOR_EXCEPTION(
99  ! mapPairsAreFitted (*A), std::invalid_argument, "Ifpack2::LocalFilter: "
100  "A's Map pairs are not fitted to each other on Process "
101  << A_->getRowMap ()->getComm ()->getRank () << " of the input matrix's "
102  "communicator. "
103  "This means that LocalFilter does not currently know how to work with A. "
104  "This will change soon. Please see discussion of Bug 5992.");
105 #endif // HAVE_IFPACK2_DEBUG
106 
107  // Build the local communicator (containing this process only).
108  RCP<const Teuchos::Comm<int> > localComm;
109 #ifdef HAVE_MPI
110  localComm = rcp (new Teuchos::MpiComm<int> (MPI_COMM_SELF));
111 #else
112  localComm = rcp (new Teuchos::SerialComm<int> ());
113 #endif // HAVE_MPI
114 
115 
116  // // FIXME (mfh 21 Nov 2013) Currently, the implementation implicitly
117  // // assumes that the range Map is fitted to the row Map. Otherwise,
118  // // it probably won't work at all.
119  // TEUCHOS_TEST_FOR_EXCEPTION(
120  // mapPairIsFitted (* (A_->getRangeMap ()), * (A_->getRowMap ())),
121  // std::logic_error, "Ifpack2::LocalFilter: Range Map of the input matrix "
122  // "is not fitted to its row Map. LocalFilter does not currently work in "
123  // "this case. See Bug 5992.");
124 
125  // Build the local row, domain, and range Maps. They both use the
126  // local communicator built above. The global indices of each are
127  // different than those of the corresponding original Map; they
128  // actually are the same as the local indices of the original Map.
129  //
130  // It's even OK if signedness of local_ordinal_type and
131  // global_ordinal_type are different. (That would be a BAD IDEA but
132  // some users seem to enjoy making trouble for themselves and us.)
133  //
134  // Both the local row and local range Maps must have the same number
135  // of entries, namely, that of the local number of entries of A's
136  // range Map.
137 
138  const size_t numRows = A_->getRangeMap()->getNodeNumElements ();
139 
140  // using std::cerr;
141  // using std::endl;
142  // cerr << "A_ has " << A_->getNodeNumRows () << " rows." << endl
143  // << "Range Map has " << A_->getRangeMap ()->getNodeNumElements () << " entries." << endl
144  // << "Row Map has " << A_->getRowMap ()->getNodeNumElements () << " entries." << endl;
145 
146  const global_ordinal_type indexBase = static_cast<global_ordinal_type> (0);
147 
148  localRowMap_ =
149  rcp (new map_type (numRows, indexBase, localComm,
150  Tpetra::GloballyDistributed));
151  // If the original matrix's range Map is not fitted to its row Map,
152  // we'll have to do an Export when applying the matrix.
153  localRangeMap_ = localRowMap_;
154 
155  // If the original matrix's domain Map is not fitted to its column
156  // Map, we'll have to do an Import when applying the matrix.
157  if (A_->getRangeMap ().getRawPtr () == A_->getDomainMap ().getRawPtr ()) {
158  // The input matrix's domain and range Maps are identical, so the
159  // locally filtered matrix's domain and range Maps can be also.
160  localDomainMap_ = localRangeMap_;
161  }
162  else {
163  const size_t numCols = A_->getDomainMap()->getNodeNumElements ();
164  localDomainMap_ =
165  rcp (new map_type (numCols, indexBase, localComm,
166  Tpetra::GloballyDistributed));
167  }
168 
169  // NodeNumEntries_ will contain the actual number of nonzeros for
170  // each localized row (that is, without external nodes, and always
171  // with the diagonal entry)
172  NumEntries_.resize (numRows);
173 
174  // tentative value for MaxNumEntries. This is the number of
175  // nonzeros in the local matrix
176  MaxNumEntries_ = A_->getNodeMaxNumRowEntries ();
177  MaxNumEntriesA_ = A_->getNodeMaxNumRowEntries ();
178 
179  // Allocate temporary arrays for getLocalRowCopy().
180  Kokkos::resize(localIndices_,MaxNumEntries_);
181  Kokkos::resize(localIndicesForGlobalCopy_,MaxNumEntries_);
182  Kokkos::resize(Values_,MaxNumEntries_);
183 
184  // now compute:
185  // - the number of nonzero per row
186  // - the total number of nonzeros
187  // - the diagonal entries
188 
189  // compute nonzeros (total and per-row), and store the
190  // diagonal entries (already modified)
191  size_t ActualMaxNumEntries = 0;
192 
193  for (size_t i = 0; i < numRows; ++i) {
194  NumEntries_[i] = 0;
195  size_t Nnz, NewNnz = 0;
196  A_->getLocalRowCopy (i, localIndices_, Values_, Nnz);
197  for (size_t j = 0; j < Nnz; ++j) {
198  // FIXME (mfh 03 Apr 2013) This assumes the following:
199  //
200  // 1. Row Map, range Map, and domain Map are all the same.
201  //
202  // 2. The column Map's list of GIDs on this process is the
203  // domain Map's list of GIDs, followed by remote GIDs. Thus,
204  // for any GID in the domain Map on this process, its LID in
205  // the domain Map (and therefore in the row Map, by (1)) is
206  // the same as its LID in the column Map. (Hence the
207  // less-than test, which if true, means that localIndices_[j]
208  // belongs to the row Map.)
209  if (static_cast<size_t> (localIndices_[j]) < numRows) {
210  ++NewNnz;
211  }
212  }
213 
214  if (NewNnz > ActualMaxNumEntries) {
215  ActualMaxNumEntries = NewNnz;
216  }
217 
218  NumNonzeros_ += NewNnz;
219  NumEntries_[i] = NewNnz;
220  }
221 
222  MaxNumEntries_ = ActualMaxNumEntries;
223 }
224 
225 
226 template<class MatrixType>
228 {}
229 
230 
231 template<class MatrixType>
232 Teuchos::RCP<const Teuchos::Comm<int> >
234 {
235  return localRowMap_->getComm ();
236 }
237 
238 
239 
240 
241 template<class MatrixType>
242 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
243  typename MatrixType::global_ordinal_type,
244  typename MatrixType::node_type> >
246 {
247  return localRowMap_;
248 }
249 
250 
251 template<class MatrixType>
252 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
253  typename MatrixType::global_ordinal_type,
254  typename MatrixType::node_type> >
256 {
257  return localRowMap_; // FIXME (mfh 20 Nov 2013)
258 }
259 
260 
261 template<class MatrixType>
262 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
263  typename MatrixType::global_ordinal_type,
264  typename MatrixType::node_type> >
266 {
267  return localDomainMap_;
268 }
269 
270 
271 template<class MatrixType>
272 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
273  typename MatrixType::global_ordinal_type,
274  typename MatrixType::node_type> >
276 {
277  return localRangeMap_;
278 }
279 
280 
281 template<class MatrixType>
282 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type,
283  typename MatrixType::global_ordinal_type,
284  typename MatrixType::node_type> >
286 {
287  // FIXME (mfh 20 Nov 2013) This is not what the documentation says
288  // this method should do! It should return the graph of the locally
289  // filtered matrix, not the original matrix's graph.
290  return A_->getGraph ();
291 }
292 
293 
294 template<class MatrixType>
296 {
297  return static_cast<global_size_t> (localRangeMap_->getNodeNumElements ());
298 }
299 
300 
301 template<class MatrixType>
303 {
304  return static_cast<global_size_t> (localDomainMap_->getNodeNumElements ());
305 }
306 
307 
308 template<class MatrixType>
310 {
311  return static_cast<size_t> (localRangeMap_->getNodeNumElements ());
312 }
313 
314 
315 template<class MatrixType>
317 {
318  return static_cast<size_t> (localDomainMap_->getNodeNumElements ());
319 }
320 
321 
322 template<class MatrixType>
323 typename MatrixType::global_ordinal_type
325 {
326  return A_->getIndexBase ();
327 }
328 
329 
330 template<class MatrixType>
332 {
333  return NumNonzeros_;
334 }
335 
336 
337 template<class MatrixType>
339 {
340  return NumNonzeros_;
341 }
342 
343 
344 template<class MatrixType>
345 size_t
347 getNumEntriesInGlobalRow (global_ordinal_type globalRow) const
348 {
349  const local_ordinal_type localRow = getRowMap ()->getLocalElement (globalRow);
350  if (localRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
351  // NOTE (mfh 26 Mar 2014) We return zero if globalRow is not in
352  // the row Map on this process, since "get the number of entries
353  // in the global row" refers only to what the calling process owns
354  // in that row. In this case, it owns no entries in that row,
355  // since it doesn't own the row.
356  return 0;
357  } else {
358  return NumEntries_[localRow];
359  }
360 }
361 
362 
363 template<class MatrixType>
364 size_t
366 getNumEntriesInLocalRow (local_ordinal_type localRow) const
367 {
368  // FIXME (mfh 07 Jul 2014) Shouldn't localRow be a local row index
369  // in the matrix's row Map, not in the LocalFilter's row Map? The
370  // latter is different; it even has different global indices!
371  // (Maybe _that_'s the bug.)
372 
373  if (getRowMap ()->isNodeLocalElement (localRow)) {
374  return NumEntries_[localRow];
375  } else {
376  // NOTE (mfh 26 Mar 2014) We return zero if localRow is not in the
377  // row Map on this process, since "get the number of entries in
378  // the local row" refers only to what the calling process owns in
379  // that row. In this case, it owns no entries in that row, since
380  // it doesn't own the row.
381  return 0;
382  }
383 }
384 
385 
386 template<class MatrixType>
388 {
389  return MaxNumEntries_;
390 }
391 
392 
393 template<class MatrixType>
395 {
396  return MaxNumEntries_;
397 }
398 
399 
400 template<class MatrixType>
402 {
403  return true;
404 }
405 
406 
407 template<class MatrixType>
409 {
410  return A_->isLocallyIndexed ();
411 }
412 
413 
414 template<class MatrixType>
416 {
417  return A_->isGloballyIndexed();
418 }
419 
420 
421 template<class MatrixType>
423 {
424  return A_->isFillComplete ();
425 }
426 
427 
428 template<class MatrixType>
429 void
431  getGlobalRowCopy (global_ordinal_type globalRow,
432  nonconst_global_inds_host_view_type &globalIndices,
433  nonconst_values_host_view_type &values,
434  size_t& numEntries) const
435 {
436  typedef local_ordinal_type LO;
437  typedef typename Teuchos::Array<LO>::size_type size_type;
438 
439  const LO localRow = this->getRowMap ()->getLocalElement (globalRow);
440  if (localRow == Teuchos::OrdinalTraits<LO>::invalid ()) {
441  // NOTE (mfh 26 Mar 2014) We return no entries if globalRow is not
442  // in the row Map on this process, since "get a copy of the
443  // entries in the global row" refers only to what the calling
444  // process owns in that row. In this case, it owns no entries in
445  // that row, since it doesn't own the row.
446  numEntries = 0;
447  }
448  else {
449  // First get a copy of the current row using local indices. Then,
450  // convert to global indices using the input matrix's column Map.
451  //
452  numEntries = this->getNumEntriesInLocalRow (localRow);
453  // FIXME (mfh 26 Mar 2014) If local_ordinal_type ==
454  // global_ordinal_type, we could just alias the input array
455  // instead of allocating a temporary array.
456 
457  // In this case, getLocalRowCopy *does* use the localIndices_, so we use a second temp array
458  this->getLocalRowCopy (localRow, localIndicesForGlobalCopy_, values, numEntries);
459 
460  const map_type& colMap = * (this->getColMap ());
461 
462  // Don't fill the output array beyond its size.
463  const size_type numEnt =
464  std::min (static_cast<size_type> (numEntries),
465  std::min ((size_type)globalIndices.size (), (size_type)values.size ()));
466  for (size_type k = 0; k < numEnt; ++k) {
467  globalIndices[k] = colMap.getGlobalElement (localIndicesForGlobalCopy_[k]);
468  }
469  }
470 }
471 
472 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
473 template<class MatrixType>
474 void
476 getGlobalRowCopy (global_ordinal_type globalRow,
477  const Teuchos::ArrayView<global_ordinal_type>& Indices,
478  const Teuchos::ArrayView<scalar_type>& Values,
479  size_t& numEntries) const {
480  using IST = typename row_matrix_type::impl_scalar_type;
481  nonconst_global_inds_host_view_type ind_in(Indices.data(),Indices.size());
482  nonconst_values_host_view_type val_in(reinterpret_cast<IST*>(Values.data()),Values.size());
483  getGlobalRowCopy(globalRow,ind_in,val_in,numEntries);
484 }
485 #endif
486 
487 template<class MatrixType>
488 void
490 getLocalRowCopy (local_ordinal_type LocalRow,
491  nonconst_local_inds_host_view_type &Indices,
492  nonconst_values_host_view_type &Values,
493  size_t& NumEntries) const
494 {
495  typedef local_ordinal_type LO;
496  typedef global_ordinal_type GO;
497 
498  if (! A_->getRowMap ()->isNodeLocalElement (LocalRow)) {
499  // The calling process owns zero entries in the row.
500  NumEntries = 0;
501  return;
502  }
503 
504  if (A_->getRowMap()->getComm()->getSize() == 1) {
505  A_->getLocalRowCopy (LocalRow, Indices, Values, NumEntries);
506  return;
507  }
508 
509 
510  const size_t numEntInLclRow = NumEntries_[LocalRow];
511  if (static_cast<size_t> (Indices.size ()) < numEntInLclRow ||
512  static_cast<size_t> (Values.size ()) < numEntInLclRow) {
513  // FIXME (mfh 07 Jul 2014) Return an error code instead of
514  // throwing. We should really attempt to fill as much space as
515  // we're given, in this case.
516  TEUCHOS_TEST_FOR_EXCEPTION(
517  true, std::runtime_error,
518  "Ifpack2::LocalFilter::getLocalRowCopy: Invalid output array length. "
519  "The output arrays must each have length at least " << numEntInLclRow
520  << " for local row " << LocalRow << " on Process "
521  << localRowMap_->getComm ()->getRank () << ".");
522  }
523  else if (numEntInLclRow == static_cast<size_t> (0)) {
524  // getNumEntriesInLocalRow() returns zero if LocalRow is not owned
525  // by the calling process. In that case, the calling process owns
526  // zero entries in the row.
527  NumEntries = 0;
528  return;
529  }
530 
531  // Always extract using the temporary arrays Values_ and
532  // localIndices_. This is because I may need more space than that
533  // given by the user. The users expects only the local (in the
534  // domain Map) column indices, but I have to extract both local and
535  // remote (not in the domain Map) column indices.
536  //
537  // FIXME (mfh 07 Jul 2014) Use of temporary local storage is not
538  // conducive to thread parallelism. A better way would be to change
539  // the interface so that it only extracts values for the "local"
540  // column indices. CrsMatrix could take a set of column indices,
541  // and return their corresponding values.
542  size_t numEntInMat = 0;
543  A_->getLocalRowCopy (LocalRow, localIndices_, Values_ , numEntInMat);
544 
545  // Fill the user's arrays with the "local" indices and values in
546  // that row. Note that the matrix might have a different column Map
547  // than the local filter.
548  const map_type& matrixDomMap = * (A_->getDomainMap ());
549  const map_type& matrixColMap = * (A_->getColMap ());
550 
551  const size_t capacity = static_cast<size_t> (std::min (Indices.size (),
552  Values.size ()));
553  NumEntries = 0;
554  const size_t numRows = localRowMap_->getNodeNumElements (); // superfluous
555  const bool buggy = true; // mfh 07 Jul 2014: See FIXME below.
556  for (size_t j = 0; j < numEntInMat; ++j) {
557  // The LocalFilter only includes entries in the domain Map on
558  // the calling process. We figure out whether an entry is in
559  // the domain Map by converting the (matrix column Map) index to
560  // a global index, and then asking whether that global index is
561  // in the domain Map.
562  const LO matrixLclCol = localIndices_[j];
563  const GO gblCol = matrixColMap.getGlobalElement (matrixLclCol);
564 
565  // FIXME (mfh 07 Jul 2014) This is the likely center of Bug 5992
566  // and perhaps other bugs, like Bug 6117. If 'buggy' is true,
567  // Ifpack2 tests pass; if 'buggy' is false, the tests don't pass.
568  // This suggests that Ifpack2 classes could be using LocalFilter
569  // incorrectly, perhaps by giving it an incorrect domain Map.
570  if (buggy) {
571  // only local indices
572  if ((size_t) localIndices_[j] < numRows) {
573  Indices[NumEntries] = localIndices_[j];
574  Values[NumEntries] = Values_[j];
575  NumEntries++;
576  }
577  } else {
578  if (matrixDomMap.isNodeGlobalElement (gblCol)) {
579  // Don't fill more space than the user gave us. It's an error
580  // for them not to give us enough space, but we still shouldn't
581  // overwrite memory that doesn't belong to us.
582  if (NumEntries < capacity) {
583  Indices[NumEntries] = matrixLclCol;
584  Values[NumEntries] = Values_[j];
585  }
586  NumEntries++;
587  }
588  }
589  }
590 }
591 
592 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
593 template<class MatrixType>
594 void
596 getLocalRowCopy (local_ordinal_type globalRow,
597  const Teuchos::ArrayView<local_ordinal_type> &Indices,
598  const Teuchos::ArrayView<scalar_type> &Values,
599  size_t &NumEntries) const
600 {
601  using IST = typename row_matrix_type::impl_scalar_type;
602  nonconst_local_inds_host_view_type ind_in(Indices.data(),Indices.size());
603  nonconst_values_host_view_type val_in(reinterpret_cast<IST*>(Values.data()),Values.size());
604  getLocalRowCopy(globalRow,ind_in,val_in,NumEntries);
605 }
606 #endif
607 
608 
609 template<class MatrixType>
610 void
612 getGlobalRowView (global_ordinal_type /*GlobalRow*/,
613  global_inds_host_view_type &/*indices*/,
614  values_host_view_type &/*values*/) const
615 {
616  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
617  "Ifpack2::LocalFilter does not implement getGlobalRowView.");
618 }
619 
620 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
621 template<class MatrixType>
622 void
624 getGlobalRowView (global_ordinal_type /* GlobalRow */,
625  Teuchos::ArrayView<const global_ordinal_type> &/* indices */,
626  Teuchos::ArrayView<const scalar_type> &/* values */) const
627 {
628  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
629  "Ifpack2::LocalFilter does not implement getGlobalRowView.");
630 }
631 #endif
632 
633 template<class MatrixType>
634 void
636 getLocalRowView (local_ordinal_type /*LocalRow*/,
637  local_inds_host_view_type &/*indices*/,
638  values_host_view_type &/*values*/) const
639 {
640  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
641  "Ifpack2::LocalFilter does not implement getLocalRowView.");
642 }
643 
644 
645 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
646 template<class MatrixType>
647 void
649 getLocalRowView (local_ordinal_type /* LocalRow */,
650  Teuchos::ArrayView<const local_ordinal_type> &/* indices */,
651  Teuchos::ArrayView<const scalar_type> &/* values */) const
652 {
653  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
654  "Ifpack2::LocalFilter does not implement getLocalRowView.");
655 }
656 #endif
657 
658 
659 template<class MatrixType>
660 void
662 getLocalDiagCopy (Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& diag) const
663 {
664  using Teuchos::RCP;
665  typedef Tpetra::Vector<scalar_type, local_ordinal_type,
666  global_ordinal_type, node_type> vector_type;
667  // This is always correct, and doesn't require a collective check
668  // for sameness of Maps.
669  RCP<vector_type> diagView = diag.offsetViewNonConst (A_->getRowMap (), 0);
670  A_->getLocalDiagCopy (*diagView);
671 }
672 
673 
674 template<class MatrixType>
675 void
677 leftScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& /* x */)
678 {
679  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
680  "Ifpack2::LocalFilter does not implement leftScale.");
681 }
682 
683 
684 template<class MatrixType>
685 void
687 rightScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& /* x */)
688 {
689  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
690  "Ifpack2::LocalFilter does not implement rightScale.");
691 }
692 
693 
694 template<class MatrixType>
695 void
697 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
698  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
699  Teuchos::ETransp mode,
700  scalar_type alpha,
701  scalar_type beta) const
702 {
703  typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV;
704  TEUCHOS_TEST_FOR_EXCEPTION(
705  X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
706  "Ifpack2::LocalFilter::apply: X and Y must have the same number of columns. "
707  "X has " << X.getNumVectors () << " columns, but Y has "
708  << Y.getNumVectors () << " columns.");
709 
710 #ifdef HAVE_IFPACK2_DEBUG
711  {
712  typedef Teuchos::ScalarTraits<magnitude_type> STM;
713  Teuchos::Array<magnitude_type> norms (X.getNumVectors ());
714  X.norm1 (norms ());
715  bool good = true;
716  for (size_t j = 0; j < X.getNumVectors (); ++j) {
717  if (STM::isnaninf (norms[j])) {
718  good = false;
719  break;
720  }
721  }
722  TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error, "Ifpack2::LocalFilter::apply: The 1-norm of the input X is NaN or Inf.");
723  }
724 #endif // HAVE_IFPACK2_DEBUG
725 
726  if (&X == &Y) {
727  // FIXME (mfh 23 Apr 2014) We have to do more work to figure out
728  // if X and Y alias one another. For example, we should check
729  // whether one is a noncontiguous view of the other.
730  //
731  // X and Y alias one another, so we have to copy X.
732  MV X_copy (X, Teuchos::Copy);
733  applyNonAliased (X_copy, Y, mode, alpha, beta);
734  } else {
735  applyNonAliased (X, Y, mode, alpha, beta);
736  }
737 
738 #ifdef HAVE_IFPACK2_DEBUG
739  {
740  typedef Teuchos::ScalarTraits<magnitude_type> STM;
741  Teuchos::Array<magnitude_type> norms (Y.getNumVectors ());
742  Y.norm1 (norms ());
743  bool good = true;
744  for (size_t j = 0; j < Y.getNumVectors (); ++j) {
745  if (STM::isnaninf (norms[j])) {
746  good = false;
747  break;
748  }
749  }
750  TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error, "Ifpack2::LocalFilter::apply: The 1-norm of the output Y is NaN or Inf.");
751  }
752 #endif // HAVE_IFPACK2_DEBUG
753 }
754 
755 template<class MatrixType>
756 void
758 applyNonAliased (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
759  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
760  Teuchos::ETransp mode,
761  scalar_type alpha,
762  scalar_type beta) const
763 {
764  using Teuchos::ArrayView;
765  using Teuchos::ArrayRCP;
766  typedef Teuchos::ScalarTraits<scalar_type> STS;
767 
768  const scalar_type zero = STS::zero ();
769  const scalar_type one = STS::one ();
770 
771  if (beta == zero) {
772  Y.putScalar (zero);
773  }
774  else if (beta != one) {
775  Y.scale (beta);
776  }
777 
778  const size_t NumVectors = Y.getNumVectors ();
779  const size_t numRows = localRowMap_->getNodeNumElements ();
780 
781  // FIXME (mfh 14 Apr 2014) At some point, we would like to
782  // parallelize this using Kokkos. This would require a
783  // Kokkos-friendly version of getLocalRowCopy, perhaps with
784  // thread-local storage.
785 
786  const bool constantStride = X.isConstantStride () && Y.isConstantStride ();
787  if (constantStride) {
788  // Since both X and Y have constant stride, we can work with "1-D"
789  // views of their data.
790  const size_t x_stride = X.getStride();
791  const size_t y_stride = Y.getStride();
792  ArrayRCP<scalar_type> y_rcp = Y.get1dViewNonConst();
793  ArrayRCP<const scalar_type> x_rcp = X.get1dView();
794  ArrayView<scalar_type> y_ptr = y_rcp();
795  ArrayView<const scalar_type> x_ptr = x_rcp();
796  for (size_t i = 0; i < numRows; ++i) {
797  size_t Nnz;
798  // Use this class's getrow to make the below code simpler
799  getLocalRowCopy (i, localIndices_ , Values_ , Nnz);
800  scalar_type* Values = reinterpret_cast<scalar_type*>(Values_.data());
801  if (mode == Teuchos::NO_TRANS) {
802  for (size_t j = 0; j < Nnz; ++j) {
803  const local_ordinal_type col = localIndices_[j];
804  for (size_t k = 0; k < NumVectors; ++k) {
805  y_ptr[i + y_stride*k] +=
806  alpha * Values[j] * x_ptr[col + x_stride*k];
807  }
808  }
809  }
810  else if (mode == Teuchos::TRANS) {
811  for (size_t j = 0; j < Nnz; ++j) {
812  const local_ordinal_type col = localIndices_[j];
813  for (size_t k = 0; k < NumVectors; ++k) {
814  y_ptr[col + y_stride*k] +=
815  alpha * Values[j] * x_ptr[i + x_stride*k];
816  }
817  }
818  }
819  else { //mode==Teuchos::CONJ_TRANS
820  for (size_t j = 0; j < Nnz; ++j) {
821  const local_ordinal_type col = localIndices_[j];
822  for (size_t k = 0; k < NumVectors; ++k) {
823  y_ptr[col + y_stride*k] +=
824  alpha * STS::conjugate (Values[j]) * x_ptr[i + x_stride*k];
825  }
826  }
827  }
828  }
829  }
830  else {
831  // At least one of X or Y does not have constant stride.
832  // This means we must work with "2-D" views of their data.
833  ArrayRCP<ArrayRCP<const scalar_type> > x_ptr = X.get2dView();
834  ArrayRCP<ArrayRCP<scalar_type> > y_ptr = Y.get2dViewNonConst();
835 
836  for (size_t i = 0; i < numRows; ++i) {
837  size_t Nnz;
838  // Use this class's getrow to make the below code simpler
839  getLocalRowCopy (i, localIndices_ , Values_ , Nnz);
840  scalar_type* Values = reinterpret_cast<scalar_type*>(Values_.data());
841  if (mode == Teuchos::NO_TRANS) {
842  for (size_t k = 0; k < NumVectors; ++k) {
843  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
844  ArrayView<scalar_type> y_local = (y_ptr())[k]();
845  for (size_t j = 0; j < Nnz; ++j) {
846  y_local[i] += alpha * Values[j] * x_local[localIndices_[j]];
847  }
848  }
849  }
850  else if (mode == Teuchos::TRANS) {
851  for (size_t k = 0; k < NumVectors; ++k) {
852  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
853  ArrayView<scalar_type> y_local = (y_ptr())[k]();
854  for (size_t j = 0; j < Nnz; ++j) {
855  y_local[localIndices_[j]] += alpha * Values[j] * x_local[i];
856  }
857  }
858  }
859  else { //mode==Teuchos::CONJ_TRANS
860  for (size_t k = 0; k < NumVectors; ++k) {
861  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
862  ArrayView<scalar_type> y_local = (y_ptr())[k]();
863  for (size_t j = 0; j < Nnz; ++j) {
864  y_local[localIndices_[j]] +=
865  alpha * STS::conjugate (Values[j]) * x_local[i];
866  }
867  }
868  }
869  }
870  }
871 }
872 
873 template<class MatrixType>
875 {
876  return true;
877 }
878 
879 
880 template<class MatrixType>
882 {
883  return false;
884 }
885 
886 
887 template<class MatrixType>
888 typename
889 LocalFilter<MatrixType>::mag_type
891 {
892 #ifdef TPETRA_HAVE_KOKKOS_REFACTOR
893  typedef Kokkos::Details::ArithTraits<scalar_type> STS;
894  typedef Kokkos::Details::ArithTraits<mag_type> STM;
895 #else
896  typedef Teuchos::ScalarTraits<scalar_type> STS;
897  typedef Teuchos::ScalarTraits<magnitude_type> STM;
898 #endif
899  typedef typename Teuchos::Array<scalar_type>::size_type size_type;
900 
901  const size_type maxNumRowEnt = getNodeMaxNumRowEntries ();
902  nonconst_local_inds_host_view_type ind ("ind",maxNumRowEnt);
903  nonconst_values_host_view_type val ("val",maxNumRowEnt);
904  const size_t numRows = static_cast<size_t> (localRowMap_->getNodeNumElements ());
905 
906  // FIXME (mfh 03 Apr 2013) Scale during sum to avoid overflow.
907  mag_type sumSquared = STM::zero ();
908  for (size_t i = 0; i < numRows; ++i) {
909  size_t numEntries = 0;
910  this->getLocalRowCopy (i, ind, val, numEntries);
911  for (size_type k = 0; k < static_cast<size_type> (numEntries); ++k) {
912  const mag_type v_k_abs = STS::magnitude (val[k]);
913  sumSquared += v_k_abs * v_k_abs;
914  }
915  }
916  return STM::squareroot (sumSquared); // Different for each process; that's OK.
917 }
918 
919 template<class MatrixType>
920 std::string
922 {
923  using Teuchos::TypeNameTraits;
924  std::ostringstream os;
925 
926  os << "Ifpack2::LocalFilter: {";
927  os << "MatrixType: " << TypeNameTraits<MatrixType>::name ();
928  if (this->getObjectLabel () != "") {
929  os << ", Label: \"" << this->getObjectLabel () << "\"";
930  }
931  os << ", Number of rows: " << getGlobalNumRows ()
932  << ", Number of columns: " << getGlobalNumCols ()
933  << "}";
934  return os.str ();
935 }
936 
937 
938 template<class MatrixType>
939 void
941 describe (Teuchos::FancyOStream &out,
942  const Teuchos::EVerbosityLevel verbLevel) const
943 {
944  using Teuchos::OSTab;
945  using Teuchos::TypeNameTraits;
946  using std::endl;
947 
948  const Teuchos::EVerbosityLevel vl =
949  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
950 
951  if (vl > Teuchos::VERB_NONE) {
952  // describe() starts with a tab, by convention.
953  OSTab tab0 (out);
954 
955  out << "Ifpack2::LocalFilter:" << endl;
956  OSTab tab1 (out);
957  out << "MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
958  if (this->getObjectLabel () != "") {
959  out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
960  }
961  out << "Number of rows: " << getGlobalNumRows () << endl
962  << "Number of columns: " << getGlobalNumCols () << endl
963  << "Number of nonzeros: " << NumNonzeros_ << endl;
964 
965  if (vl > Teuchos::VERB_LOW) {
966  out << "Row Map:" << endl;
967  localRowMap_->describe (out, vl);
968  out << "Domain Map:" << endl;
969  localDomainMap_->describe (out, vl);
970  out << "Range Map:" << endl;
971  localRangeMap_->describe (out, vl);
972  }
973  }
974 }
975 
976 template<class MatrixType>
977 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
978  typename MatrixType::local_ordinal_type,
979  typename MatrixType::global_ordinal_type,
980  typename MatrixType::node_type> >
982 {
983  return A_;
984 }
985 
986 
987 } // namespace Ifpack2
988 
989 #define IFPACK2_LOCALFILTER_INSTANT(S,LO,GO,N) \
990  template class Ifpack2::LocalFilter< Tpetra::RowMatrix<S, LO, GO, N> >;
991 
992 #endif
virtual void getGlobalRowCopy(global_ordinal_type GlobalRow, nonconst_global_inds_host_view_type &Indices, nonconst_values_host_view_type &Values, size_t &NumEntries) const
Get the entries in the given row, using global indices.
Definition: Ifpack2_LocalFilter_def.hpp:431
virtual bool isFillComplete() const
Returns true if fillComplete() has been called.
Definition: Ifpack2_LocalFilter_def.hpp:422
virtual void rightScale(const Tpetra::Vector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &x)
Scales the RowMatrix on the right with the Vector x.
Definition: Ifpack2_LocalFilter_def.hpp:687
virtual mag_type getFrobeniusNorm() const
The Frobenius norm of the (locally filtered) matrix.
Definition: Ifpack2_LocalFilter_def.hpp:890
virtual Teuchos::RCP< const row_matrix_type > getUnderlyingMatrix() const
Return matrix that LocalFilter was built on.
Definition: Ifpack2_LocalFilter_def.hpp:981
virtual global_size_t getGlobalNumCols() const
The number of global columns in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:302
virtual bool isLocallyIndexed() const
Whether the underlying sparse matrix is locally (opposite of globally) indexed.
Definition: Ifpack2_LocalFilter_def.hpp:408
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object to the given output stream.
Definition: Ifpack2_LocalFilter_def.hpp:941
virtual size_t getNodeNumCols() const
The number of columns in the (locally filtered) matrix.
Definition: Ifpack2_LocalFilter_def.hpp:316
virtual global_size_t getGlobalNumRows() const
The number of global rows in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:295
virtual Teuchos::RCP< const map_type > getRangeMap() const
Returns the Map that describes the range distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:275
Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Type of the Tpetra::Map specialization that this class uses.
Definition: Ifpack2_LocalFilter_decl.hpp:216
virtual bool supportsRowViews() const
Returns true if RowViews are supported.
Definition: Ifpack2_LocalFilter_def.hpp:881
virtual Teuchos::RCP< const map_type > getColMap() const
Returns the Map that describes the column distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:255
virtual void leftScale(const Tpetra::Vector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &x)
Scales the RowMatrix on the left with the Vector x.
Definition: Ifpack2_LocalFilter_def.hpp:677
virtual Teuchos::RCP< const Tpetra::RowGraph< local_ordinal_type, global_ordinal_type, node_type > > getGraph() const
The (locally filtered) matrix&#39;s graph.
Definition: Ifpack2_LocalFilter_def.hpp:285
virtual size_t getNodeNumRows() const
The number of rows owned on the calling process.
Definition: Ifpack2_LocalFilter_def.hpp:309
virtual void getLocalRowCopy(local_ordinal_type LocalRow, nonconst_local_inds_host_view_type &Indices, nonconst_values_host_view_type &Values, size_t &NumEntries) const
Get the entries in the given row, using local indices.
Definition: Ifpack2_LocalFilter_def.hpp:490
virtual void apply(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Compute Y = beta*Y + alpha*A_local*X.
Definition: Ifpack2_LocalFilter_def.hpp:697
virtual global_ordinal_type getIndexBase() const
Returns the index base for global indices for this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:324
virtual bool hasColMap() const
Whether this matrix has a well-defined column Map.
Definition: Ifpack2_LocalFilter_def.hpp:401
virtual size_t getGlobalMaxNumRowEntries() const
The maximum number of entries across all rows/columns on all processes.
Definition: Ifpack2_LocalFilter_def.hpp:387
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the communicator.
Definition: Ifpack2_LocalFilter_def.hpp:233
virtual std::string description() const
A one-line description of this object.
Definition: Ifpack2_LocalFilter_def.hpp:921
virtual size_t getNodeMaxNumRowEntries() const
The maximum number of entries across all rows/columns on this process.
Definition: Ifpack2_LocalFilter_def.hpp:394
LocalFilter(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_LocalFilter_def.hpp:88
virtual ~LocalFilter()
Destructor.
Definition: Ifpack2_LocalFilter_def.hpp:227
virtual global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:331
virtual Teuchos::RCP< const map_type > getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:245
virtual size_t getNodeNumEntries() const
Returns the local number of entries in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:338
virtual bool isGloballyIndexed() const
Whether the underlying sparse matrix is globally (opposite of locally) indexed.
Definition: Ifpack2_LocalFilter_def.hpp:415
virtual void getLocalDiagCopy(Tpetra::Vector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &diag) const
Get the diagonal entries of the (locally filtered) matrix.
Definition: Ifpack2_LocalFilter_def.hpp:662
Access only local rows and columns of a sparse matrix.
Definition: Ifpack2_LocalFilter_decl.hpp:160
virtual Teuchos::RCP< const map_type > getDomainMap() const
Returns the Map that describes the domain distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:265
virtual size_t getNumEntriesInGlobalRow(global_ordinal_type globalRow) const
The current number of entries on this node in the specified global row.
Definition: Ifpack2_LocalFilter_def.hpp:347
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
virtual void getLocalRowView(local_ordinal_type LocalRow, local_inds_host_view_type &indices, values_host_view_type &values) const
Extract a const, non-persisting view of local indices in a specified row of the matrix.
Definition: Ifpack2_LocalFilter_def.hpp:636
virtual bool hasTransposeApply() const
Whether this operator supports applying the transpose or conjugate transpose.
Definition: Ifpack2_LocalFilter_def.hpp:874
virtual void getGlobalRowView(global_ordinal_type GlobalRow, global_inds_host_view_type &indices, values_host_view_type &values) const
Extract a const, non-persisting view of global indices in a specified row of the matrix.
Definition: Ifpack2_LocalFilter_def.hpp:612
virtual size_t getNumEntriesInLocalRow(local_ordinal_type localRow) const
The current number of entries on this node in the specified local row.
Definition: Ifpack2_LocalFilter_def.hpp:366