Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_LocalFilter_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Tempated 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
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 Tpetra::Details::isLocallyFitted (map1, map2);
83 }
84 
85 
86 template<class MatrixType>
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
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, A_->getNode ()));
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, A_->getNode ()));
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  localIndices_.resize (MaxNumEntries_);
181  Values_.resize (MaxNumEntries_);
182 
183  // now compute:
184  // - the number of nonzero per row
185  // - the total number of nonzeros
186  // - the diagonal entries
187 
188  // compute nonzeros (total and per-row), and store the
189  // diagonal entries (already modified)
190  size_t ActualMaxNumEntries = 0;
191 
192  for (size_t i = 0; i < numRows; ++i) {
193  NumEntries_[i] = 0;
194  size_t Nnz, NewNnz = 0;
195  A_->getLocalRowCopy (i, localIndices_, Values_, Nnz);
196  for (size_t j = 0; j < Nnz; ++j) {
197  // FIXME (mfh 03 Apr 2013) This assumes the following:
198  //
199  // 1. Row Map, range Map, and domain Map are all the same.
200  //
201  // 2. The column Map's list of GIDs on this process is the
202  // domain Map's list of GIDs, followed by remote GIDs. Thus,
203  // for any GID in the domain Map on this process, its LID in
204  // the domain Map (and therefore in the row Map, by (1)) is
205  // the same as its LID in the column Map. (Hence the
206  // less-than test, which if true, means that localIndices_[j]
207  // belongs to the row Map.)
208  if (static_cast<size_t> (localIndices_[j]) < numRows) {
209  ++NewNnz;
210  }
211  }
212 
213  if (NewNnz > ActualMaxNumEntries) {
214  ActualMaxNumEntries = NewNnz;
215  }
216 
217  NumNonzeros_ += NewNnz;
218  NumEntries_[i] = NewNnz;
219  }
220 
221  MaxNumEntries_ = ActualMaxNumEntries;
222 }
223 
224 
225 template<class MatrixType>
227 {}
228 
229 
230 template<class MatrixType>
233 {
234  return localRowMap_->getComm ();
235 }
236 
237 
238 template<class MatrixType>
241 {
242  return A_->getNode ();
243 }
244 
245 
246 template<class MatrixType>
247 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
248  typename MatrixType::global_ordinal_type,
249  typename MatrixType::node_type> >
251 {
252  return localRowMap_;
253 }
254 
255 
256 template<class MatrixType>
257 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
258  typename MatrixType::global_ordinal_type,
259  typename MatrixType::node_type> >
261 {
262  return localRowMap_; // FIXME (mfh 20 Nov 2013)
263 }
264 
265 
266 template<class MatrixType>
267 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
268  typename MatrixType::global_ordinal_type,
269  typename MatrixType::node_type> >
271 {
272  return localDomainMap_;
273 }
274 
275 
276 template<class MatrixType>
277 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
278  typename MatrixType::global_ordinal_type,
279  typename MatrixType::node_type> >
281 {
282  return localRangeMap_;
283 }
284 
285 
286 template<class MatrixType>
287 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type,
288  typename MatrixType::global_ordinal_type,
289  typename MatrixType::node_type> >
291 {
292  // FIXME (mfh 20 Nov 2013) This is not what the documentation says
293  // this method should do! It should return the graph of the locally
294  // filtered matrix, not the original matrix's graph.
295  return A_->getGraph ();
296 }
297 
298 
299 template<class MatrixType>
301 {
302  return static_cast<global_size_t> (localRangeMap_->getNodeNumElements ());
303 }
304 
305 
306 template<class MatrixType>
308 {
309  return static_cast<global_size_t> (localDomainMap_->getNodeNumElements ());
310 }
311 
312 
313 template<class MatrixType>
315 {
316  return static_cast<size_t> (localRangeMap_->getNodeNumElements ());
317 }
318 
319 
320 template<class MatrixType>
322 {
323  return static_cast<size_t> (localDomainMap_->getNodeNumElements ());
324 }
325 
326 
327 template<class MatrixType>
328 typename MatrixType::global_ordinal_type
330 {
331  return A_->getIndexBase ();
332 }
333 
334 
335 template<class MatrixType>
337 {
338  return NumNonzeros_;
339 }
340 
341 
342 template<class MatrixType>
344 {
345  return NumNonzeros_;
346 }
347 
348 
349 template<class MatrixType>
350 size_t
353 {
354  const local_ordinal_type localRow = getRowMap ()->getLocalElement (globalRow);
356  // NOTE (mfh 26 Mar 2014) We return zero if globalRow is not in
357  // the row Map on this process, since "get the number of entries
358  // in the global row" refers only to what the calling process owns
359  // in that row. In this case, it owns no entries in that row,
360  // since it doesn't own the row.
361  return 0;
362  } else {
363  return NumEntries_[localRow];
364  }
365 }
366 
367 
368 template<class MatrixType>
369 size_t
372 {
373  // FIXME (mfh 07 Jul 2014) Shouldn't localRow be a local row index
374  // in the matrix's row Map, not in the LocalFilter's row Map? The
375  // latter is different; it even has different global indices!
376  // (Maybe _that_'s the bug.)
377 
378  if (getRowMap ()->isNodeLocalElement (localRow)) {
379  return NumEntries_[localRow];
380  } else {
381  // NOTE (mfh 26 Mar 2014) We return zero if localRow is not in the
382  // row Map on this process, since "get the number of entries in
383  // the local row" refers only to what the calling process owns in
384  // that row. In this case, it owns no entries in that row, since
385  // it doesn't own the row.
386  return 0;
387  }
388 }
389 
390 
391 template<class MatrixType>
393 {
394  return A_->getGlobalNumDiags ();
395 }
396 
397 
398 template<class MatrixType>
400 {
401  return A_->getNodeNumDiags ();
402 }
403 
404 
405 template<class MatrixType>
407 {
408  return MaxNumEntries_;
409 }
410 
411 
412 template<class MatrixType>
414 {
415  return MaxNumEntries_;
416 }
417 
418 
419 template<class MatrixType>
421 {
422  return true;
423 }
424 
425 
426 template<class MatrixType>
428 {
429  return A_->isLowerTriangular();
430 }
431 
432 
433 template<class MatrixType>
435 {
436  return A_->isUpperTriangular();
437 }
438 
439 
440 template<class MatrixType>
442 {
443  return A_->isLocallyIndexed ();
444 }
445 
446 
447 template<class MatrixType>
449 {
450  return A_->isGloballyIndexed();
451 }
452 
453 
454 template<class MatrixType>
456 {
457  return A_->isFillComplete ();
458 }
459 
460 
461 template<class MatrixType>
462 void
465  const Teuchos::ArrayView<global_ordinal_type>& globalIndices,
466  const Teuchos::ArrayView<scalar_type>& values,
467  size_t& numEntries) const
468 {
469  typedef local_ordinal_type LO;
470  typedef typename Teuchos::Array<LO>::size_type size_type;
471 
472  const LO localRow = this->getRowMap ()->getLocalElement (globalRow);
473  if (localRow == Teuchos::OrdinalTraits<LO>::invalid ()) {
474  // NOTE (mfh 26 Mar 2014) We return no entries if globalRow is not
475  // in the row Map on this process, since "get a copy of the
476  // entries in the global row" refers only to what the calling
477  // process owns in that row. In this case, it owns no entries in
478  // that row, since it doesn't own the row.
479  numEntries = 0;
480  }
481  else {
482  // First get a copy of the current row using local indices. Then,
483  // convert to global indices using the input matrix's column Map.
484  //
485  numEntries = this->getNumEntriesInLocalRow (localRow);
486  // FIXME (mfh 26 Mar 2014) If local_ordinal_type ==
487  // global_ordinal_type, we could just alias the input array
488  // instead of allocating a temporary array.
489  Teuchos::Array<LO> localIndices (numEntries);
490  this->getLocalRowCopy (localRow, localIndices (), values, numEntries);
491 
492  const map_type& colMap = * (this->getColMap ());
493 
494  // Don't fill the output array beyond its size.
495  const size_type numEnt =
496  std::min (static_cast<size_type> (numEntries),
497  std::min (globalIndices.size (), values.size ()));
498  for (size_type k = 0; k < numEnt; ++k) {
499  globalIndices[k] = colMap.getGlobalElement (localIndices[k]);
500  }
501  }
502 }
503 
504 
505 template<class MatrixType>
506 void
510  const Teuchos::ArrayView<scalar_type> &Values,
511  size_t &NumEntries) const
512 {
513  typedef local_ordinal_type LO;
514  typedef global_ordinal_type GO;
515 
516  if (! A_->getRowMap ()->isNodeLocalElement (LocalRow)) {
517  // The calling process owns zero entries in the row.
518  NumEntries = 0;
519  return;
520  }
521 
522  if (A_->getRowMap()->getComm()->getSize() == 1) {
523  A_->getLocalRowCopy (LocalRow, Indices (), Values (), NumEntries);
524  return;
525  }
526 
527 
528  const size_t numEntInLclRow = NumEntries_[LocalRow];
529  if (static_cast<size_t> (Indices.size ()) < numEntInLclRow ||
530  static_cast<size_t> (Values.size ()) < numEntInLclRow) {
531  // FIXME (mfh 07 Jul 2014) Return an error code instead of
532  // throwing. We should really attempt to fill as much space as
533  // we're given, in this case.
535  true, std::runtime_error,
536  "Ifpack2::LocalFilter::getLocalRowCopy: Invalid output array length. "
537  "The output arrays must each have length at least " << numEntInLclRow
538  << " for local row " << LocalRow << " on Process "
539  << localRowMap_->getComm ()->getRank () << ".");
540  }
541  else if (numEntInLclRow == static_cast<size_t> (0)) {
542  // getNumEntriesInLocalRow() returns zero if LocalRow is not owned
543  // by the calling process. In that case, the calling process owns
544  // zero entries in the row.
545  NumEntries = 0;
546  return;
547  }
548 
549  // Always extract using the temporary arrays Values_ and
550  // localIndices_. This is because I may need more space than that
551  // given by the user. The users expects only the local (in the
552  // domain Map) column indices, but I have to extract both local and
553  // remote (not in the domain Map) column indices.
554  //
555  // FIXME (mfh 07 Jul 2014) Use of temporary local storage is not
556  // conducive to thread parallelism. A better way would be to change
557  // the interface so that it only extracts values for the "local"
558  // column indices. CrsMatrix could take a set of column indices,
559  // and return their corresponding values.
560  size_t numEntInMat = 0;
561  A_->getLocalRowCopy (LocalRow, localIndices_ (), Values_ (), numEntInMat);
562 
563  // Fill the user's arrays with the "local" indices and values in
564  // that row. Note that the matrix might have a different column Map
565  // than the local filter.
566  const map_type& matrixDomMap = * (A_->getDomainMap ());
567  const map_type& matrixColMap = * (A_->getColMap ());
568 
569  const size_t capacity = static_cast<size_t> (std::min (Indices.size (),
570  Values.size ()));
571  NumEntries = 0;
572  const size_t numRows = localRowMap_->getNodeNumElements (); // superfluous
573  const bool buggy = true; // mfh 07 Jul 2014: See FIXME below.
574  for (size_t j = 0; j < numEntInMat; ++j) {
575  // The LocalFilter only includes entries in the domain Map on
576  // the calling process. We figure out whether an entry is in
577  // the domain Map by converting the (matrix column Map) index to
578  // a global index, and then asking whether that global index is
579  // in the domain Map.
580  const LO matrixLclCol = localIndices_[j];
581  const GO gblCol = matrixColMap.getGlobalElement (matrixLclCol);
582 
583  // FIXME (mfh 07 Jul 2014) This is the likely center of Bug 5992
584  // and perhaps other bugs, like Bug 6117. If 'buggy' is true,
585  // Ifpack2 tests pass; if 'buggy' is false, the tests don't pass.
586  // This suggests that Ifpack2 classes could be using LocalFilter
587  // incorrectly, perhaps by giving it an incorrect domain Map.
588  if (buggy) {
589  // only local indices
590  if ((size_t) localIndices_[j] < numRows) {
591  Indices[NumEntries] = localIndices_[j];
592  Values[NumEntries] = Values_[j];
593  NumEntries++;
594  }
595  } else {
596  if (matrixDomMap.isNodeGlobalElement (gblCol)) {
597  // Don't fill more space than the user gave us. It's an error
598  // for them not to give us enough space, but we still shouldn't
599  // overwrite memory that doesn't belong to us.
600  if (NumEntries < capacity) {
601  Indices[NumEntries] = matrixLclCol;
602  Values[NumEntries] = Values_[j];
603  }
604  NumEntries++;
605  }
606  }
607  }
608 }
609 
610 
611 template<class MatrixType>
612 void
617 {
618  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
619  "Ifpack2::LocalFilter does not implement getGlobalRowView.");
620 }
621 
622 
623 template<class MatrixType>
624 void
629 {
630  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
631  "Ifpack2::LocalFilter does not implement getLocalRowView.");
632 }
633 
634 
635 template<class MatrixType>
636 void
638 getLocalDiagCopy (Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& diag) const
639 {
640  using Teuchos::RCP;
641  typedef Tpetra::Vector<scalar_type, local_ordinal_type,
642  global_ordinal_type, node_type> vector_type;
643  // This is always correct, and doesn't require a collective check
644  // for sameness of Maps.
645  RCP<vector_type> diagView = diag.offsetViewNonConst (A_->getRowMap (), 0);
646  A_->getLocalDiagCopy (*diagView);
647 }
648 
649 
650 template<class MatrixType>
651 void
653 leftScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
654 {
655  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
656  "Ifpack2::LocalFilter does not implement leftScale.");
657 }
658 
659 
660 template<class MatrixType>
661 void
663 rightScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
664 {
665  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
666  "Ifpack2::LocalFilter does not implement rightScale.");
667 }
668 
669 
670 template<class MatrixType>
671 void
673 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
674  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
675  Teuchos::ETransp mode,
676  scalar_type alpha,
677  scalar_type beta) const
678 {
679  typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV;
681  X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
682  "Ifpack2::LocalFilter::apply: X and Y must have the same number of columns. "
683  "X has " << X.getNumVectors () << " columns, but Y has "
684  << Y.getNumVectors () << " columns.");
685 
686 #ifdef HAVE_IFPACK2_DEBUG
687  {
689  Teuchos::Array<magnitude_type> norms (X.getNumVectors ());
690  X.norm1 (norms ());
691  bool good = true;
692  for (size_t j = 0; j < X.getNumVectors (); ++j) {
693  if (STM::isnaninf (norms[j])) {
694  good = false;
695  break;
696  }
697  }
698  TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error, "Ifpack2::LocalFilter::apply: The 1-norm of the input X is NaN or Inf.");
699  }
700 #endif // HAVE_IFPACK2_DEBUG
701 
702  if (&X == &Y) {
703  // FIXME (mfh 23 Apr 2014) We have to do more work to figure out
704  // if X and Y alias one another. For example, we should check
705  // whether one is a noncontiguous view of the other.
706  //
707  // X and Y alias one another, so we have to copy X.
708  MV X_copy (X, Teuchos::Copy);
709  applyNonAliased (X_copy, Y, mode, alpha, beta);
710  } else {
711  applyNonAliased (X, Y, mode, alpha, beta);
712  }
713 
714 #ifdef HAVE_IFPACK2_DEBUG
715  {
717  Teuchos::Array<magnitude_type> norms (Y.getNumVectors ());
718  Y.norm1 (norms ());
719  bool good = true;
720  for (size_t j = 0; j < Y.getNumVectors (); ++j) {
721  if (STM::isnaninf (norms[j])) {
722  good = false;
723  break;
724  }
725  }
726  TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error, "Ifpack2::LocalFilter::apply: The 1-norm of the output Y is NaN or Inf.");
727  }
728 #endif // HAVE_IFPACK2_DEBUG
729 }
730 
731 template<class MatrixType>
732 void
734 applyNonAliased (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
735  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
736  Teuchos::ETransp mode,
737  scalar_type alpha,
738  scalar_type beta) const
739 {
740  using Teuchos::ArrayView;
741  using Teuchos::ArrayRCP;
743 
744  const scalar_type zero = STS::zero ();
745  const scalar_type one = STS::one ();
746 
747  if (beta == zero) {
748  Y.putScalar (zero);
749  }
750  else if (beta != one) {
751  Y.scale (beta);
752  }
753 
754  const size_t NumVectors = Y.getNumVectors ();
755  const size_t numRows = localRowMap_->getNodeNumElements ();
756 
757  // FIXME (mfh 14 Apr 2014) At some point, we would like to
758  // parallelize this using Kokkos. This would require a
759  // Kokkos-friendly version of getLocalRowCopy, perhaps with
760  // thread-local storage.
761 
762  const bool constantStride = X.isConstantStride () && Y.isConstantStride ();
763  if (constantStride) {
764  // Since both X and Y have constant stride, we can work with "1-D"
765  // views of their data.
766  const size_t x_stride = X.getStride();
767  const size_t y_stride = Y.getStride();
768  ArrayRCP<scalar_type> y_rcp = Y.get1dViewNonConst();
769  ArrayRCP<const scalar_type> x_rcp = X.get1dView();
770  ArrayView<scalar_type> y_ptr = y_rcp();
771  ArrayView<const scalar_type> x_ptr = x_rcp();
772  for (size_t i = 0; i < numRows; ++i) {
773  size_t Nnz;
774  // Use this class's getrow to make the below code simpler
775  getLocalRowCopy (i, localIndices_ (), Values_ (), Nnz);
776  if (mode == Teuchos::NO_TRANS) {
777  for (size_t j = 0; j < Nnz; ++j) {
778  const local_ordinal_type col = localIndices_[j];
779  for (size_t k = 0; k < NumVectors; ++k) {
780  y_ptr[i + y_stride*k] +=
781  alpha * Values_[j] * x_ptr[col + x_stride*k];
782  }
783  }
784  }
785  else if (mode == Teuchos::TRANS) {
786  for (size_t j = 0; j < Nnz; ++j) {
787  const local_ordinal_type col = localIndices_[j];
788  for (size_t k = 0; k < NumVectors; ++k) {
789  y_ptr[col + y_stride*k] +=
790  alpha * Values_[j] * x_ptr[i + x_stride*k];
791  }
792  }
793  }
794  else { //mode==Teuchos::CONJ_TRANS
795  for (size_t j = 0; j < Nnz; ++j) {
796  const local_ordinal_type col = localIndices_[j];
797  for (size_t k = 0; k < NumVectors; ++k) {
798  y_ptr[col + y_stride*k] +=
799  alpha * STS::conjugate (Values_[j]) * x_ptr[i + x_stride*k];
800  }
801  }
802  }
803  }
804  }
805  else {
806  // At least one of X or Y does not have constant stride.
807  // This means we must work with "2-D" views of their data.
808  ArrayRCP<ArrayRCP<const scalar_type> > x_ptr = X.get2dView();
809  ArrayRCP<ArrayRCP<scalar_type> > y_ptr = Y.get2dViewNonConst();
810 
811  for (size_t i = 0; i < numRows; ++i) {
812  size_t Nnz;
813  // Use this class's getrow to make the below code simpler
814  getLocalRowCopy (i, localIndices_ (), Values_ (), Nnz);
815  if (mode == Teuchos::NO_TRANS) {
816  for (size_t k = 0; k < NumVectors; ++k) {
817  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
818  ArrayView<scalar_type> y_local = (y_ptr())[k]();
819  for (size_t j = 0; j < Nnz; ++j) {
820  y_local[i] += alpha * Values_[j] * x_local[localIndices_[j]];
821  }
822  }
823  }
824  else if (mode == Teuchos::TRANS) {
825  for (size_t k = 0; k < NumVectors; ++k) {
826  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
827  ArrayView<scalar_type> y_local = (y_ptr())[k]();
828  for (size_t j = 0; j < Nnz; ++j) {
829  y_local[localIndices_[j]] += alpha * Values_[j] * x_local[i];
830  }
831  }
832  }
833  else { //mode==Teuchos::CONJ_TRANS
834  for (size_t k = 0; k < NumVectors; ++k) {
835  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
836  ArrayView<scalar_type> y_local = (y_ptr())[k]();
837  for (size_t j = 0; j < Nnz; ++j) {
838  y_local[localIndices_[j]] +=
839  alpha * STS::conjugate (Values_[j]) * x_local[i];
840  }
841  }
842  }
843  }
844  }
845 }
846 
847 template<class MatrixType>
849 {
850  return true;
851 }
852 
853 
854 template<class MatrixType>
856 {
857  return false;
858 }
859 
860 
861 template<class MatrixType>
862 typename
863 LocalFilter<MatrixType>::mag_type
865 {
866 #ifdef TPETRA_HAVE_KOKKOS_REFACTOR
867  typedef Kokkos::Details::ArithTraits<scalar_type> STS;
868  typedef Kokkos::Details::ArithTraits<mag_type> STM;
869 #else
872 #endif
873  typedef typename Teuchos::Array<scalar_type>::size_type size_type;
874 
875  const size_type maxNumRowEnt = getNodeMaxNumRowEntries ();
876  Teuchos::Array<local_ordinal_type> ind (maxNumRowEnt);
877  Teuchos::Array<scalar_type> val (maxNumRowEnt);
878  const size_t numRows = static_cast<size_t> (localRowMap_->getNodeNumElements ());
879 
880  // FIXME (mfh 03 Apr 2013) Scale during sum to avoid overflow.
881  mag_type sumSquared = STM::zero ();
882  for (size_t i = 0; i < numRows; ++i) {
883  size_t numEntries = 0;
884  this->getLocalRowCopy (i, ind (), val (), numEntries);
885  for (size_type k = 0; k < static_cast<size_type> (numEntries); ++k) {
886  const mag_type v_k_abs = STS::magnitude (val[k]);
887  sumSquared += v_k_abs * v_k_abs;
888  }
889  }
890  return STM::squareroot (sumSquared); // Different for each process; that's OK.
891 }
892 
893 template<class MatrixType>
894 std::string
896 {
898  std::ostringstream os;
899 
900  os << "Ifpack2::LocalFilter: {";
901  os << "MatrixType: " << TypeNameTraits<MatrixType>::name ();
902  if (this->getObjectLabel () != "") {
903  os << ", Label: \"" << this->getObjectLabel () << "\"";
904  }
905  os << ", Number of rows: " << getGlobalNumRows ()
906  << ", Number of columns: " << getGlobalNumCols ()
907  << "}";
908  return os.str ();
909 }
910 
911 
912 template<class MatrixType>
913 void
916  const Teuchos::EVerbosityLevel verbLevel) const
917 {
918  using Teuchos::OSTab;
920  using std::endl;
921 
922  const Teuchos::EVerbosityLevel vl =
923  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
924 
925  if (vl > Teuchos::VERB_NONE) {
926  // describe() starts with a tab, by convention.
927  OSTab tab0 (out);
928 
929  out << "Ifpack2::LocalFilter:" << endl;
930  OSTab tab1 (out);
931  out << "MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
932  if (this->getObjectLabel () != "") {
933  out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
934  }
935  out << "Number of rows: " << getGlobalNumRows () << endl
936  << "Number of columns: " << getGlobalNumCols () << endl
937  << "Number of nonzeros: " << NumNonzeros_ << endl;
938 
939  if (vl > Teuchos::VERB_LOW) {
940  out << "Row Map:" << endl;
941  localRowMap_->describe (out, vl);
942  out << "Domain Map:" << endl;
943  localDomainMap_->describe (out, vl);
944  out << "Range Map:" << endl;
945  localRangeMap_->describe (out, vl);
946  }
947  }
948 }
949 
950 template<class MatrixType>
951 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
952  typename MatrixType::local_ordinal_type,
953  typename MatrixType::global_ordinal_type,
954  typename MatrixType::node_type> >
956 {
957  return A_;
958 }
959 
960 
961 } // namespace Ifpack2
962 
963 #define IFPACK2_LOCALFILTER_INSTANT(S,LO,GO,N) \
964  template class Ifpack2::LocalFilter< Tpetra::RowMatrix<S, LO, GO, N> >;
965 
966 #endif
virtual void getLocalRowView(local_ordinal_type LocalRow, Teuchos::ArrayView< const local_ordinal_type > &indices, Teuchos::ArrayView< const scalar_type > &values) const
Extract a const, non-persisting view of local indices in a specified row of the matrix.
Definition: Ifpack2_LocalFilter_def.hpp:626
virtual bool isFillComplete() const
Returns true if fillComplete() has been called.
Definition: Ifpack2_LocalFilter_def.hpp:455
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:663
virtual mag_type getFrobeniusNorm() const
The Frobenius norm of the (locally filtered) matrix.
Definition: Ifpack2_LocalFilter_def.hpp:864
VERB_DEFAULT
basic_OSTab< char > OSTab
virtual Teuchos::RCP< const row_matrix_type > getUnderlyingMatrix() const
Return matrix that LocalFilter was built on.
Definition: Ifpack2_LocalFilter_def.hpp:955
VERB_LOW
VERB_NONE
virtual global_size_t getGlobalNumCols() const
The number of global columns in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:307
virtual bool isUpperTriangular() const
Whether this matrix is upper triangular.
Definition: Ifpack2_LocalFilter_def.hpp:434
virtual bool isLocallyIndexed() const
Whether the underlying sparse matrix is locally (opposite of globally) indexed.
Definition: Ifpack2_LocalFilter_def.hpp:441
virtual size_t getNodeNumDiags() const
The number of local diagonal entries, based on global row/column index comparisons.
Definition: Ifpack2_LocalFilter_def.hpp:399
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:915
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type size() const
virtual size_t getNodeNumCols() const
The number of columns in the (locally filtered) matrix.
Definition: Ifpack2_LocalFilter_def.hpp:321
virtual void getGlobalRowView(global_ordinal_type GlobalRow, Teuchos::ArrayView< const global_ordinal_type > &indices, Teuchos::ArrayView< const scalar_type > &values) const
Extract a const, non-persisting view of global indices in a specified row of the matrix.
Definition: Ifpack2_LocalFilter_def.hpp:614
virtual global_size_t getGlobalNumDiags() const
The number of global diagonal entries, based on global row/column index comparisons.
Definition: Ifpack2_LocalFilter_def.hpp:392
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_LocalFilter_decl.hpp:190
virtual global_size_t getGlobalNumRows() const
The number of global rows in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:300
virtual Teuchos::RCP< const map_type > getRangeMap() const
Returns the Map that describes the range distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:280
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:207
virtual bool supportsRowViews() const
Returns true if RowViews are supported.
Definition: Ifpack2_LocalFilter_def.hpp:855
Ordinal size_type
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_LocalFilter_decl.hpp:187
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual Teuchos::RCP< const map_type > getColMap() const
Returns the Map that describes the column distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:260
virtual void getGlobalRowCopy(global_ordinal_type GlobalRow, const Teuchos::ArrayView< global_ordinal_type > &Indices, const Teuchos::ArrayView< scalar_type > &Values, size_t &NumEntries) const
Get the entries in the given row, using global indices.
Definition: Ifpack2_LocalFilter_def.hpp:464
virtual void getLocalRowCopy(local_ordinal_type LocalRow, const Teuchos::ArrayView< local_ordinal_type > &Indices, const Teuchos::ArrayView< scalar_type > &Values, size_t &NumEntries) const
Get the entries in the given row, using local indices.
Definition: Ifpack2_LocalFilter_def.hpp:508
virtual bool isLowerTriangular() const
Whether this matrix is lower triangular.
Definition: Ifpack2_LocalFilter_def.hpp:427
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:653
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:290
virtual size_t getNodeNumRows() const
The number of rows owned on the calling process.
Definition: Ifpack2_LocalFilter_def.hpp:314
LinearOp zero(const VectorSpace &vs)
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:673
virtual global_ordinal_type getIndexBase() const
Returns the index base for global indices for this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:329
virtual bool hasColMap() const
Whether this matrix has a well-defined column Map.
Definition: Ifpack2_LocalFilter_def.hpp:420
void resize(size_type new_size, const value_type &x=value_type())
virtual size_t getGlobalMaxNumRowEntries() const
The maximum number of entries across all rows/columns on all processes.
Definition: Ifpack2_LocalFilter_def.hpp:406
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the communicator.
Definition: Ifpack2_LocalFilter_def.hpp:232
virtual std::string description() const
A one-line description of this object.
Definition: Ifpack2_LocalFilter_def.hpp:895
virtual size_t getNodeMaxNumRowEntries() const
The maximum number of entries across all rows/columns on this process.
Definition: Ifpack2_LocalFilter_def.hpp:413
LocalFilter(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_LocalFilter_def.hpp:88
virtual ~LocalFilter()
Destructor.
Definition: Ifpack2_LocalFilter_def.hpp:226
virtual global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:336
virtual Teuchos::RCP< const map_type > getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:250
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_LocalFilter_decl.hpp:184
virtual size_t getNodeNumEntries() const
Returns the local number of entries in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:343
virtual bool isGloballyIndexed() const
Whether the underlying sparse matrix is globally (opposite of locally) indexed.
Definition: Ifpack2_LocalFilter_def.hpp:448
T * getRawPtr() const
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:638
virtual Teuchos::RCP< node_type > getNode() const
Returns the underlying Node object.
Definition: Ifpack2_LocalFilter_def.hpp:240
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:270
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:352
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
virtual bool hasTransposeApply() const
Whether this operator supports applying the transpose or conjugate transpose.
Definition: Ifpack2_LocalFilter_def.hpp:848
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:371
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_LocalFilter_decl.hpp:193