Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_OverlappingRowMatrix_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_OVERLAPPINGROWMATRIX_DEF_HPP
44 #define IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
45 
46 #include <sstream>
47 
48 #include <Ifpack2_Details_OverlappingRowGraph.hpp>
49 #include <Tpetra_CrsMatrix.hpp>
50 #include <Tpetra_Import.hpp>
51 #include "Tpetra_Map.hpp"
52 #include <Teuchos_CommHelpers.hpp>
53 
54 namespace Ifpack2 {
55 
56 template<class MatrixType>
58 OverlappingRowMatrix (const Teuchos::RCP<const row_matrix_type>& A,
59  const int overlapLevel) :
60  A_ (Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A, true)),
61  OverlapLevel_ (overlapLevel)
62 {
63  using Teuchos::RCP;
64  using Teuchos::rcp;
65  using Teuchos::Array;
66  using Teuchos::outArg;
67  using Teuchos::REDUCE_SUM;
68  using Teuchos::reduceAll;
69  typedef Tpetra::global_size_t GST;
70  typedef Tpetra::CrsGraph<local_ordinal_type,
71  global_ordinal_type, node_type> crs_graph_type;
72  TEUCHOS_TEST_FOR_EXCEPTION(
73  OverlapLevel_ <= 0, std::runtime_error,
74  "Ifpack2::OverlappingRowMatrix: OverlapLevel must be > 0.");
75  TEUCHOS_TEST_FOR_EXCEPTION
76  (A_.is_null (), std::runtime_error,
77  "Ifpack2::OverlappingRowMatrix: The input matrix must be a "
78  "Tpetra::CrsMatrix with the same scalar_type, local_ordinal_type, "
79  "global_ordinal_type, and device_type typedefs as MatrixType.");
80  TEUCHOS_TEST_FOR_EXCEPTION(
81  A_->getComm()->getSize() == 1, std::runtime_error,
82  "Ifpack2::OverlappingRowMatrix: Matrix must be "
83  "distributed over more than one MPI process.");
84 
85  RCP<const crs_graph_type> A_crsGraph = A_->getCrsGraph ();
86  const size_t numMyRowsA = A_->getNodeNumRows ();
87  const global_ordinal_type global_invalid =
88  Teuchos::OrdinalTraits<global_ordinal_type>::invalid ();
89 
90  // Temp arrays
91  Array<global_ordinal_type> ExtElements;
92  RCP<map_type> TmpMap;
93  RCP<crs_graph_type> TmpGraph;
94  RCP<import_type> TmpImporter;
95  RCP<const map_type> RowMap, ColMap;
96  ExtHaloStarts_.resize(OverlapLevel_+1);
97 
98  // The big import loop
99  for (int overlap = 0 ; overlap < OverlapLevel_ ; ++overlap) {
100  ExtHaloStarts_[overlap] = (size_t) ExtElements.size();
101 
102  // Get the current maps
103  if (overlap == 0) {
104  RowMap = A_->getRowMap ();
105  ColMap = A_->getColMap ();
106  }
107  else {
108  RowMap = TmpGraph->getRowMap ();
109  ColMap = TmpGraph->getColMap ();
110  }
111 
112  const size_t size = ColMap->getNodeNumElements () - RowMap->getNodeNumElements ();
113  Array<global_ordinal_type> mylist (size);
114  size_t count = 0;
115 
116  // define the set of rows that are in ColMap but not in RowMap
117  for (local_ordinal_type i = 0 ; (size_t) i < ColMap->getNodeNumElements() ; ++i) {
118  const global_ordinal_type GID = ColMap->getGlobalElement (i);
119  if (A_->getRowMap ()->getLocalElement (GID) == global_invalid) {
120  typedef typename Array<global_ordinal_type>::iterator iter_type;
121  const iter_type end = ExtElements.end ();
122  const iter_type pos = std::find (ExtElements.begin (), end, GID);
123  if (pos == end) {
124  ExtElements.push_back (GID);
125  mylist[count] = GID;
126  ++count;
127  }
128  }
129  }
130 
131  // mfh 24 Nov 2013: We don't need TmpMap, TmpGraph, or
132  // TmpImporter after this loop, so we don't have to construct them
133  // on the last round.
134  if (overlap + 1 < OverlapLevel_) {
135  // Allocate & import new matrices, maps, etc.
136  //
137  // FIXME (mfh 24 Nov 2013) Do we always want to use index base
138  // zero? It doesn't really matter, since the actual index base
139  // (in the current implementation of Map) will always be the
140  // globally least GID.
141  TmpMap = rcp (new map_type (global_invalid, mylist (0, count),
142  Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
143  A_->getComm ()));
144  TmpGraph = rcp (new crs_graph_type (TmpMap, 0));
145  TmpImporter = rcp (new import_type (A_->getRowMap (), TmpMap));
146 
147  TmpGraph->doImport (*A_crsGraph, *TmpImporter, Tpetra::INSERT);
148  TmpGraph->fillComplete (A_->getDomainMap (), TmpMap);
149  }
150  } // end overlap loop
151  ExtHaloStarts_[OverlapLevel_] = (size_t) ExtElements.size();
152 
153 
154  // build the map containing all the nodes (original
155  // matrix + extended matrix)
156  Array<global_ordinal_type> mylist (numMyRowsA + ExtElements.size ());
157  for (local_ordinal_type i = 0; (size_t)i < numMyRowsA; ++i) {
158  mylist[i] = A_->getRowMap ()->getGlobalElement (i);
159  }
160  for (local_ordinal_type i = 0; i < ExtElements.size (); ++i) {
161  mylist[i + numMyRowsA] = ExtElements[i];
162  }
163 
164  RowMap_ = rcp (new map_type (global_invalid, mylist (),
165  Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
166  A_->getComm ()));
167  Importer_ = rcp (new import_type (A_->getRowMap (), RowMap_));
168  ColMap_ = RowMap_;
169 
170  // now build the map corresponding to all the external nodes
171  // (with respect to A().RowMatrixRowMap().
172  ExtMap_ = rcp (new map_type (global_invalid, ExtElements (),
173  Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
174  A_->getComm ()));
175  ExtImporter_ = rcp (new import_type (A_->getRowMap (), ExtMap_));
176 
177  {
178  RCP<crs_matrix_type> ExtMatrix_nc =
179  rcp (new crs_matrix_type (ExtMap_, ColMap_, 0));
180  ExtMatrix_nc->doImport (*A_, *ExtImporter_, Tpetra::INSERT);
181  ExtMatrix_nc->fillComplete (A_->getDomainMap (), RowMap_);
182  ExtMatrix_ = ExtMatrix_nc; // we only need the const version after here
183  }
184 
185  // fix indices for overlapping matrix
186  const size_t numMyRowsB = ExtMatrix_->getNodeNumRows ();
187 
188  GST NumMyNonzeros_tmp = A_->getNodeNumEntries () + ExtMatrix_->getNodeNumEntries ();
189  GST NumMyRows_tmp = numMyRowsA + numMyRowsB;
190  {
191  GST inArray[2], outArray[2];
192  inArray[0] = NumMyNonzeros_tmp;
193  inArray[1] = NumMyRows_tmp;
194  outArray[0] = 0;
195  outArray[1] = 0;
196  reduceAll<int, GST> (* (A_->getComm ()), REDUCE_SUM, 2, inArray, outArray);
197  NumGlobalNonzeros_ = outArray[0];
198  NumGlobalRows_ = outArray[1];
199  }
200  // reduceAll<int, GST> (* (A_->getComm ()), REDUCE_SUM, NumMyNonzeros_tmp,
201  // outArg (NumGlobalNonzeros_));
202  // reduceAll<int, GST> (* (A_->getComm ()), REDUCE_SUM, NumMyRows_tmp,
203  // outArg (NumGlobalRows_));
204 
205  MaxNumEntries_ = A_->getNodeMaxNumRowEntries ();
206  if (MaxNumEntries_ < ExtMatrix_->getNodeMaxNumRowEntries ()) {
207  MaxNumEntries_ = ExtMatrix_->getNodeMaxNumRowEntries ();
208  }
209 
210  // Create the graph (returned by getGraph()).
211  typedef Details::OverlappingRowGraph<row_graph_type> row_graph_impl_type;
212  RCP<row_graph_impl_type> graph =
213  rcp (new row_graph_impl_type (A_->getGraph (),
214  ExtMatrix_->getGraph (),
215  RowMap_,
216  ColMap_,
217  NumGlobalRows_,
218  NumGlobalRows_, // # global cols == # global rows
219  NumGlobalNonzeros_,
220  MaxNumEntries_,
221  Importer_,
222  ExtImporter_));
223  graph_ = Teuchos::rcp_const_cast<const row_graph_type>
224  (Teuchos::rcp_implicit_cast<row_graph_type> (graph));
225  // Resize temp arrays
226  Kokkos::resize(Indices_,MaxNumEntries_);
227  Kokkos::resize(Values_,MaxNumEntries_);
228 }
229 
230 
231 template<class MatrixType>
232 Teuchos::RCP<const Teuchos::Comm<int> >
234 {
235  return A_->getComm ();
236 }
237 
238 
239 
240 
241 template<class MatrixType>
242 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
244 {
245  // FIXME (mfh 12 July 2013) Is this really the right Map to return?
246  return RowMap_;
247 }
248 
249 
250 template<class MatrixType>
251 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
253 {
254  // FIXME (mfh 12 July 2013) Is this really the right Map to return?
255  return ColMap_;
256 }
257 
258 
259 template<class MatrixType>
260 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
262 {
263  // The original matrix's domain map is irrelevant; we want the map associated
264  // with the overlap. This can then be used by LocalFilter, for example, while
265  // letting LocalFilter still filter based on domain and range maps (instead of
266  // column and row maps).
267  // FIXME Ideally, this would be the same map but restricted to a local
268  // communicator. If replaceCommWithSubset were free, that would be the way to
269  // go. That would require a new Map ctor. For now, we'll stick with ColMap_'s
270  // global communicator.
271  return ColMap_;
272 }
273 
274 
275 template<class MatrixType>
276 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
278 {
279  return RowMap_;
280 }
281 
282 
283 template<class MatrixType>
284 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
286 {
287  return graph_;
288 }
289 
290 
291 template<class MatrixType>
293 {
294  return NumGlobalRows_;
295 }
296 
297 
298 template<class MatrixType>
300 {
301  return NumGlobalRows_;
302 }
303 
304 
305 template<class MatrixType>
307 {
308  return A_->getNodeNumRows () + ExtMatrix_->getNodeNumRows ();
309 }
310 
311 
312 template<class MatrixType>
314 {
315  return this->getNodeNumRows ();
316 }
317 
318 
319 template<class MatrixType>
320 typename MatrixType::global_ordinal_type
322 {
323  return A_->getIndexBase();
324 }
325 
326 
327 template<class MatrixType>
329 {
330  return NumGlobalNonzeros_;
331 }
332 
333 
334 template<class MatrixType>
336 {
337  return A_->getNodeNumEntries () + ExtMatrix_->getNodeNumEntries ();
338 }
339 
340 
341 template<class MatrixType>
342 size_t
344 getNumEntriesInGlobalRow (global_ordinal_type globalRow) const
345 {
346  const local_ordinal_type localRow = RowMap_->getLocalElement (globalRow);
347  if (localRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
348  return Teuchos::OrdinalTraits<size_t>::invalid();
349  } else {
350  return getNumEntriesInLocalRow (localRow);
351  }
352 }
353 
354 
355 template<class MatrixType>
356 size_t
358 getNumEntriesInLocalRow (local_ordinal_type localRow) const
359 {
360  using Teuchos::as;
361  const size_t numMyRowsA = A_->getNodeNumRows ();
362  if (as<size_t> (localRow) < numMyRowsA) {
363  return A_->getNumEntriesInLocalRow (localRow);
364  } else {
365  return ExtMatrix_->getNumEntriesInLocalRow (as<local_ordinal_type> (localRow - numMyRowsA));
366  }
367 }
368 
369 
370 template<class MatrixType>
372 {
373  throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getGlobalMaxNumRowEntries() not supported.");
374 }
375 
376 
377 template<class MatrixType>
379 {
380  return MaxNumEntries_;
381 }
382 
383 
384 template<class MatrixType>
386 {
387  return true;
388 }
389 
390 
391 template<class MatrixType>
393 {
394  return true;
395 }
396 
397 
398 template<class MatrixType>
400 {
401  return false;
402 }
403 
404 
405 template<class MatrixType>
407 {
408  return true;
409 }
410 
411 
412 template<class MatrixType>
413 void
415  getGlobalRowCopy (global_ordinal_type GlobalRow,
416  nonconst_global_inds_host_view_type &Indices,
417  nonconst_values_host_view_type &Values,
418  size_t& NumEntries) const
419 {
420  const local_ordinal_type LocalRow = RowMap_->getLocalElement (GlobalRow);
421  if (LocalRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
422  NumEntries = Teuchos::OrdinalTraits<size_t>::invalid ();
423  } else {
424  if (Teuchos::as<size_t> (LocalRow) < A_->getNodeNumRows ()) {
425  A_->getGlobalRowCopy (GlobalRow, Indices, Values, NumEntries);
426  } else {
427  ExtMatrix_->getGlobalRowCopy (GlobalRow, Indices, Values, NumEntries);
428  }
429  }
430 }
431 
432 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
433 template<class MatrixType>
435 getGlobalRowCopy (global_ordinal_type GlobalRow,
436  const Teuchos::ArrayView<global_ordinal_type> &Indices,
437  const Teuchos::ArrayView<scalar_type> &Values,
438  size_t &NumEntries) const {
439  using IST = typename row_matrix_type::impl_scalar_type;
440  nonconst_global_inds_host_view_type ind_in(Indices.data(),Indices.size());
441  nonconst_values_host_view_type val_in(reinterpret_cast<IST*>(Values.data()),Values.size());
442  getGlobalRowCopy(GlobalRow,ind_in,val_in,NumEntries);
443 }
444 #endif
445 
446 template<class MatrixType>
447 void
449  getLocalRowCopy (local_ordinal_type LocalRow,
450  nonconst_local_inds_host_view_type &Indices,
451  nonconst_values_host_view_type &Values,
452  size_t& NumEntries) const
453 {
454  using Teuchos::as;
455  const size_t numMyRowsA = A_->getNodeNumRows ();
456  if (as<size_t> (LocalRow) < numMyRowsA) {
457  A_->getLocalRowCopy (LocalRow, Indices, Values, NumEntries);
458  } else {
459  ExtMatrix_->getLocalRowCopy (LocalRow - as<local_ordinal_type> (numMyRowsA),
460  Indices, Values, NumEntries);
461  }
462 }
463 
464 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
465 template<class MatrixType>
466 void
468 getLocalRowCopy (local_ordinal_type LocalRow,
469  const Teuchos::ArrayView<local_ordinal_type> &Indices,
470  const Teuchos::ArrayView<scalar_type> &Values,
471  size_t &NumEntries) const
472 {
473  using IST = typename row_matrix_type::impl_scalar_type;
474  nonconst_local_inds_host_view_type ind_in(Indices.data(),Indices.size());
475  nonconst_values_host_view_type val_in(reinterpret_cast<IST*>(Values.data()),Values.size());
476  getLocalRowCopy(LocalRow,ind_in,val_in,NumEntries);
477 }
478 #endif
479 
480 template<class MatrixType>
481 void
483 getGlobalRowView (global_ordinal_type GlobalRow,
484  global_inds_host_view_type &indices,
485  values_host_view_type &values) const {
486  const local_ordinal_type LocalRow = RowMap_->getLocalElement (GlobalRow);
487  if (LocalRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid()) {
488  indices = global_inds_host_view_type();
489  values = values_host_view_type();
490  } else {
491  if (Teuchos::as<size_t> (LocalRow) < A_->getNodeNumRows ()) {
492  A_->getGlobalRowView (GlobalRow, indices, values);
493  } else {
494  ExtMatrix_->getGlobalRowView (GlobalRow, indices, values);
495  }
496  }
497 }
498 
499 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
500 template<class MatrixType>
501 void
503 getGlobalRowView (global_ordinal_type GlobalRow,
504  Teuchos::ArrayView<const global_ordinal_type>& indices,
505  Teuchos::ArrayView<const scalar_type>& values) const
506 {
507  const local_ordinal_type LocalRow = RowMap_->getLocalElement (GlobalRow);
508  if (LocalRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid()) {
509  indices = Teuchos::null;
510  values = Teuchos::null;
511  } else {
512  if (Teuchos::as<size_t> (LocalRow) < A_->getNodeNumRows ()) {
513  A_->getGlobalRowView (GlobalRow, indices, values);
514  } else {
515  ExtMatrix_->getGlobalRowView (GlobalRow, indices, values);
516  }
517  }
518 }
519 #endif
520 
521 template<class MatrixType>
522 void
524  getLocalRowView (local_ordinal_type LocalRow,
525  local_inds_host_view_type & indices,
526  values_host_view_type & values) const {
527  using Teuchos::as;
528  const size_t numMyRowsA = A_->getNodeNumRows ();
529  if (as<size_t> (LocalRow) < numMyRowsA) {
530  A_->getLocalRowView (LocalRow, indices, values);
531  } else {
532  ExtMatrix_->getLocalRowView (LocalRow - as<local_ordinal_type> (numMyRowsA),
533  indices, values);
534  }
535 
536 }
537 
538 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
539 template<class MatrixType>
540 void
542 getLocalRowView (local_ordinal_type LocalRow,
543  Teuchos::ArrayView<const local_ordinal_type>& indices,
544  Teuchos::ArrayView<const scalar_type>& values) const
545 {
546  using Teuchos::as;
547  const size_t numMyRowsA = A_->getNodeNumRows ();
548  if (as<size_t> (LocalRow) < numMyRowsA) {
549  A_->getLocalRowView (LocalRow, indices, values);
550  } else {
551  ExtMatrix_->getLocalRowView (LocalRow - as<local_ordinal_type> (numMyRowsA),
552  indices, values);
553  }
554 }
555 #endif
556 
557 
558 template<class MatrixType>
559 void
561 getLocalDiagCopy (Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& diag) const
562 {
563  using Teuchos::Array;
564 
565  //extract diagonal of original matrix
566  vector_type baseDiag(A_->getRowMap()); // diagonal of original matrix A_
567  A_->getLocalDiagCopy(baseDiag);
568  Array<scalar_type> baseDiagVals(baseDiag.getLocalLength());
569  baseDiag.get1dCopy(baseDiagVals());
570  //extra diagonal of ghost matrix
571  vector_type extDiag(ExtMatrix_->getRowMap());
572  ExtMatrix_->getLocalDiagCopy(extDiag);
573  Array<scalar_type> extDiagVals(extDiag.getLocalLength());
574  extDiag.get1dCopy(extDiagVals());
575 
576  Teuchos::ArrayRCP<scalar_type> allDiagVals = diag.getDataNonConst();
577  if (allDiagVals.size() != baseDiagVals.size() + extDiagVals.size()) {
578  std::ostringstream errStr;
579  errStr << "Ifpack2::OverlappingRowMatrix::getLocalDiagCopy : Mismatch in diagonal lengths, "
580  << allDiagVals.size() << " != " << baseDiagVals.size() << "+" << extDiagVals.size();
581  throw std::runtime_error(errStr.str());
582  }
583  for (Teuchos::Ordinal i=0; i<baseDiagVals.size(); ++i)
584  allDiagVals[i] = baseDiagVals[i];
585  Teuchos_Ordinal offset=baseDiagVals.size();
586  for (Teuchos::Ordinal i=0; i<extDiagVals.size(); ++i)
587  allDiagVals[i+offset] = extDiagVals[i];
588 }
589 
590 
591 template<class MatrixType>
592 void
594 leftScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& /* x */)
595 {
596  throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
597 }
598 
599 
600 template<class MatrixType>
601 void
603 rightScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& /* x */)
604 {
605  throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
606 }
607 
608 
609 template<class MatrixType>
610 typename OverlappingRowMatrix<MatrixType>::mag_type
612 {
613  throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support getFrobeniusNorm.");
614 }
615 
616 
617 template<class MatrixType>
618 void
620 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
621  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
622  Teuchos::ETransp mode,
623  scalar_type alpha,
624  scalar_type beta) const
625 {
626  using MV = Tpetra::MultiVector<scalar_type, local_ordinal_type,
627  global_ordinal_type, node_type>;
628  TEUCHOS_TEST_FOR_EXCEPTION
629  (X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
630  "Ifpack2::OverlappingRowMatrix::apply: X.getNumVectors() = "
631  << X.getNumVectors() << " != Y.getNumVectors() = " << Y.getNumVectors()
632  << ".");
633  // If X aliases Y, we'll need to copy X.
634  bool aliases = X.aliases(Y);
635  if (aliases) {
636  MV X_copy (X, Teuchos::Copy);
637  this->apply (X_copy, Y, mode, alpha, beta);
638  return;
639  }
640 
641  const auto& rowMap0 = * (A_->getRowMap ());
642  const auto& colMap0 = * (A_->getColMap ());
643  MV X_0 (X, mode == Teuchos::NO_TRANS ? colMap0 : rowMap0, 0);
644  MV Y_0 (Y, mode == Teuchos::NO_TRANS ? rowMap0 : colMap0, 0);
645  A_->localApply (X_0, Y_0, mode, alpha, beta);
646 
647  const auto& rowMap1 = * (ExtMatrix_->getRowMap ());
648  const auto& colMap1 = * (ExtMatrix_->getColMap ());
649  MV X_1 (X, mode == Teuchos::NO_TRANS ? colMap1 : rowMap1, 0);
650  MV Y_1 (Y, mode == Teuchos::NO_TRANS ? rowMap1 : colMap1, A_->getNodeNumRows ());
651  ExtMatrix_->localApply (X_1, Y_1, mode, alpha, beta);
652 }
653 
654 
655 template<class MatrixType>
656 void
658 importMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
659  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
660  Tpetra::CombineMode CM)
661 {
662  OvX.doImport (X, *Importer_, CM);
663 }
664 
665 
666 template<class MatrixType>
667 void
668 OverlappingRowMatrix<MatrixType>::
669 exportMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
670  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
671  Tpetra::CombineMode CM)
672 {
673  X.doExport (OvX, *Importer_, CM);
674 }
675 
676 
677 template<class MatrixType>
679 {
680  return true;
681 }
682 
683 
684 template<class MatrixType>
686 {
687  return false;
688 }
689 
690 template<class MatrixType>
692 {
693  std::ostringstream oss;
694  if (isFillComplete()) {
695  oss << "{ isFillComplete: true"
696  << ", global rows: " << getGlobalNumRows()
697  << ", global columns: " << getGlobalNumCols()
698  << ", global entries: " << getGlobalNumEntries()
699  << " }";
700  }
701  else {
702  oss << "{ isFillComplete: false"
703  << ", global rows: " << getGlobalNumRows()
704  << " }";
705  }
706  return oss.str();
707 }
708 
709 template<class MatrixType>
710 void OverlappingRowMatrix<MatrixType>::describe(Teuchos::FancyOStream &out,
711  const Teuchos::EVerbosityLevel verbLevel) const
712 {
713  using std::endl;
714  using std::setw;
715  using Teuchos::as;
716  using Teuchos::VERB_DEFAULT;
717  using Teuchos::VERB_NONE;
718  using Teuchos::VERB_LOW;
719  using Teuchos::VERB_MEDIUM;
720  using Teuchos::VERB_HIGH;
721  using Teuchos::VERB_EXTREME;
722  using Teuchos::RCP;
723  using Teuchos::null;
724  using Teuchos::ArrayView;
725 
726  Teuchos::EVerbosityLevel vl = verbLevel;
727  if (vl == VERB_DEFAULT) {
728  vl = VERB_LOW;
729  }
730  RCP<const Teuchos::Comm<int> > comm = this->getComm();
731  const int myRank = comm->getRank();
732  const int numProcs = comm->getSize();
733  size_t width = 1;
734  for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
735  ++width;
736  }
737  width = std::max<size_t> (width, as<size_t> (11)) + 2;
738  Teuchos::OSTab tab(out);
739  // none: print nothing
740  // low: print O(1) info from node 0
741  // medium: print O(P) info, num entries per process
742  // high: print O(N) info, num entries per row
743  // extreme: print O(NNZ) info: print indices and values
744  //
745  // for medium and higher, print constituent objects at specified verbLevel
746  if (vl != VERB_NONE) {
747  if (myRank == 0) {
748  out << this->description() << std::endl;
749  }
750  // O(1) globals, minus what was already printed by description()
751  //if (isFillComplete() && myRank == 0) {
752  // out << "Global max number of entries in a row: " << getGlobalMaxNumRowEntries() << std::endl;
753  //}
754  // constituent objects
755  if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
756  if (myRank == 0) {
757  out << endl << "Row map:" << endl;
758  }
759  getRowMap()->describe(out,vl);
760  //
761  if (getColMap() != null) {
762  if (getColMap() == getRowMap()) {
763  if (myRank == 0) {
764  out << endl << "Column map is row map.";
765  }
766  }
767  else {
768  if (myRank == 0) {
769  out << endl << "Column map:" << endl;
770  }
771  getColMap()->describe(out,vl);
772  }
773  }
774  if (getDomainMap() != null) {
775  if (getDomainMap() == getRowMap()) {
776  if (myRank == 0) {
777  out << endl << "Domain map is row map.";
778  }
779  }
780  else if (getDomainMap() == getColMap()) {
781  if (myRank == 0) {
782  out << endl << "Domain map is column map.";
783  }
784  }
785  else {
786  if (myRank == 0) {
787  out << endl << "Domain map:" << endl;
788  }
789  getDomainMap()->describe(out,vl);
790  }
791  }
792  if (getRangeMap() != null) {
793  if (getRangeMap() == getDomainMap()) {
794  if (myRank == 0) {
795  out << endl << "Range map is domain map." << endl;
796  }
797  }
798  else if (getRangeMap() == getRowMap()) {
799  if (myRank == 0) {
800  out << endl << "Range map is row map." << endl;
801  }
802  }
803  else {
804  if (myRank == 0) {
805  out << endl << "Range map: " << endl;
806  }
807  getRangeMap()->describe(out,vl);
808  }
809  }
810  if (myRank == 0) {
811  out << endl;
812  }
813  }
814  // O(P) data
815  if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
816  for (int curRank = 0; curRank < numProcs; ++curRank) {
817  if (myRank == curRank) {
818  out << "Process rank: " << curRank << std::endl;
819  out << " Number of entries: " << getNodeNumEntries() << std::endl;
820  out << " Max number of entries per row: " << getNodeMaxNumRowEntries() << std::endl;
821  }
822  comm->barrier();
823  comm->barrier();
824  comm->barrier();
825  }
826  }
827  // O(N) and O(NNZ) data
828  if (vl == VERB_HIGH || vl == VERB_EXTREME) {
829  for (int curRank = 0; curRank < numProcs; ++curRank) {
830  if (myRank == curRank) {
831  out << std::setw(width) << "Proc Rank"
832  << std::setw(width) << "Global Row"
833  << std::setw(width) << "Num Entries";
834  if (vl == VERB_EXTREME) {
835  out << std::setw(width) << "(Index,Value)";
836  }
837  out << endl;
838  for (size_t r = 0; r < getNodeNumRows (); ++r) {
839  const size_t nE = getNumEntriesInLocalRow(r);
840  typename MatrixType::global_ordinal_type gid = getRowMap()->getGlobalElement(r);
841  out << std::setw(width) << myRank
842  << std::setw(width) << gid
843  << std::setw(width) << nE;
844  if (vl == VERB_EXTREME) {
845  if (isGloballyIndexed()) {
846  global_inds_host_view_type rowinds;
847  values_host_view_type rowvals;
848  getGlobalRowView (gid, rowinds, rowvals);
849  for (size_t j = 0; j < nE; ++j) {
850  out << " (" << rowinds[j]
851  << ", " << rowvals[j]
852  << ") ";
853  }
854  }
855  else if (isLocallyIndexed()) {
856  local_inds_host_view_type rowinds;
857  values_host_view_type rowvals;
858  getLocalRowView (r, rowinds, rowvals);
859  for (size_t j=0; j < nE; ++j) {
860  out << " (" << getColMap()->getGlobalElement(rowinds[j])
861  << ", " << rowvals[j]
862  << ") ";
863  }
864  } // globally or locally indexed
865  } // vl == VERB_EXTREME
866  out << endl;
867  } // for each row r on this process
868 
869  } // if (myRank == curRank)
870  comm->barrier();
871  comm->barrier();
872  comm->barrier();
873  }
874 
875  out.setOutputToRootOnly(0);
876  out << "===========\nlocal matrix\n=================" << std::endl;
877  out.setOutputToRootOnly(-1);
878  A_->describe(out,Teuchos::VERB_EXTREME);
879  out.setOutputToRootOnly(0);
880  out << "===========\nend of local matrix\n=================" << std::endl;
881  comm->barrier();
882  out.setOutputToRootOnly(0);
883  out << "=================\nghost matrix\n=================" << std::endl;
884  out.setOutputToRootOnly(-1);
885  ExtMatrix_->describe(out,Teuchos::VERB_EXTREME);
886  out.setOutputToRootOnly(0);
887  out << "===========\nend of ghost matrix\n=================" << std::endl;
888  }
889  }
890 }
891 
892 template<class MatrixType>
893 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
894 OverlappingRowMatrix<MatrixType>::getUnderlyingMatrix() const
895 {
896  return A_;
897 }
898 
899 template<class MatrixType>
900 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
901 OverlappingRowMatrix<MatrixType>::getExtMatrix() const
902 {
903  return ExtMatrix_;
904 }
905 
906 template<class MatrixType>
907 Teuchos::ArrayView<const size_t> OverlappingRowMatrix<MatrixType>::getExtHaloStarts() const
908 {
909  return ExtHaloStarts_();
910 }
911 
912 
913 } // namespace Ifpack2
914 
915 #define IFPACK2_OVERLAPPINGROWMATRIX_INSTANT(S,LO,GO,N) \
916  template class Ifpack2::OverlappingRowMatrix< Tpetra::RowMatrix<S, LO, GO, N> >;
917 
918 #endif // IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
virtual Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getDomainMap() const
The Map that describes the domain of this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:261
virtual Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getRangeMap() const
The Map that describes the range of this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:277
virtual Teuchos::RCP< const Tpetra::RowGraph< local_ordinal_type, global_ordinal_type, node_type > > getGraph() const
This matrix&#39;s graph.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:285
virtual Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getColMap() const
The Map that describes the distribution of columns over processes.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:252
virtual bool hasTransposeApply() const
Whether this operator&#39;s apply() method can apply the adjoint (transpose).
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:678
virtual mag_type getFrobeniusNorm() const
Returns the Frobenius norm of the matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:611
virtual size_t getNumEntriesInGlobalRow(global_ordinal_type globalRow) const
The number of entries in the given global row that are owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:344
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:233
virtual bool isLocallyIndexed() const
Whether this matrix is locally indexed.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:392
virtual void getGlobalRowView(global_ordinal_type GlobalRow, global_inds_host_view_type &indices, values_host_view_type &values) const
Extract a const, non-persisting view of global indices in a specified row of the matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:483
virtual global_size_t getGlobalNumCols() const
The global number of columns in this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:299
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_OverlappingRowMatrix_def.hpp:594
virtual bool isFillComplete() const
true if fillComplete() has been called, else false.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:406
virtual void getLocalRowCopy(local_ordinal_type LocalRow, nonconst_local_inds_host_view_type &Indices, nonconst_values_host_view_type &Values, size_t &NumEntries) const
Extract a list of entries in a specified local row of the graph. Put into storage allocated by callin...
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:449
virtual bool isGloballyIndexed() const
Whether this matrix is globally indexed.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:399
virtual Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getRowMap() const
The Map that describes the distribution of rows over processes.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:243
virtual bool hasColMap() const
Whether this matrix has a column Map.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:385
virtual size_t getNodeNumEntries() const
The number of entries in this matrix owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:335
virtual void getLocalRowView(local_ordinal_type LocalRow, local_inds_host_view_type &indices, values_host_view_type &values) const
Extract a const, non-persisting view of local indices in a specified row of the matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:524
virtual void getGlobalRowCopy(global_ordinal_type GlobalRow, nonconst_global_inds_host_view_type &Indices, nonconst_values_host_view_type &Values, size_t &NumEntries) const
Extract a list of entries in a specified global row of this matrix. Put into pre-allocated storage...
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:415
virtual bool supportsRowViews() const
true if row views are supported, else false.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:685
Sparse graph (Tpetra::RowGraph subclass) with ghost rows.
Definition: Ifpack2_Details_OverlappingRowGraph_decl.hpp:65
virtual size_t getNodeNumCols() const
The number of columns owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:313
virtual size_t getGlobalMaxNumRowEntries() const
The maximum number of entries in any row on any process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:371
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_OverlappingRowMatrix_def.hpp:603
Definition: Ifpack2_Container_decl.hpp:576
Sparse matrix (Tpetra::RowMatrix subclass) with ghost rows.
Definition: Ifpack2_OverlappingRowMatrix_decl.hpp:58
virtual size_t getNodeMaxNumRowEntries() const
The maximum number of entries in any row on the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:378
OverlappingRowMatrix(const Teuchos::RCP< const row_matrix_type > &A, const int overlapLevel)
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:58
virtual size_t getNodeNumRows() const
The number of rows owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:306
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
Computes the operator-multivector application.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:620
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
virtual global_size_t getGlobalNumEntries() const
The global number of entries in this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:328
virtual size_t getNumEntriesInLocalRow(local_ordinal_type localRow) const
The number of entries in the given local row that are owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:358
virtual global_ordinal_type getIndexBase() const
The index base for global indices for this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:321
virtual global_size_t getGlobalNumRows() const
The global number of rows in this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:292
virtual void getLocalDiagCopy(Tpetra::Vector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &diag) const
Get a copy of the diagonal entries owned by this node, with local row indices.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:561