Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
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
55namespace Ifpack2 {
56
57template<class MatrixType>
59OverlappingRowMatrix (const Teuchos::RCP<const row_matrix_type>& A,
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;
73 TEUCHOS_TEST_FOR_EXCEPTION(
74 OverlapLevel_ <= 0, std::runtime_error,
75 "Ifpack2::OverlappingRowMatrix: OverlapLevel must be > 0.");
76 TEUCHOS_TEST_FOR_EXCEPTION
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.");
81 TEUCHOS_TEST_FOR_EXCEPTION(
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 =
90 Teuchos::OrdinalTraits<global_ordinal_type>::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),
143 Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
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 (),
169 Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
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 (),
177 Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
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
233template<class MatrixType>
234Teuchos::RCP<const Teuchos::Comm<int> >
236{
237 return A_->getComm ();
238}
239
240
241
242
243template<class MatrixType>
244Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
246{
247 // FIXME (mfh 12 July 2013) Is this really the right Map to return?
248 return RowMap_;
249}
250
251
252template<class MatrixType>
253Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
255{
256 // FIXME (mfh 12 July 2013) Is this really the right Map to return?
257 return ColMap_;
258}
259
260
261template<class MatrixType>
262Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
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
277template<class MatrixType>
278Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
280{
281 return RowMap_;
282}
283
284
285template<class MatrixType>
286Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
288{
289 return graph_;
290}
291
292
293template<class MatrixType>
295{
296 return NumGlobalRows_;
297}
298
299
300template<class MatrixType>
302{
303 return NumGlobalRows_;
304}
305
306
307template<class MatrixType>
309{
310 return A_->getLocalNumRows () + ExtMatrix_->getLocalNumRows ();
311}
312
313
314template<class MatrixType>
316{
317 return this->getLocalNumRows ();
318}
319
320
321template<class MatrixType>
322typename MatrixType::global_ordinal_type
324{
325 return A_->getIndexBase();
326}
327
328
329template<class MatrixType>
331{
332 return NumGlobalNonzeros_;
333}
334
335
336template<class MatrixType>
338{
339 return A_->getLocalNumEntries () + ExtMatrix_->getLocalNumEntries ();
340}
341
342
343template<class MatrixType>
344size_t
346getNumEntriesInGlobalRow (global_ordinal_type globalRow) const
347{
348 const local_ordinal_type localRow = RowMap_->getLocalElement (globalRow);
349 if (localRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
350 return Teuchos::OrdinalTraits<size_t>::invalid();
351 } else {
352 return getNumEntriesInLocalRow (localRow);
353 }
354}
355
356
357template<class MatrixType>
358size_t
360getNumEntriesInLocalRow (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
372template<class MatrixType>
374{
375 return std::max<size_t>(A_->getGlobalMaxNumRowEntries(), ExtMatrix_->getGlobalMaxNumRowEntries());
376}
377
378
379template<class MatrixType>
381{
382 return MaxNumEntries_;
383}
384
385
386template<class MatrixType>
388{
389 return true;
390}
391
392
393template<class MatrixType>
395{
396 return true;
397}
398
399
400template<class MatrixType>
402{
403 return false;
404}
405
406
407template<class MatrixType>
409{
410 return true;
411}
412
413
414template<class MatrixType>
415void
417 getGlobalRowCopy (global_ordinal_type GlobalRow,
418 nonconst_global_inds_host_view_type &Indices,
419 nonconst_values_host_view_type &Values,
420 size_t& NumEntries) const
421{
422 throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getGlobalRowCopy() not supported.");
423}
424
425template<class MatrixType>
426void
428 getLocalRowCopy (local_ordinal_type LocalRow,
429 nonconst_local_inds_host_view_type &Indices,
430 nonconst_values_host_view_type &Values,
431 size_t& NumEntries) const
432{
433 using Teuchos::as;
434 const size_t numMyRowsA = A_->getLocalNumRows ();
435 if (as<size_t> (LocalRow) < numMyRowsA) {
436 A_->getLocalRowCopy (LocalRow, Indices, Values, NumEntries);
437 } else {
438 ExtMatrix_->getLocalRowCopy (LocalRow - as<local_ordinal_type> (numMyRowsA),
439 Indices, Values, NumEntries);
440 }
441}
442
443template<class MatrixType>
444void
446getGlobalRowView (global_ordinal_type GlobalRow,
447 global_inds_host_view_type &indices,
448 values_host_view_type &values) const {
449 const local_ordinal_type LocalRow = RowMap_->getLocalElement (GlobalRow);
450 if (LocalRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid()) {
451 indices = global_inds_host_view_type();
452 values = values_host_view_type();
453 } else {
454 if (Teuchos::as<size_t> (LocalRow) < A_->getLocalNumRows ()) {
455 A_->getGlobalRowView (GlobalRow, indices, values);
456 } else {
457 ExtMatrix_->getGlobalRowView (GlobalRow, indices, values);
458 }
459 }
460}
461
462template<class MatrixType>
463void
465 getLocalRowView (local_ordinal_type LocalRow,
466 local_inds_host_view_type & indices,
467 values_host_view_type & values) const {
468 using Teuchos::as;
469 const size_t numMyRowsA = A_->getLocalNumRows ();
470 if (as<size_t> (LocalRow) < numMyRowsA) {
471 A_->getLocalRowView (LocalRow, indices, values);
472 } else {
473 ExtMatrix_->getLocalRowView (LocalRow - as<local_ordinal_type> (numMyRowsA),
474 indices, values);
475 }
476
477}
478
479template<class MatrixType>
480void
482getLocalDiagCopy (Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& diag) const
483{
484 using Teuchos::Array;
485
486 //extract diagonal of original matrix
487 vector_type baseDiag(A_->getRowMap()); // diagonal of original matrix A_
488 A_->getLocalDiagCopy(baseDiag);
489 Array<scalar_type> baseDiagVals(baseDiag.getLocalLength());
490 baseDiag.get1dCopy(baseDiagVals());
491 //extra diagonal of ghost matrix
492 vector_type extDiag(ExtMatrix_->getRowMap());
493 ExtMatrix_->getLocalDiagCopy(extDiag);
494 Array<scalar_type> extDiagVals(extDiag.getLocalLength());
495 extDiag.get1dCopy(extDiagVals());
496
497 Teuchos::ArrayRCP<scalar_type> allDiagVals = diag.getDataNonConst();
498 if (allDiagVals.size() != baseDiagVals.size() + extDiagVals.size()) {
499 std::ostringstream errStr;
500 errStr << "Ifpack2::OverlappingRowMatrix::getLocalDiagCopy : Mismatch in diagonal lengths, "
501 << allDiagVals.size() << " != " << baseDiagVals.size() << "+" << extDiagVals.size();
502 throw std::runtime_error(errStr.str());
503 }
504 for (Teuchos::Ordinal i=0; i<baseDiagVals.size(); ++i)
505 allDiagVals[i] = baseDiagVals[i];
506 Teuchos_Ordinal offset=baseDiagVals.size();
507 for (Teuchos::Ordinal i=0; i<extDiagVals.size(); ++i)
508 allDiagVals[i+offset] = extDiagVals[i];
509}
510
511
512template<class MatrixType>
513void
515leftScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& /* x */)
516{
517 throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
518}
519
520
521template<class MatrixType>
522void
524rightScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& /* x */)
525{
526 throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
527}
528
529
530template<class MatrixType>
531typename OverlappingRowMatrix<MatrixType>::mag_type
533{
534 throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support getFrobeniusNorm.");
535}
536
537
538template<class MatrixType>
539void
541apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
542 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
543 Teuchos::ETransp mode,
544 scalar_type alpha,
545 scalar_type beta) const
546{
547 using MV = Tpetra::MultiVector<scalar_type, local_ordinal_type,
548 global_ordinal_type, node_type>;
549 TEUCHOS_TEST_FOR_EXCEPTION
550 (X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
551 "Ifpack2::OverlappingRowMatrix::apply: X.getNumVectors() = "
552 << X.getNumVectors() << " != Y.getNumVectors() = " << Y.getNumVectors()
553 << ".");
554 // If X aliases Y, we'll need to copy X.
555 bool aliases = X.aliases(Y);
556 if (aliases) {
557 MV X_copy (X, Teuchos::Copy);
558 this->apply (X_copy, Y, mode, alpha, beta);
559 return;
560 }
561
562 const auto& rowMap0 = * (A_->getRowMap ());
563 const auto& colMap0 = * (A_->getColMap ());
564 MV X_0 (X, mode == Teuchos::NO_TRANS ? colMap0 : rowMap0, 0);
565 MV Y_0 (Y, mode == Teuchos::NO_TRANS ? rowMap0 : colMap0, 0);
566 A_->localApply (X_0, Y_0, mode, alpha, beta);
567
568 const auto& rowMap1 = * (ExtMatrix_->getRowMap ());
569 const auto& colMap1 = * (ExtMatrix_->getColMap ());
570 MV X_1 (X, mode == Teuchos::NO_TRANS ? colMap1 : rowMap1, 0);
571 MV Y_1 (Y, mode == Teuchos::NO_TRANS ? rowMap1 : colMap1, A_->getLocalNumRows ());
572 ExtMatrix_->localApply (X_1, Y_1, mode, alpha, beta);
573}
574
575
576template<class MatrixType>
577void
579importMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
580 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
581 Tpetra::CombineMode CM)
582{
583 OvX.doImport (X, *Importer_, CM);
584}
585
586
587template<class MatrixType>
588void
589OverlappingRowMatrix<MatrixType>::
590exportMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
591 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
592 Tpetra::CombineMode CM)
593{
594 X.doExport (OvX, *Importer_, CM);
595}
596
597
598template<class MatrixType>
600{
601 return true;
602}
603
604
605template<class MatrixType>
607{
608 return false;
609}
610
611template<class MatrixType>
613{
614 std::ostringstream oss;
615 if (isFillComplete()) {
616 oss << "{ isFillComplete: true"
617 << ", global rows: " << getGlobalNumRows()
618 << ", global columns: " << getGlobalNumCols()
619 << ", global entries: " << getGlobalNumEntries()
620 << " }";
621 }
622 else {
623 oss << "{ isFillComplete: false"
624 << ", global rows: " << getGlobalNumRows()
625 << " }";
626 }
627 return oss.str();
628}
629
630template<class MatrixType>
631void OverlappingRowMatrix<MatrixType>::describe(Teuchos::FancyOStream &out,
632 const Teuchos::EVerbosityLevel verbLevel) const
633{
634 using std::endl;
635 using std::setw;
636 using Teuchos::as;
637 using Teuchos::VERB_DEFAULT;
638 using Teuchos::VERB_NONE;
639 using Teuchos::VERB_LOW;
640 using Teuchos::VERB_MEDIUM;
641 using Teuchos::VERB_HIGH;
642 using Teuchos::VERB_EXTREME;
643 using Teuchos::RCP;
644 using Teuchos::null;
645 using Teuchos::ArrayView;
646
647 Teuchos::EVerbosityLevel vl = verbLevel;
648 if (vl == VERB_DEFAULT) {
649 vl = VERB_LOW;
650 }
651 RCP<const Teuchos::Comm<int> > comm = this->getComm();
652 const int myRank = comm->getRank();
653 const int numProcs = comm->getSize();
654 size_t width = 1;
655 for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
656 ++width;
657 }
658 width = std::max<size_t> (width, as<size_t> (11)) + 2;
659 Teuchos::OSTab tab(out);
660 // none: print nothing
661 // low: print O(1) info from node 0
662 // medium: print O(P) info, num entries per process
663 // high: print O(N) info, num entries per row
664 // extreme: print O(NNZ) info: print indices and values
665 //
666 // for medium and higher, print constituent objects at specified verbLevel
667 if (vl != VERB_NONE) {
668 if (myRank == 0) {
669 out << this->description() << std::endl;
670 }
671 // O(1) globals, minus what was already printed by description()
672 //if (isFillComplete() && myRank == 0) {
673 // out << "Global max number of entries in a row: " << getGlobalMaxNumRowEntries() << std::endl;
674 //}
675 // constituent objects
676 if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
677 if (myRank == 0) {
678 out << endl << "Row map:" << endl;
679 }
680 getRowMap()->describe(out,vl);
681 //
682 if (getColMap() != null) {
683 if (getColMap() == getRowMap()) {
684 if (myRank == 0) {
685 out << endl << "Column map is row map.";
686 }
687 }
688 else {
689 if (myRank == 0) {
690 out << endl << "Column map:" << endl;
691 }
692 getColMap()->describe(out,vl);
693 }
694 }
695 if (getDomainMap() != null) {
696 if (getDomainMap() == getRowMap()) {
697 if (myRank == 0) {
698 out << endl << "Domain map is row map.";
699 }
700 }
701 else if (getDomainMap() == getColMap()) {
702 if (myRank == 0) {
703 out << endl << "Domain map is column map.";
704 }
705 }
706 else {
707 if (myRank == 0) {
708 out << endl << "Domain map:" << endl;
709 }
710 getDomainMap()->describe(out,vl);
711 }
712 }
713 if (getRangeMap() != null) {
714 if (getRangeMap() == getDomainMap()) {
715 if (myRank == 0) {
716 out << endl << "Range map is domain map." << endl;
717 }
718 }
719 else if (getRangeMap() == getRowMap()) {
720 if (myRank == 0) {
721 out << endl << "Range map is row map." << endl;
722 }
723 }
724 else {
725 if (myRank == 0) {
726 out << endl << "Range map: " << endl;
727 }
728 getRangeMap()->describe(out,vl);
729 }
730 }
731 if (myRank == 0) {
732 out << endl;
733 }
734 }
735 // O(P) data
736 if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
737 for (int curRank = 0; curRank < numProcs; ++curRank) {
738 if (myRank == curRank) {
739 out << "Process rank: " << curRank << std::endl;
740 out << " Number of entries: " << getLocalNumEntries() << std::endl;
741 out << " Max number of entries per row: " << getLocalMaxNumRowEntries() << std::endl;
742 }
743 comm->barrier();
744 comm->barrier();
745 comm->barrier();
746 }
747 }
748 // O(N) and O(NNZ) data
749 if (vl == VERB_HIGH || vl == VERB_EXTREME) {
750 for (int curRank = 0; curRank < numProcs; ++curRank) {
751 if (myRank == curRank) {
752 out << std::setw(width) << "Proc Rank"
753 << std::setw(width) << "Global Row"
754 << std::setw(width) << "Num Entries";
755 if (vl == VERB_EXTREME) {
756 out << std::setw(width) << "(Index,Value)";
757 }
758 out << endl;
759 for (size_t r = 0; r < getLocalNumRows (); ++r) {
760 const size_t nE = getNumEntriesInLocalRow(r);
761 typename MatrixType::global_ordinal_type gid = getRowMap()->getGlobalElement(r);
762 out << std::setw(width) << myRank
763 << std::setw(width) << gid
764 << std::setw(width) << nE;
765 if (vl == VERB_EXTREME) {
766 if (isGloballyIndexed()) {
767 global_inds_host_view_type rowinds;
768 values_host_view_type rowvals;
769 getGlobalRowView (gid, rowinds, rowvals);
770 for (size_t j = 0; j < nE; ++j) {
771 out << " (" << rowinds[j]
772 << ", " << rowvals[j]
773 << ") ";
774 }
775 }
776 else if (isLocallyIndexed()) {
777 local_inds_host_view_type rowinds;
778 values_host_view_type rowvals;
779 getLocalRowView (r, rowinds, rowvals);
780 for (size_t j=0; j < nE; ++j) {
781 out << " (" << getColMap()->getGlobalElement(rowinds[j])
782 << ", " << rowvals[j]
783 << ") ";
784 }
785 } // globally or locally indexed
786 } // vl == VERB_EXTREME
787 out << endl;
788 } // for each row r on this process
789
790 } // if (myRank == curRank)
791 comm->barrier();
792 comm->barrier();
793 comm->barrier();
794 }
795
796 out.setOutputToRootOnly(0);
797 out << "===========\nlocal matrix\n=================" << std::endl;
798 out.setOutputToRootOnly(-1);
799 A_->describe(out,Teuchos::VERB_EXTREME);
800 out.setOutputToRootOnly(0);
801 out << "===========\nend of local matrix\n=================" << std::endl;
802 comm->barrier();
803 out.setOutputToRootOnly(0);
804 out << "=================\nghost matrix\n=================" << std::endl;
805 out.setOutputToRootOnly(-1);
806 ExtMatrix_->describe(out,Teuchos::VERB_EXTREME);
807 out.setOutputToRootOnly(0);
808 out << "===========\nend of ghost matrix\n=================" << std::endl;
809 }
810 }
811}
812
813template<class MatrixType>
814Teuchos::RCP<const Tpetra::CrsMatrix<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
815OverlappingRowMatrix<MatrixType>::getUnderlyingMatrix() const
816{
817 return A_;
818}
819
820template<class MatrixType>
821Teuchos::RCP<const Tpetra::CrsMatrix<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
822OverlappingRowMatrix<MatrixType>::getExtMatrix() const
823{
824 return ExtMatrix_;
825}
826
827template<class MatrixType>
828Kokkos::View<size_t*, typename OverlappingRowMatrix<MatrixType>::device_type>
829OverlappingRowMatrix<MatrixType>::getExtHaloStarts() const
830{
831 return ExtHaloStarts_;
832}
833
834template<class MatrixType>
835typename Kokkos::View<size_t*, typename OverlappingRowMatrix<MatrixType>::device_type>::HostMirror
836OverlappingRowMatrix<MatrixType>::getExtHaloStartsHost() const
837{
838 return ExtHaloStarts_h;
839}
840
841template<class MatrixType>
842void OverlappingRowMatrix<MatrixType>::doExtImport()
843{
844 //TODO: CrsMatrix can't doImport after resumeFill (see #9720). Ideally, this import could
845 //happen using combine mode REPLACE without reconstructing the matrix.
846 //Maybe even without another fillComplete since this doesn't change structure - see #9655.
847 ExtMatrix_ = rcp (new crs_matrix_type (ExtMap_, ColMap_, 0));
848 ExtMatrix_->doImport (*A_, *ExtImporter_, Tpetra::INSERT);
849 ExtMatrix_->fillComplete (A_->getDomainMap (), RowMap_);
850}
851
852} // namespace Ifpack2
853
854#define IFPACK2_OVERLAPPINGROWMATRIX_INSTANT(S,LO,GO,N) \
855 template class Ifpack2::OverlappingRowMatrix< Tpetra::RowMatrix<S, LO, GO, N> >;
856
857#endif // IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
Sparse graph (Tpetra::RowGraph subclass) with ghost rows.
Definition Ifpack2_Details_OverlappingRowGraph_decl.hpp:66
Sparse matrix (Tpetra::RowMatrix subclass) with ghost rows.
Definition Ifpack2_OverlappingRowMatrix_decl.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:541
virtual bool hasTransposeApply() const
Whether this operator's apply() method can apply the adjoint (transpose).
Definition Ifpack2_OverlappingRowMatrix_def.hpp:599
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 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:428
virtual bool isLocallyIndexed() const
Whether this matrix is locally indexed.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:394
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 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:465
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:417
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 size_t getGlobalMaxNumRowEntries() const
The maximum number of entries in any row on any process.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:373
virtual global_size_t getGlobalNumCols() const
The global number of columns in this matrix.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:301
virtual bool isGloballyIndexed() const
Whether this matrix is globally indexed.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:401
virtual global_ordinal_type getIndexBase() const
The index base for global indices for this matrix.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:323
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::RowGraph< local_ordinal_type, global_ordinal_type, node_type > > getGraph() const
This matrix's graph.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:287
virtual size_t getLocalNumCols() const
The number of columns owned by the calling process.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:315
virtual size_t getLocalNumRows() const
The number of rows owned by the calling process.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:308
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:446
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:524
virtual bool hasColMap() const
Whether this matrix has a column Map.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:387
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:235
virtual global_size_t getGlobalNumRows() const
The global number of rows in this matrix.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:294
virtual bool supportsRowViews() const
true if row views are supported, else false.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:606
virtual global_size_t getGlobalNumEntries() const
The global number of entries in this matrix.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:330
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 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:515
OverlappingRowMatrix(const Teuchos::RCP< const row_matrix_type > &A, const int overlapLevel)
Definition Ifpack2_OverlappingRowMatrix_def.hpp:59
virtual size_t getLocalMaxNumRowEntries() const
The maximum number of entries in any row on the calling process.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:380
virtual bool isFillComplete() const
true if fillComplete() has been called, else false.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:408
virtual mag_type getFrobeniusNorm() const
Returns the Frobenius norm of the matrix.
Definition Ifpack2_OverlappingRowMatrix_def.hpp:532
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:482
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
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74