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
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>
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
352 getNumEntriesInGlobalRow (global_ordinal_type globalRow) const
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 MaxNumEntries_;
395 }
396 
397 
398 template<class MatrixType>
400 {
401  return MaxNumEntries_;
402 }
403 
404 
405 template<class MatrixType>
407 {
408  return true;
409 }
410 
411 
412 template<class MatrixType>
414 {
415  return A_->isLocallyIndexed ();
416 }
417 
418 
419 template<class MatrixType>
421 {
422  return A_->isGloballyIndexed();
423 }
424 
425 
426 template<class MatrixType>
428 {
429  return A_->isFillComplete ();
430 }
431 
432 
433 template<class MatrixType>
434 void
436 getGlobalRowCopy (global_ordinal_type globalRow,
437  const Teuchos::ArrayView<global_ordinal_type>& globalIndices,
438  const Teuchos::ArrayView<scalar_type>& values,
439  size_t& numEntries) const
440 {
441  typedef local_ordinal_type LO;
442  typedef typename Teuchos::Array<LO>::size_type size_type;
443 
444  const LO localRow = this->getRowMap ()->getLocalElement (globalRow);
445  if (localRow == Teuchos::OrdinalTraits<LO>::invalid ()) {
446  // NOTE (mfh 26 Mar 2014) We return no entries if globalRow is not
447  // in the row Map on this process, since "get a copy of the
448  // entries in the global row" refers only to what the calling
449  // process owns in that row. In this case, it owns no entries in
450  // that row, since it doesn't own the row.
451  numEntries = 0;
452  }
453  else {
454  // First get a copy of the current row using local indices. Then,
455  // convert to global indices using the input matrix's column Map.
456  //
457  numEntries = this->getNumEntriesInLocalRow (localRow);
458  // FIXME (mfh 26 Mar 2014) If local_ordinal_type ==
459  // global_ordinal_type, we could just alias the input array
460  // instead of allocating a temporary array.
461  Teuchos::Array<LO> localIndices (numEntries);
462  this->getLocalRowCopy (localRow, localIndices (), values, numEntries);
463 
464  const map_type& colMap = * (this->getColMap ());
465 
466  // Don't fill the output array beyond its size.
467  const size_type numEnt =
468  std::min (static_cast<size_type> (numEntries),
469  std::min (globalIndices.size (), values.size ()));
470  for (size_type k = 0; k < numEnt; ++k) {
471  globalIndices[k] = colMap.getGlobalElement (localIndices[k]);
472  }
473  }
474 }
475 
476 
477 template<class MatrixType>
478 void
482  const Teuchos::ArrayView<scalar_type> &Values,
483  size_t &NumEntries) const
484 {
485  typedef local_ordinal_type LO;
486  typedef global_ordinal_type GO;
487 
488  if (! A_->getRowMap ()->isNodeLocalElement (LocalRow)) {
489  // The calling process owns zero entries in the row.
490  NumEntries = 0;
491  return;
492  }
493 
494  if (A_->getRowMap()->getComm()->getSize() == 1) {
495  A_->getLocalRowCopy (LocalRow, Indices (), Values (), NumEntries);
496  return;
497  }
498 
499 
500  const size_t numEntInLclRow = NumEntries_[LocalRow];
501  if (static_cast<size_t> (Indices.size ()) < numEntInLclRow ||
502  static_cast<size_t> (Values.size ()) < numEntInLclRow) {
503  // FIXME (mfh 07 Jul 2014) Return an error code instead of
504  // throwing. We should really attempt to fill as much space as
505  // we're given, in this case.
507  true, std::runtime_error,
508  "Ifpack2::LocalFilter::getLocalRowCopy: Invalid output array length. "
509  "The output arrays must each have length at least " << numEntInLclRow
510  << " for local row " << LocalRow << " on Process "
511  << localRowMap_->getComm ()->getRank () << ".");
512  }
513  else if (numEntInLclRow == static_cast<size_t> (0)) {
514  // getNumEntriesInLocalRow() returns zero if LocalRow is not owned
515  // by the calling process. In that case, the calling process owns
516  // zero entries in the row.
517  NumEntries = 0;
518  return;
519  }
520 
521  // Always extract using the temporary arrays Values_ and
522  // localIndices_. This is because I may need more space than that
523  // given by the user. The users expects only the local (in the
524  // domain Map) column indices, but I have to extract both local and
525  // remote (not in the domain Map) column indices.
526  //
527  // FIXME (mfh 07 Jul 2014) Use of temporary local storage is not
528  // conducive to thread parallelism. A better way would be to change
529  // the interface so that it only extracts values for the "local"
530  // column indices. CrsMatrix could take a set of column indices,
531  // and return their corresponding values.
532  size_t numEntInMat = 0;
533  A_->getLocalRowCopy (LocalRow, localIndices_ (), Values_ (), numEntInMat);
534 
535  // Fill the user's arrays with the "local" indices and values in
536  // that row. Note that the matrix might have a different column Map
537  // than the local filter.
538  const map_type& matrixDomMap = * (A_->getDomainMap ());
539  const map_type& matrixColMap = * (A_->getColMap ());
540 
541  const size_t capacity = static_cast<size_t> (std::min (Indices.size (),
542  Values.size ()));
543  NumEntries = 0;
544  const size_t numRows = localRowMap_->getNodeNumElements (); // superfluous
545  const bool buggy = true; // mfh 07 Jul 2014: See FIXME below.
546  for (size_t j = 0; j < numEntInMat; ++j) {
547  // The LocalFilter only includes entries in the domain Map on
548  // the calling process. We figure out whether an entry is in
549  // the domain Map by converting the (matrix column Map) index to
550  // a global index, and then asking whether that global index is
551  // in the domain Map.
552  const LO matrixLclCol = localIndices_[j];
553  const GO gblCol = matrixColMap.getGlobalElement (matrixLclCol);
554 
555  // FIXME (mfh 07 Jul 2014) This is the likely center of Bug 5992
556  // and perhaps other bugs, like Bug 6117. If 'buggy' is true,
557  // Ifpack2 tests pass; if 'buggy' is false, the tests don't pass.
558  // This suggests that Ifpack2 classes could be using LocalFilter
559  // incorrectly, perhaps by giving it an incorrect domain Map.
560  if (buggy) {
561  // only local indices
562  if ((size_t) localIndices_[j] < numRows) {
563  Indices[NumEntries] = localIndices_[j];
564  Values[NumEntries] = Values_[j];
565  NumEntries++;
566  }
567  } else {
568  if (matrixDomMap.isNodeGlobalElement (gblCol)) {
569  // Don't fill more space than the user gave us. It's an error
570  // for them not to give us enough space, but we still shouldn't
571  // overwrite memory that doesn't belong to us.
572  if (NumEntries < capacity) {
573  Indices[NumEntries] = matrixLclCol;
574  Values[NumEntries] = Values_[j];
575  }
576  NumEntries++;
577  }
578  }
579  }
580 }
581 
582 
583 template<class MatrixType>
584 void
586 getGlobalRowView (global_ordinal_type /* GlobalRow */,
588  Teuchos::ArrayView<const scalar_type> &/* values */) const
589 {
590  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
591  "Ifpack2::LocalFilter does not implement getGlobalRowView.");
592 }
593 
594 
595 template<class MatrixType>
596 void
600  Teuchos::ArrayView<const scalar_type> &/* values */) const
601 {
602  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
603  "Ifpack2::LocalFilter does not implement getLocalRowView.");
604 }
605 
606 
607 template<class MatrixType>
608 void
610 getLocalDiagCopy (Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& diag) const
611 {
612  using Teuchos::RCP;
613  typedef Tpetra::Vector<scalar_type, local_ordinal_type,
614  global_ordinal_type, node_type> vector_type;
615  // This is always correct, and doesn't require a collective check
616  // for sameness of Maps.
617  RCP<vector_type> diagView = diag.offsetViewNonConst (A_->getRowMap (), 0);
618  A_->getLocalDiagCopy (*diagView);
619 }
620 
621 
622 template<class MatrixType>
623 void
625 leftScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& /* x */)
626 {
627  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
628  "Ifpack2::LocalFilter does not implement leftScale.");
629 }
630 
631 
632 template<class MatrixType>
633 void
635 rightScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& /* x */)
636 {
637  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
638  "Ifpack2::LocalFilter does not implement rightScale.");
639 }
640 
641 
642 template<class MatrixType>
643 void
645 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
646  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
647  Teuchos::ETransp mode,
648  scalar_type alpha,
649  scalar_type beta) const
650 {
651  typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV;
653  X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
654  "Ifpack2::LocalFilter::apply: X and Y must have the same number of columns. "
655  "X has " << X.getNumVectors () << " columns, but Y has "
656  << Y.getNumVectors () << " columns.");
657 
658 #ifdef HAVE_IFPACK2_DEBUG
659  {
661  Teuchos::Array<magnitude_type> norms (X.getNumVectors ());
662  X.norm1 (norms ());
663  bool good = true;
664  for (size_t j = 0; j < X.getNumVectors (); ++j) {
665  if (STM::isnaninf (norms[j])) {
666  good = false;
667  break;
668  }
669  }
670  TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error, "Ifpack2::LocalFilter::apply: The 1-norm of the input X is NaN or Inf.");
671  }
672 #endif // HAVE_IFPACK2_DEBUG
673 
674  if (&X == &Y) {
675  // FIXME (mfh 23 Apr 2014) We have to do more work to figure out
676  // if X and Y alias one another. For example, we should check
677  // whether one is a noncontiguous view of the other.
678  //
679  // X and Y alias one another, so we have to copy X.
680  MV X_copy (X, Teuchos::Copy);
681  applyNonAliased (X_copy, Y, mode, alpha, beta);
682  } else {
683  applyNonAliased (X, Y, mode, alpha, beta);
684  }
685 
686 #ifdef HAVE_IFPACK2_DEBUG
687  {
689  Teuchos::Array<magnitude_type> norms (Y.getNumVectors ());
690  Y.norm1 (norms ());
691  bool good = true;
692  for (size_t j = 0; j < Y.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 output Y is NaN or Inf.");
699  }
700 #endif // HAVE_IFPACK2_DEBUG
701 }
702 
703 template<class MatrixType>
704 void
706 applyNonAliased (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
707  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
708  Teuchos::ETransp mode,
709  scalar_type alpha,
710  scalar_type beta) const
711 {
712  using Teuchos::ArrayView;
713  using Teuchos::ArrayRCP;
715 
716  const scalar_type zero = STS::zero ();
717  const scalar_type one = STS::one ();
718 
719  if (beta == zero) {
720  Y.putScalar (zero);
721  }
722  else if (beta != one) {
723  Y.scale (beta);
724  }
725 
726  const size_t NumVectors = Y.getNumVectors ();
727  const size_t numRows = localRowMap_->getNodeNumElements ();
728 
729  // FIXME (mfh 14 Apr 2014) At some point, we would like to
730  // parallelize this using Kokkos. This would require a
731  // Kokkos-friendly version of getLocalRowCopy, perhaps with
732  // thread-local storage.
733 
734  const bool constantStride = X.isConstantStride () && Y.isConstantStride ();
735  if (constantStride) {
736  // Since both X and Y have constant stride, we can work with "1-D"
737  // views of their data.
738  const size_t x_stride = X.getStride();
739  const size_t y_stride = Y.getStride();
740  ArrayRCP<scalar_type> y_rcp = Y.get1dViewNonConst();
741  ArrayRCP<const scalar_type> x_rcp = X.get1dView();
742  ArrayView<scalar_type> y_ptr = y_rcp();
743  ArrayView<const scalar_type> x_ptr = x_rcp();
744  for (size_t i = 0; i < numRows; ++i) {
745  size_t Nnz;
746  // Use this class's getrow to make the below code simpler
747  getLocalRowCopy (i, localIndices_ (), Values_ (), Nnz);
748  if (mode == Teuchos::NO_TRANS) {
749  for (size_t j = 0; j < Nnz; ++j) {
750  const local_ordinal_type col = localIndices_[j];
751  for (size_t k = 0; k < NumVectors; ++k) {
752  y_ptr[i + y_stride*k] +=
753  alpha * Values_[j] * x_ptr[col + x_stride*k];
754  }
755  }
756  }
757  else if (mode == Teuchos::TRANS) {
758  for (size_t j = 0; j < Nnz; ++j) {
759  const local_ordinal_type col = localIndices_[j];
760  for (size_t k = 0; k < NumVectors; ++k) {
761  y_ptr[col + y_stride*k] +=
762  alpha * Values_[j] * x_ptr[i + x_stride*k];
763  }
764  }
765  }
766  else { //mode==Teuchos::CONJ_TRANS
767  for (size_t j = 0; j < Nnz; ++j) {
768  const local_ordinal_type col = localIndices_[j];
769  for (size_t k = 0; k < NumVectors; ++k) {
770  y_ptr[col + y_stride*k] +=
771  alpha * STS::conjugate (Values_[j]) * x_ptr[i + x_stride*k];
772  }
773  }
774  }
775  }
776  }
777  else {
778  // At least one of X or Y does not have constant stride.
779  // This means we must work with "2-D" views of their data.
780  ArrayRCP<ArrayRCP<const scalar_type> > x_ptr = X.get2dView();
781  ArrayRCP<ArrayRCP<scalar_type> > y_ptr = Y.get2dViewNonConst();
782 
783  for (size_t i = 0; i < numRows; ++i) {
784  size_t Nnz;
785  // Use this class's getrow to make the below code simpler
786  getLocalRowCopy (i, localIndices_ (), Values_ (), Nnz);
787  if (mode == Teuchos::NO_TRANS) {
788  for (size_t k = 0; k < NumVectors; ++k) {
789  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
790  ArrayView<scalar_type> y_local = (y_ptr())[k]();
791  for (size_t j = 0; j < Nnz; ++j) {
792  y_local[i] += alpha * Values_[j] * x_local[localIndices_[j]];
793  }
794  }
795  }
796  else if (mode == Teuchos::TRANS) {
797  for (size_t k = 0; k < NumVectors; ++k) {
798  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
799  ArrayView<scalar_type> y_local = (y_ptr())[k]();
800  for (size_t j = 0; j < Nnz; ++j) {
801  y_local[localIndices_[j]] += alpha * Values_[j] * x_local[i];
802  }
803  }
804  }
805  else { //mode==Teuchos::CONJ_TRANS
806  for (size_t k = 0; k < NumVectors; ++k) {
807  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
808  ArrayView<scalar_type> y_local = (y_ptr())[k]();
809  for (size_t j = 0; j < Nnz; ++j) {
810  y_local[localIndices_[j]] +=
811  alpha * STS::conjugate (Values_[j]) * x_local[i];
812  }
813  }
814  }
815  }
816  }
817 }
818 
819 template<class MatrixType>
821 {
822  return true;
823 }
824 
825 
826 template<class MatrixType>
828 {
829  return false;
830 }
831 
832 
833 template<class MatrixType>
834 typename
835 LocalFilter<MatrixType>::mag_type
837 {
838 #ifdef TPETRA_HAVE_KOKKOS_REFACTOR
839  typedef Kokkos::Details::ArithTraits<scalar_type> STS;
840  typedef Kokkos::Details::ArithTraits<mag_type> STM;
841 #else
844 #endif
845  typedef typename Teuchos::Array<scalar_type>::size_type size_type;
846 
847  const size_type maxNumRowEnt = getNodeMaxNumRowEntries ();
848  Teuchos::Array<local_ordinal_type> ind (maxNumRowEnt);
849  Teuchos::Array<scalar_type> val (maxNumRowEnt);
850  const size_t numRows = static_cast<size_t> (localRowMap_->getNodeNumElements ());
851 
852  // FIXME (mfh 03 Apr 2013) Scale during sum to avoid overflow.
853  mag_type sumSquared = STM::zero ();
854  for (size_t i = 0; i < numRows; ++i) {
855  size_t numEntries = 0;
856  this->getLocalRowCopy (i, ind (), val (), numEntries);
857  for (size_type k = 0; k < static_cast<size_type> (numEntries); ++k) {
858  const mag_type v_k_abs = STS::magnitude (val[k]);
859  sumSquared += v_k_abs * v_k_abs;
860  }
861  }
862  return STM::squareroot (sumSquared); // Different for each process; that's OK.
863 }
864 
865 template<class MatrixType>
866 std::string
868 {
870  std::ostringstream os;
871 
872  os << "Ifpack2::LocalFilter: {";
873  os << "MatrixType: " << TypeNameTraits<MatrixType>::name ();
874  if (this->getObjectLabel () != "") {
875  os << ", Label: \"" << this->getObjectLabel () << "\"";
876  }
877  os << ", Number of rows: " << getGlobalNumRows ()
878  << ", Number of columns: " << getGlobalNumCols ()
879  << "}";
880  return os.str ();
881 }
882 
883 
884 template<class MatrixType>
885 void
888  const Teuchos::EVerbosityLevel verbLevel) const
889 {
890  using Teuchos::OSTab;
892  using std::endl;
893 
894  const Teuchos::EVerbosityLevel vl =
895  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
896 
897  if (vl > Teuchos::VERB_NONE) {
898  // describe() starts with a tab, by convention.
899  OSTab tab0 (out);
900 
901  out << "Ifpack2::LocalFilter:" << endl;
902  OSTab tab1 (out);
903  out << "MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
904  if (this->getObjectLabel () != "") {
905  out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
906  }
907  out << "Number of rows: " << getGlobalNumRows () << endl
908  << "Number of columns: " << getGlobalNumCols () << endl
909  << "Number of nonzeros: " << NumNonzeros_ << endl;
910 
911  if (vl > Teuchos::VERB_LOW) {
912  out << "Row Map:" << endl;
913  localRowMap_->describe (out, vl);
914  out << "Domain Map:" << endl;
915  localDomainMap_->describe (out, vl);
916  out << "Range Map:" << endl;
917  localRangeMap_->describe (out, vl);
918  }
919  }
920 }
921 
922 template<class MatrixType>
923 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
924  typename MatrixType::local_ordinal_type,
925  typename MatrixType::global_ordinal_type,
926  typename MatrixType::node_type> >
928 {
929  return A_;
930 }
931 
932 
933 } // namespace Ifpack2
934 
935 #define IFPACK2_LOCALFILTER_INSTANT(S,LO,GO,N) \
936  template class Ifpack2::LocalFilter< Tpetra::RowMatrix<S, LO, GO, N> >;
937 
938 #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:598
virtual bool isFillComplete() const
Returns true if fillComplete() has been called.
Definition: Ifpack2_LocalFilter_def.hpp:427
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:635
virtual mag_type getFrobeniusNorm() const
The Frobenius norm of the (locally filtered) matrix.
Definition: Ifpack2_LocalFilter_def.hpp:836
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:927
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 isLocallyIndexed() const
Whether the underlying sparse matrix is locally (opposite of globally) indexed.
Definition: Ifpack2_LocalFilter_def.hpp:413
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:887
#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:586
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:204
virtual bool supportsRowViews() const
Returns true if RowViews are supported.
Definition: Ifpack2_LocalFilter_def.hpp:827
Ordinal size_type
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_LocalFilter_decl.hpp:184
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:436
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:480
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:625
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:645
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:406
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:392
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:867
virtual size_t getNodeMaxNumRowEntries() const
The maximum number of entries across all rows/columns on this process.
Definition: Ifpack2_LocalFilter_def.hpp:399
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:181
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:420
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:610
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:820
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