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