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