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
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  using Teuchos::ArrayView;
83  typedef global_ordinal_type GO; // a handy abbreviation
84  typedef typename ArrayView<const GO>::size_type size_type;
85 
86  bool fitted = true;
87  if (&map1 == &map2) {
88  fitted = true;
89  }
90  else if (map1.isContiguous () && map2.isContiguous () &&
91  map1.getMinGlobalIndex () == map2.getMinGlobalIndex () &&
92  map1.getMaxGlobalIndex () <= map2.getMaxGlobalIndex ()) {
93  // Special case where both Maps are contiguous.
94  fitted = true;
95  }
96  else {
97  ArrayView<const GO> inds_map2 = map2.getNodeElementList ();
98  const size_type numInds_map1 = static_cast<size_type> (map1.getNodeNumElements ());
99 
100  if (map1.isContiguous ()) {
101  // Avoid calling getNodeElementList() on the always one-to-one
102  // Map, if it is contiguous (a common case). When called on a
103  // contiguous Map, getNodeElementList() causes allocation of an
104  // array that sticks around, even though the array isn't needed.
105  // (The Map is contiguous, so you can compute the entries; you
106  // don't need to store them.)
107  if (numInds_map1 > inds_map2.size ()) {
108  // There are fewer indices in map1 on this process than in
109  // map2. This case might be impossible.
110  fitted = false;
111  }
112  else {
113  // Do all the map1 indices match the initial map2 indices?
114  const GO minInd_map1 = map1.getMinGlobalIndex ();
115  for (size_type k = 0; k < numInds_map1; ++k) {
116  const GO inds_map1_k = static_cast<GO> (k) + minInd_map1;
117  if (inds_map1_k != inds_map2[k]) {
118  fitted = false;
119  break;
120  }
121  }
122  }
123  }
124  else { // map1 is not contiguous.
125  // Get index lists from both Maps, and compare their indices.
126  ArrayView<const GO> inds_map1 = map1.getNodeElementList ();
127  if (numInds_map1 > inds_map2.size ()) {
128  // There are fewer indices in map1 on this process than in
129  // map2. This case might be impossible.
130  fitted = false;
131  }
132  else {
133  // Do all the map1 indices match those in map2?
134  for (size_type k = 0; k < numInds_map1; ++k) {
135  if (inds_map1[k] != inds_map2[k]) {
136  fitted = false;
137  break;
138  }
139  }
140  }
141  }
142  }
143  return fitted;
144 }
145 
146 
147 template<class MatrixType>
149 LocalFilter (const Teuchos::RCP<const row_matrix_type>& A) :
150  A_ (A),
151  NumNonzeros_ (0),
152  MaxNumEntries_ (0),
153  MaxNumEntriesA_ (0)
154 {
155  using Teuchos::RCP;
156  using Teuchos::rcp;
157 
158 #ifdef HAVE_IFPACK2_DEBUG
159  TEUCHOS_TEST_FOR_EXCEPTION(
160  ! mapPairsAreFitted (*A), std::invalid_argument, "Ifpack2::LocalFilter: "
161  "A's Map pairs are not fitted to each other on Process "
162  << A_->getRowMap ()->getComm ()->getRank () << " of the input matrix's "
163  "communicator. "
164  "This means that LocalFilter does not currently know how to work with A. "
165  "This will change soon. Please see discussion of Bug 5992.");
166 #endif // HAVE_IFPACK2_DEBUG
167 
168  // Build the local communicator (containing this process only).
169  RCP<const Teuchos::Comm<int> > localComm;
170 #ifdef HAVE_MPI
171  localComm = rcp (new Teuchos::MpiComm<int> (MPI_COMM_SELF));
172 #else
173  localComm = rcp (new Teuchos::SerialComm<int> ());
174 #endif // HAVE_MPI
175 
176 
177  // // FIXME (mfh 21 Nov 2013) Currently, the implementation implicitly
178  // // assumes that the range Map is fitted to the row Map. Otherwise,
179  // // it probably won't work at all.
180  // TEUCHOS_TEST_FOR_EXCEPTION(
181  // mapPairIsFitted (* (A_->getRangeMap ()), * (A_->getRowMap ())),
182  // std::logic_error, "Ifpack2::LocalFilter: Range Map of the input matrix "
183  // "is not fitted to its row Map. LocalFilter does not currently work in "
184  // "this case. See Bug 5992.");
185 
186  // Build the local row, domain, and range Maps. They both use the
187  // local communicator built above. The global indices of each are
188  // different than those of the corresponding original Map; they
189  // actually are the same as the local indices of the original Map.
190  //
191  // It's even OK if signedness of local_ordinal_type and
192  // global_ordinal_type are different. (That would be a BAD IDEA but
193  // some users seem to enjoy making trouble for themselves and us.)
194  //
195  // Both the local row and local range Maps must have the same number
196  // of entries, namely, that of the local number of entries of A's
197  // range Map.
198 
199  // FIXME (mfh 21 Nov 2013) For some reason, we have to use this,
200  // even if it differs from the number of entries in the range Map.
201  // Otherwise, AdditiveSchwarz Test1 fails, down in the local solve,
202  // where the matrix has 8 columns but the local part of the vector
203  // only has five rows.
204  const size_t numRows = A_->getNodeNumRows ();
205 
206  // using std::cerr;
207  // using std::endl;
208  // cerr << "A_ has " << A_->getNodeNumRows () << " rows." << endl
209  // << "Range Map has " << A_->getRangeMap ()->getNodeNumElements () << " entries." << endl
210  // << "Row Map has " << A_->getRowMap ()->getNodeNumElements () << " entries." << endl;
211 
212  const global_ordinal_type indexBase = static_cast<global_ordinal_type> (0);
213 
214  localRowMap_ =
215  rcp (new map_type (numRows, indexBase, localComm,
216  Tpetra::GloballyDistributed, A_->getNode ()));
217  // If the original matrix's range Map is not fitted to its row Map,
218  // we'll have to do an Export when applying the matrix.
219  localRangeMap_ = localRowMap_;
220 
221  // If the original matrix's domain Map is not fitted to its column
222  // Map, we'll have to do an Import when applying the matrix.
223  const size_t numCols = A_->getDomainMap ()->getNodeNumElements ();
224  if (A_->getRangeMap ().getRawPtr () == A_->getDomainMap ().getRawPtr ()) {
225  // The input matrix's domain and range Maps are identical, so the
226  // locally filtered matrix's domain and range Maps can be also.
227  localDomainMap_ = localRangeMap_;
228  }
229  else {
230  localDomainMap_ =
231  rcp (new map_type (numCols, indexBase, localComm,
232  Tpetra::GloballyDistributed, A_->getNode ()));
233  }
234 
235  // NodeNumEntries_ will contain the actual number of nonzeros for
236  // each localized row (that is, without external nodes, and always
237  // with the diagonal entry)
238  NumEntries_.resize (numRows);
239 
240  // tentative value for MaxNumEntries. This is the number of
241  // nonzeros in the local matrix
242  MaxNumEntries_ = A_->getNodeMaxNumRowEntries ();
243  MaxNumEntriesA_ = A_->getNodeMaxNumRowEntries ();
244 
245  // Allocate temporary arrays for getLocalRowCopy().
246  localIndices_.resize (MaxNumEntries_);
247  Values_.resize (MaxNumEntries_);
248 
249  // now compute:
250  // - the number of nonzero per row
251  // - the total number of nonzeros
252  // - the diagonal entries
253 
254  // compute nonzeros (total and per-row), and store the
255  // diagonal entries (already modified)
256  size_t ActualMaxNumEntries = 0;
257 
258  for (size_t i = 0; i < numRows; ++i) {
259  NumEntries_[i] = 0;
260  size_t Nnz, NewNnz = 0;
261  A_->getLocalRowCopy (i, localIndices_, Values_, Nnz);
262  for (size_t j = 0; j < Nnz; ++j) {
263  // FIXME (mfh 03 Apr 2013) This assumes the following:
264  //
265  // 1. Row Map, range Map, and domain Map are all the same.
266  //
267  // 2. The column Map's list of GIDs on this process is the
268  // domain Map's list of GIDs, followed by remote GIDs. Thus,
269  // for any GID in the domain Map on this process, its LID in
270  // the domain Map (and therefore in the row Map, by (1)) is
271  // the same as its LID in the column Map. (Hence the
272  // less-than test, which if true, means that localIndices_[j]
273  // belongs to the row Map.)
274  if (static_cast<size_t> (localIndices_[j]) < numRows) {
275  ++NewNnz;
276  }
277  }
278 
279  if (NewNnz > ActualMaxNumEntries) {
280  ActualMaxNumEntries = NewNnz;
281  }
282 
283  NumNonzeros_ += NewNnz;
284  NumEntries_[i] = NewNnz;
285  }
286 
287  MaxNumEntries_ = ActualMaxNumEntries;
288 }
289 
290 
291 template<class MatrixType>
293 {}
294 
295 
296 template<class MatrixType>
297 Teuchos::RCP<const Teuchos::Comm<int> >
299 {
300  return localRowMap_->getComm ();
301 }
302 
303 
304 template<class MatrixType>
305 Teuchos::RCP<typename MatrixType::node_type>
307 {
308  return A_->getNode ();
309 }
310 
311 
312 template<class MatrixType>
313 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
314  typename MatrixType::global_ordinal_type,
315  typename MatrixType::node_type> >
317 {
318  return localRowMap_;
319 }
320 
321 
322 template<class MatrixType>
323 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
324  typename MatrixType::global_ordinal_type,
325  typename MatrixType::node_type> >
327 {
328  return localRowMap_; // FIXME (mfh 20 Nov 2013)
329 }
330 
331 
332 template<class MatrixType>
333 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
334  typename MatrixType::global_ordinal_type,
335  typename MatrixType::node_type> >
337 {
338  return localDomainMap_;
339 }
340 
341 
342 template<class MatrixType>
343 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
344  typename MatrixType::global_ordinal_type,
345  typename MatrixType::node_type> >
347 {
348  return localRangeMap_;
349 }
350 
351 
352 template<class MatrixType>
353 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type,
354  typename MatrixType::global_ordinal_type,
355  typename MatrixType::node_type> >
357 {
358  // FIXME (mfh 20 Nov 2013) This is not what the documentation says
359  // this method should do! It should return the graph of the locally
360  // filtered matrix, not the original matrix's graph.
361  return A_->getGraph ();
362 }
363 
364 
365 template<class MatrixType>
367 {
368  return static_cast<global_size_t> (localRangeMap_->getNodeNumElements ());
369 }
370 
371 
372 template<class MatrixType>
374 {
375  return static_cast<global_size_t> (localDomainMap_->getNodeNumElements ());
376 }
377 
378 
379 template<class MatrixType>
381 {
382  return static_cast<size_t> (localRangeMap_->getNodeNumElements ());
383 }
384 
385 
386 template<class MatrixType>
388 {
389  return static_cast<size_t> (localDomainMap_->getNodeNumElements ());
390 }
391 
392 
393 template<class MatrixType>
394 typename MatrixType::global_ordinal_type
396 {
397  return A_->getIndexBase ();
398 }
399 
400 
401 template<class MatrixType>
403 {
404  return NumNonzeros_;
405 }
406 
407 
408 template<class MatrixType>
410 {
411  return NumNonzeros_;
412 }
413 
414 
415 template<class MatrixType>
416 size_t
419 {
420  const local_ordinal_type localRow = getRowMap ()->getLocalElement (globalRow);
421  if (localRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
422  // NOTE (mfh 26 Mar 2014) We return zero if globalRow is not in
423  // the row Map on this process, since "get the number of entries
424  // in the global row" refers only to what the calling process owns
425  // in that row. In this case, it owns no entries in that row,
426  // since it doesn't own the row.
427  return 0;
428  } else {
429  return NumEntries_[localRow];
430  }
431 }
432 
433 
434 template<class MatrixType>
435 size_t
438 {
439  // FIXME (mfh 07 Jul 2014) Shouldn't localRow be a local row index
440  // in the matrix's row Map, not in the LocalFilter's row Map? The
441  // latter is different; it even has different global indices!
442  // (Maybe _that_'s the bug.)
443 
444  if (getRowMap ()->isNodeLocalElement (localRow)) {
445  return NumEntries_[localRow];
446  } else {
447  // NOTE (mfh 26 Mar 2014) We return zero if localRow is not in the
448  // row Map on this process, since "get the number of entries in
449  // the local row" refers only to what the calling process owns in
450  // that row. In this case, it owns no entries in that row, since
451  // it doesn't own the row.
452  return 0;
453  }
454 }
455 
456 
457 template<class MatrixType>
459 {
460  return A_->getGlobalNumDiags ();
461 }
462 
463 
464 template<class MatrixType>
466 {
467  return A_->getNodeNumDiags ();
468 }
469 
470 
471 template<class MatrixType>
473 {
474  return MaxNumEntries_;
475 }
476 
477 
478 template<class MatrixType>
480 {
481  return MaxNumEntries_;
482 }
483 
484 
485 template<class MatrixType>
487 {
488  return true;
489 }
490 
491 
492 template<class MatrixType>
494 {
495  return A_->isLowerTriangular();
496 }
497 
498 
499 template<class MatrixType>
501 {
502  return A_->isUpperTriangular();
503 }
504 
505 
506 template<class MatrixType>
508 {
509  return A_->isLocallyIndexed ();
510 }
511 
512 
513 template<class MatrixType>
515 {
516  return A_->isGloballyIndexed();
517 }
518 
519 
520 template<class MatrixType>
522 {
523  return A_->isFillComplete ();
524 }
525 
526 
527 template<class MatrixType>
528 void
531  const Teuchos::ArrayView<global_ordinal_type>& globalIndices,
532  const Teuchos::ArrayView<scalar_type>& values,
533  size_t& numEntries) const
534 {
535  typedef local_ordinal_type LO;
536  typedef typename Teuchos::Array<LO>::size_type size_type;
537 
538  const LO localRow = this->getRowMap ()->getLocalElement (globalRow);
539  if (localRow == Teuchos::OrdinalTraits<LO>::invalid ()) {
540  // NOTE (mfh 26 Mar 2014) We return no entries if globalRow is not
541  // in the row Map on this process, since "get a copy of the
542  // entries in the global row" refers only to what the calling
543  // process owns in that row. In this case, it owns no entries in
544  // that row, since it doesn't own the row.
545  numEntries = 0;
546  }
547  else {
548  // First get a copy of the current row using local indices. Then,
549  // convert to global indices using the input matrix's column Map.
550  //
551  numEntries = this->getNumEntriesInLocalRow (localRow);
552  // FIXME (mfh 26 Mar 2014) If local_ordinal_type ==
553  // global_ordinal_type, we could just alias the input array
554  // instead of allocating a temporary array.
555  Teuchos::Array<LO> localIndices (numEntries);
556  this->getLocalRowCopy (localRow, localIndices (), values, numEntries);
557 
558  const map_type& colMap = * (this->getColMap ());
559 
560  // Don't fill the output array beyond its size.
561  const size_type numEnt =
562  std::min (static_cast<size_type> (numEntries),
563  std::min (globalIndices.size (), values.size ()));
564  for (size_type k = 0; k < numEnt; ++k) {
565  globalIndices[k] = colMap.getGlobalElement (localIndices[k]);
566  }
567  }
568 }
569 
570 
571 template<class MatrixType>
572 void
575  const Teuchos::ArrayView<local_ordinal_type> &Indices,
576  const Teuchos::ArrayView<scalar_type> &Values,
577  size_t &NumEntries) const
578 {
579  typedef local_ordinal_type LO;
580  typedef global_ordinal_type GO;
581 
582  if (! A_->getRowMap ()->isNodeLocalElement (LocalRow)) {
583  // The calling process owns zero entries in the row.
584  NumEntries = 0;
585  return;
586  }
587 
588  const size_t numEntInLclRow = NumEntries_[LocalRow];
589  if (static_cast<size_t> (Indices.size ()) < numEntInLclRow ||
590  static_cast<size_t> (Values.size ()) < numEntInLclRow) {
591  // FIXME (mfh 07 Jul 2014) Return an error code instead of
592  // throwing. We should really attempt to fill as much space as
593  // we're given, in this case.
594  TEUCHOS_TEST_FOR_EXCEPTION(
595  true, std::runtime_error,
596  "Ifpack2::LocalFilter::getLocalRowCopy: Invalid output array length. "
597  "The output arrays must each have length at least " << numEntInLclRow
598  << " for local row " << LocalRow << " on Process "
599  << localRowMap_->getComm ()->getRank () << ".");
600  }
601  else if (numEntInLclRow == static_cast<size_t> (0)) {
602  // getNumEntriesInLocalRow() returns zero if LocalRow is not owned
603  // by the calling process. In that case, the calling process owns
604  // zero entries in the row.
605  NumEntries = 0;
606  return;
607  }
608 
609  // Always extract using the temporary arrays Values_ and
610  // localIndices_. This is because I may need more space than that
611  // given by the user. The users expects only the local (in the
612  // domain Map) column indices, but I have to extract both local and
613  // remote (not in the domain Map) column indices.
614  //
615  // FIXME (mfh 07 Jul 2014) Use of temporary local storage is not
616  // conducive to thread parallelism. A better way would be to change
617  // the interface so that it only extracts values for the "local"
618  // column indices. CrsMatrix could take a set of column indices,
619  // and return their corresponding values.
620  size_t numEntInMat = 0;
621  A_->getLocalRowCopy (LocalRow, localIndices_ (), Values_ (), numEntInMat);
622 
623  // Fill the user's arrays with the "local" indices and values in
624  // that row. Note that the matrix might have a different column Map
625  // than the local filter.
626  const map_type& matrixDomMap = * (A_->getDomainMap ());
627  const map_type& matrixColMap = * (A_->getColMap ());
628 
629  const size_t capacity = static_cast<size_t> (std::min (Indices.size (),
630  Values.size ()));
631  NumEntries = 0;
632  const size_t numRows = localRowMap_->getNodeNumElements (); // superfluous
633  const bool buggy = true; // mfh 07 Jul 2014: See FIXME below.
634  for (size_t j = 0; j < numEntInMat; ++j) {
635  // The LocalFilter only includes entries in the domain Map on
636  // the calling process. We figure out whether an entry is in
637  // the domain Map by converting the (matrix column Map) index to
638  // a global index, and then asking whether that global index is
639  // in the domain Map.
640  const LO matrixLclCol = localIndices_[j];
641  const GO gblCol = matrixColMap.getGlobalElement (matrixLclCol);
642 
643  // FIXME (mfh 07 Jul 2014) This is the likely center of Bug 5992
644  // and perhaps other bugs, like Bug 6117. If 'buggy' is true,
645  // Ifpack2 tests pass; if 'buggy' is false, the tests don't pass.
646  // This suggests that Ifpack2 classes could be using LocalFilter
647  // incorrectly, perhaps by giving it an incorrect domain Map.
648  if (buggy) {
649  // only local indices
650  if ((size_t) localIndices_[j] < numRows) {
651  Indices[NumEntries] = localIndices_[j];
652  Values[NumEntries] = Values_[j];
653  NumEntries++;
654  }
655  } else {
656  if (matrixDomMap.isNodeGlobalElement (gblCol)) {
657  // Don't fill more space than the user gave us. It's an error
658  // for them not to give us enough space, but we still shouldn't
659  // overwrite memory that doesn't belong to us.
660  if (NumEntries < capacity) {
661  Indices[NumEntries] = matrixLclCol;
662  Values[NumEntries] = Values_[j];
663  }
664  NumEntries++;
665  }
666  }
667  }
668 }
669 
670 
671 template<class MatrixType>
672 void
675  Teuchos::ArrayView<const global_ordinal_type> &indices,
676  Teuchos::ArrayView<const scalar_type> &values) const
677 {
678  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
679  "Ifpack2::LocalFilter does not implement getGlobalRowView.");
680 }
681 
682 
683 template<class MatrixType>
684 void
687  Teuchos::ArrayView<const local_ordinal_type> &indices,
688  Teuchos::ArrayView<const scalar_type> &values) const
689 {
690  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
691  "Ifpack2::LocalFilter does not implement getLocalRowView.");
692 }
693 
694 
695 template<class MatrixType>
696 void
698 getLocalDiagCopy (Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& diag) const
699 {
700  using Teuchos::RCP;
701  typedef Tpetra::Vector<scalar_type, local_ordinal_type,
702  global_ordinal_type, node_type> vector_type;
703  // This is always correct, and doesn't require a collective check
704  // for sameness of Maps.
705  RCP<vector_type> diagView = diag.offsetViewNonConst (A_->getRowMap (), 0);
706  A_->getLocalDiagCopy (*diagView);
707 }
708 
709 
710 template<class MatrixType>
711 void
713 leftScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
714 {
715  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
716  "Ifpack2::LocalFilter does not implement leftScale.");
717 }
718 
719 
720 template<class MatrixType>
721 void
723 rightScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
724 {
725  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
726  "Ifpack2::LocalFilter does not implement rightScale.");
727 }
728 
729 
730 template<class MatrixType>
731 void
733 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
734  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
735  Teuchos::ETransp mode,
736  scalar_type alpha,
737  scalar_type beta) const
738 {
739  typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV;
740  TEUCHOS_TEST_FOR_EXCEPTION(
741  X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
742  "Ifpack2::LocalFilter::apply: X and Y must have the same number of columns. "
743  "X has " << X.getNumVectors () << " columns, but Y has "
744  << Y.getNumVectors () << " columns.");
745 
746 #ifdef HAVE_IFPACK2_DEBUG
747  {
748  typedef Teuchos::ScalarTraits<magnitude_type> STM;
749  Teuchos::Array<magnitude_type> norms (X.getNumVectors ());
750  X.norm1 (norms ());
751  bool good = true;
752  for (size_t j = 0; j < X.getNumVectors (); ++j) {
753  if (STM::isnaninf (norms[j])) {
754  good = false;
755  break;
756  }
757  }
758  TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error, "Ifpack2::LocalFilter::apply: The 1-norm of the input X is NaN or Inf.");
759  }
760 #endif // HAVE_IFPACK2_DEBUG
761 
762  if (&X == &Y) {
763  // FIXME (mfh 23 Apr 2014) We have to do more work to figure out
764  // if X and Y alias one another. For example, we should check
765  // whether one is a noncontiguous view of the other.
766  //
767  // X and Y alias one another, so we have to copy X.
768  MV X_copy (X, Teuchos::Copy);
769  applyNonAliased (X_copy, Y, mode, alpha, beta);
770  } else {
771  applyNonAliased (X, Y, mode, alpha, beta);
772  }
773 
774 #ifdef HAVE_IFPACK2_DEBUG
775  {
776  typedef Teuchos::ScalarTraits<magnitude_type> STM;
777  Teuchos::Array<magnitude_type> norms (Y.getNumVectors ());
778  Y.norm1 (norms ());
779  bool good = true;
780  for (size_t j = 0; j < Y.getNumVectors (); ++j) {
781  if (STM::isnaninf (norms[j])) {
782  good = false;
783  break;
784  }
785  }
786  TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error, "Ifpack2::LocalFilter::apply: The 1-norm of the output Y is NaN or Inf.");
787  }
788 #endif // HAVE_IFPACK2_DEBUG
789 }
790 
791 template<class MatrixType>
792 void
794 applyNonAliased (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
795  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
796  Teuchos::ETransp mode,
797  scalar_type alpha,
798  scalar_type beta) const
799 {
800  using Teuchos::ArrayView;
801  using Teuchos::ArrayRCP;
802  typedef Teuchos::ScalarTraits<scalar_type> STS;
803 
804  const scalar_type zero = STS::zero ();
805  const scalar_type one = STS::one ();
806 
807  if (beta == zero) {
808  Y.putScalar (zero);
809  }
810  else if (beta != one) {
811  Y.scale (beta);
812  }
813 
814  const size_t NumVectors = Y.getNumVectors ();
815  const size_t numRows = localRowMap_->getNodeNumElements ();
816 
817  // FIXME (mfh 14 Apr 2014) At some point, we would like to
818  // parallelize this using Kokkos. This would require a
819  // Kokkos-friendly version of getLocalRowCopy, perhaps with
820  // thread-local storage.
821 
822  const bool constantStride = X.isConstantStride () && Y.isConstantStride ();
823  if (constantStride) {
824  // Since both X and Y have constant stride, we can work with "1-D"
825  // views of their data.
826  const size_t x_stride = X.getStride();
827  const size_t y_stride = Y.getStride();
828  ArrayRCP<scalar_type> y_rcp = Y.get1dViewNonConst();
829  ArrayRCP<const scalar_type> x_rcp = X.get1dView();
830  ArrayView<scalar_type> y_ptr = y_rcp();
831  ArrayView<const scalar_type> x_ptr = x_rcp();
832  for (size_t i = 0; i < numRows; ++i) {
833  size_t Nnz;
834  // Use this class's getrow to make the below code simpler
835  getLocalRowCopy (i, localIndices_ (), Values_ (), Nnz);
836  if (mode == Teuchos::NO_TRANS) {
837  for (size_t j = 0; j < Nnz; ++j) {
838  const local_ordinal_type col = localIndices_[j];
839  for (size_t k = 0; k < NumVectors; ++k) {
840  y_ptr[i + y_stride*k] +=
841  alpha * Values_[j] * x_ptr[col + x_stride*k];
842  }
843  }
844  }
845  else if (mode == Teuchos::TRANS) {
846  for (size_t j = 0; j < Nnz; ++j) {
847  const local_ordinal_type col = localIndices_[j];
848  for (size_t k = 0; k < NumVectors; ++k) {
849  y_ptr[col + y_stride*k] +=
850  alpha * Values_[j] * x_ptr[i + x_stride*k];
851  }
852  }
853  }
854  else { //mode==Teuchos::CONJ_TRANS
855  for (size_t j = 0; j < Nnz; ++j) {
856  const local_ordinal_type col = localIndices_[j];
857  for (size_t k = 0; k < NumVectors; ++k) {
858  y_ptr[col + y_stride*k] +=
859  alpha * STS::conjugate (Values_[j]) * x_ptr[i + x_stride*k];
860  }
861  }
862  }
863  }
864  }
865  else {
866  // At least one of X or Y does not have constant stride.
867  // This means we must work with "2-D" views of their data.
868  ArrayRCP<ArrayRCP<const scalar_type> > x_ptr = X.get2dView();
869  ArrayRCP<ArrayRCP<scalar_type> > y_ptr = Y.get2dViewNonConst();
870 
871  for (size_t i = 0; i < numRows; ++i) {
872  size_t Nnz;
873  // Use this class's getrow to make the below code simpler
874  getLocalRowCopy (i, localIndices_ (), Values_ (), Nnz);
875  if (mode == Teuchos::NO_TRANS) {
876  for (size_t k = 0; k < NumVectors; ++k) {
877  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
878  ArrayView<scalar_type> y_local = (y_ptr())[k]();
879  for (size_t j = 0; j < Nnz; ++j) {
880  y_local[i] += alpha * Values_[j] * x_local[localIndices_[j]];
881  }
882  }
883  }
884  else if (mode == Teuchos::TRANS) {
885  for (size_t k = 0; k < NumVectors; ++k) {
886  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
887  ArrayView<scalar_type> y_local = (y_ptr())[k]();
888  for (size_t j = 0; j < Nnz; ++j) {
889  y_local[localIndices_[j]] += alpha * Values_[j] * x_local[i];
890  }
891  }
892  }
893  else { //mode==Teuchos::CONJ_TRANS
894  for (size_t k = 0; k < NumVectors; ++k) {
895  ArrayView<const scalar_type> x_local = (x_ptr())[k]();
896  ArrayView<scalar_type> y_local = (y_ptr())[k]();
897  for (size_t j = 0; j < Nnz; ++j) {
898  y_local[localIndices_[j]] +=
899  alpha * STS::conjugate (Values_[j]) * x_local[i];
900  }
901  }
902  }
903  }
904  }
905 }
906 
907 template<class MatrixType>
909 {
910  return true;
911 }
912 
913 
914 template<class MatrixType>
916 {
917  return false;
918 }
919 
920 
921 template<class MatrixType>
922 typename
923 LocalFilter<MatrixType>::mag_type
925 {
926 #ifdef TPETRA_HAVE_KOKKOS_REFACTOR
927  typedef Kokkos::Details::ArithTraits<scalar_type> STS;
928  typedef Kokkos::Details::ArithTraits<mag_type> STM;
929 #else
930  typedef Teuchos::ScalarTraits<scalar_type> STS;
931  typedef Teuchos::ScalarTraits<magnitude_type> STM;
932 #endif
933  typedef typename Teuchos::Array<scalar_type>::size_type size_type;
934 
935  const size_type maxNumRowEnt = getNodeMaxNumRowEntries ();
936  Teuchos::Array<local_ordinal_type> ind (maxNumRowEnt);
937  Teuchos::Array<scalar_type> val (maxNumRowEnt);
938  const size_t numRows = static_cast<size_t> (localRowMap_->getNodeNumElements ());
939 
940  // FIXME (mfh 03 Apr 2013) Scale during sum to avoid overflow.
941  mag_type sumSquared = STM::zero ();
942  for (size_t i = 0; i < numRows; ++i) {
943  size_t numEntries = 0;
944  this->getLocalRowCopy (i, ind (), val (), numEntries);
945  for (size_type k = 0; k < static_cast<size_type> (numEntries); ++k) {
946  const mag_type v_k_abs = STS::magnitude (val[k]);
947  sumSquared += v_k_abs * v_k_abs;
948  }
949  }
950  return STM::squareroot (sumSquared); // Different for each process; that's OK.
951 }
952 
953 template<class MatrixType>
954 std::string
956 {
957  using Teuchos::TypeNameTraits;
958  std::ostringstream os;
959 
960  os << "Ifpack2::LocalFilter: {";
961  os << "MatrixType: " << TypeNameTraits<MatrixType>::name ();
962  if (this->getObjectLabel () != "") {
963  os << ", Label: \"" << this->getObjectLabel () << "\"";
964  }
965  os << ", Number of rows: " << getGlobalNumRows ()
966  << ", Number of columns: " << getGlobalNumCols ()
967  << "}";
968  return os.str ();
969 }
970 
971 
972 template<class MatrixType>
973 void
975 describe (Teuchos::FancyOStream &out,
976  const Teuchos::EVerbosityLevel verbLevel) const
977 {
978  using Teuchos::OSTab;
979  using Teuchos::TypeNameTraits;
980  using std::endl;
981 
982  const Teuchos::EVerbosityLevel vl =
983  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
984 
985  if (vl > Teuchos::VERB_NONE) {
986  // describe() starts with a tab, by convention.
987  OSTab tab0 (out);
988 
989  out << "Ifpack2::LocalFilter:" << endl;
990  OSTab tab1 (out);
991  out << "MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
992  if (this->getObjectLabel () != "") {
993  out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
994  }
995  out << "Number of rows: " << getGlobalNumRows () << endl
996  << "Number of columns: " << getGlobalNumCols () << endl
997  << "Number of nonzeros: " << NumNonzeros_ << endl;
998 
999  if (vl > Teuchos::VERB_LOW) {
1000  out << "Row Map:" << endl;
1001  localRowMap_->describe (out, vl);
1002  out << "Domain Map:" << endl;
1003  localDomainMap_->describe (out, vl);
1004  out << "Range Map:" << endl;
1005  localRangeMap_->describe (out, vl);
1006  }
1007  }
1008 }
1009 
1010 template<class MatrixType>
1011 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
1012  typename MatrixType::local_ordinal_type,
1013  typename MatrixType::global_ordinal_type,
1014  typename MatrixType::node_type> >
1016 {
1017  return A_;
1018 }
1019 
1020 
1021 } // namespace Ifpack2
1022 
1023 #define IFPACK2_LOCALFILTER_INSTANT(S,LO,GO,N) \
1024  template class Ifpack2::LocalFilter< Tpetra::RowMatrix<S, LO, GO, N> >;
1025 
1026 #endif
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:723
virtual size_t getNodeNumDiags() const
The number of local diagonal entries, based on global row/column index comparisons.
Definition: Ifpack2_LocalFilter_def.hpp:465
virtual Teuchos::RCP< const map_type > getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:316
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:975
virtual Teuchos::RCP< node_type > getNode() const
Returns the underlying Node object.
Definition: Ifpack2_LocalFilter_def.hpp:306
virtual size_t getNodeNumRows() const
The number of rows owned on the calling process.
Definition: Ifpack2_LocalFilter_def.hpp:380
virtual Teuchos::RCP< const map_type > getDomainMap() const
Returns the Map that describes the domain distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:336
virtual bool isFillComplete() const
Returns true if fillComplete() has been called.
Definition: Ifpack2_LocalFilter_def.hpp:521
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the communicator.
Definition: Ifpack2_LocalFilter_def.hpp:298
virtual size_t getNodeMaxNumRowEntries() const
The maximum number of entries across all rows/columns on this process.
Definition: Ifpack2_LocalFilter_def.hpp:479
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:733
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_LocalFilter_decl.hpp:190
virtual bool hasTransposeApply() const
Whether this operator supports applying the transpose or conjugate transpose.
Definition: Ifpack2_LocalFilter_def.hpp:908
virtual bool isLocallyIndexed() const
Whether the underlying sparse matrix is locally (opposite of globally) indexed.
Definition: Ifpack2_LocalFilter_def.hpp:507
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:674
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:698
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 mag_type getFrobeniusNorm() const
The Frobenius norm of the (locally filtered) matrix.
Definition: Ifpack2_LocalFilter_def.hpp:924
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:356
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_LocalFilter_decl.hpp:187
virtual bool supportsRowViews() const
Returns true if RowViews are supported.
Definition: Ifpack2_LocalFilter_def.hpp:915
virtual global_size_t getGlobalNumDiags() const
The number of global diagonal entries, based on global row/column index comparisons.
Definition: Ifpack2_LocalFilter_def.hpp:458
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:418
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:713
virtual global_ordinal_type getIndexBase() const
Returns the index base for global indices for this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:395
virtual bool isGloballyIndexed() const
Whether the underlying sparse matrix is globally (opposite of locally) indexed.
Definition: Ifpack2_LocalFilter_def.hpp:514
virtual bool isUpperTriangular() const
Whether this matrix is upper triangular.
Definition: Ifpack2_LocalFilter_def.hpp:500
virtual size_t getGlobalMaxNumRowEntries() const
The maximum number of entries across all rows/columns on all processes.
Definition: Ifpack2_LocalFilter_def.hpp:472
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:574
virtual size_t getNodeNumCols() const
The number of columns in the (locally filtered) matrix.
Definition: Ifpack2_LocalFilter_def.hpp:387
LocalFilter(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_LocalFilter_def.hpp:149
virtual ~LocalFilter()
Destructor.
Definition: Ifpack2_LocalFilter_def.hpp:292
virtual global_size_t getGlobalNumRows() const
The number of global rows in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:366
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_LocalFilter_decl.hpp:184
virtual std::string description() const
A one-line description of this object.
Definition: Ifpack2_LocalFilter_def.hpp:955
virtual size_t getNodeNumEntries() const
Returns the local number of entries in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:409
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:437
virtual Teuchos::RCP< const map_type > getRangeMap() const
Returns the Map that describes the range distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:346
virtual global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:402
virtual Teuchos::RCP< const row_matrix_type > getUnderlyingMatrix() const
Return matrix that LocalFilter was built on.
Definition: Ifpack2_LocalFilter_def.hpp:1015
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:530
Access only local rows and columns of a sparse matrix.
Definition: Ifpack2_LocalFilter_decl.hpp:160
virtual global_size_t getGlobalNumCols() const
The number of global columns in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:373
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
virtual bool isLowerTriangular() const
Whether this matrix is lower triangular.
Definition: Ifpack2_LocalFilter_def.hpp:493
virtual Teuchos::RCP< const map_type > getColMap() const
Returns the Map that describes the column distribution in this matrix.
Definition: Ifpack2_LocalFilter_def.hpp:326
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:686
virtual bool hasColMap() const
Whether this matrix has a well-defined column Map.
Definition: Ifpack2_LocalFilter_def.hpp:486
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_LocalFilter_decl.hpp:193