Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_OverlappingRowMatrix_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #ifndef IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
44 #define IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
45 
46 #include <Ifpack2_OverlappingRowMatrix_decl.hpp>
47 #include <Ifpack2_Details_OverlappingRowGraph.hpp>
48 #include <Tpetra_CrsMatrix.hpp>
49 #include <Teuchos_CommHelpers.hpp>
50 
51 namespace Ifpack2 {
52 
53 template<class MatrixType>
55 OverlappingRowMatrix (const Teuchos::RCP<const row_matrix_type>& A,
56  const int overlapLevel) :
57  A_ (A),
58  OverlapLevel_ (overlapLevel)
59 {
60  using Teuchos::RCP;
61  using Teuchos::rcp;
62  using Teuchos::Array;
63  using Teuchos::outArg;
64  using Teuchos::rcp_const_cast;
65  using Teuchos::rcp_dynamic_cast;
66  using Teuchos::rcp_implicit_cast;
67  using Teuchos::REDUCE_SUM;
68  using Teuchos::reduceAll;
69  typedef Tpetra::global_size_t GST;
70  typedef Tpetra::CrsGraph<local_ordinal_type,
71  global_ordinal_type, node_type> crs_graph_type;
72  TEUCHOS_TEST_FOR_EXCEPTION(
73  OverlapLevel_ <= 0, std::runtime_error,
74  "Ifpack2::OverlappingRowMatrix: OverlapLevel must be > 0.");
75  TEUCHOS_TEST_FOR_EXCEPTION(
76  A_->getComm()->getSize() == 1, std::runtime_error,
77  "Ifpack2::OverlappingRowMatrix: Matrix must be "
78  "distributed over more than one MPI process.");
79 
80  RCP<const crs_matrix_type> ACRS =
81  rcp_dynamic_cast<const crs_matrix_type, const row_matrix_type> (A_);
82  TEUCHOS_TEST_FOR_EXCEPTION(
83  ACRS.is_null (), std::runtime_error,
84  "Ifpack2::OverlappingRowMatrix: The input matrix must be a Tpetra::"
85  "CrsMatrix with matching template parameters. This class currently "
86  "requires that CrsMatrix's fifth template parameter be the default.");
87  RCP<const crs_graph_type> A_crsGraph = ACRS->getCrsGraph ();
88 
89  const size_t numMyRowsA = A_->getNodeNumRows ();
90  const global_ordinal_type global_invalid =
91  Teuchos::OrdinalTraits<global_ordinal_type>::invalid ();
92 
93  // Temp arrays
94  Array<global_ordinal_type> ExtElements;
95  RCP<map_type> TmpMap;
96  RCP<crs_graph_type> TmpGraph;
97  RCP<import_type> TmpImporter;
98  RCP<const map_type> RowMap, ColMap;
99 
100  // The big import loop
101  for (int overlap = 0 ; overlap < OverlapLevel_ ; ++overlap) {
102  // Get the current maps
103  if (overlap == 0) {
104  RowMap = A_->getRowMap ();
105  ColMap = A_->getColMap ();
106  }
107  else {
108  RowMap = TmpGraph->getRowMap ();
109  ColMap = TmpGraph->getColMap ();
110  }
111 
112  const size_t size = ColMap->getNodeNumElements () - RowMap->getNodeNumElements ();
113  Array<global_ordinal_type> mylist (size);
114  size_t count = 0;
115 
116  // define the set of rows that are in ColMap but not in RowMap
117  for (local_ordinal_type i = 0 ; (size_t) i < ColMap->getNodeNumElements() ; ++i) {
118  const global_ordinal_type GID = ColMap->getGlobalElement (i);
119  if (A_->getRowMap ()->getLocalElement (GID) == global_invalid) {
120  typedef typename Array<global_ordinal_type>::iterator iter_type;
121  const iter_type end = ExtElements.end ();
122  const iter_type pos = std::find (ExtElements.begin (), end, GID);
123  if (pos == end) {
124  ExtElements.push_back (GID);
125  mylist[count] = GID;
126  ++count;
127  }
128  }
129  }
130 
131  // mfh 24 Nov 2013: We don't need TmpMap, TmpGraph, or
132  // TmpImporter after this loop, so we don't have to construct them
133  // on the last round.
134  if (overlap + 1 < OverlapLevel_) {
135  // Allocate & import new matrices, maps, etc.
136  //
137  // FIXME (mfh 24 Nov 2013) Do we always want to use index base
138  // zero? It doesn't really matter, since the actual index base
139  // (in the current implementation of Map) will always be the
140  // globally least GID.
141  TmpMap = rcp (new map_type (global_invalid, mylist (0, count),
142  Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
143  A_->getComm (), A_->getNode ()));
144  TmpGraph = rcp (new crs_graph_type (TmpMap, 0));
145  TmpImporter = rcp (new import_type (A_->getRowMap (), TmpMap));
146 
147  TmpGraph->doImport (*A_crsGraph, *TmpImporter, Tpetra::INSERT);
148  TmpGraph->fillComplete (A_->getDomainMap (), TmpMap);
149  }
150  }
151 
152  // build the map containing all the nodes (original
153  // matrix + extended matrix)
154  Array<global_ordinal_type> mylist (numMyRowsA + ExtElements.size ());
155  for (local_ordinal_type i = 0; (size_t)i < numMyRowsA; ++i) {
156  mylist[i] = A_->getRowMap ()->getGlobalElement (i);
157  }
158  for (local_ordinal_type i = 0; i < ExtElements.size (); ++i) {
159  mylist[i + numMyRowsA] = ExtElements[i];
160  }
161 
162  RowMap_ = rcp (new map_type (global_invalid, mylist (),
163  Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
164  A_->getComm (), A_->getNode ()));
165  ColMap_ = RowMap_;
166 
167  // now build the map corresponding to all the external nodes
168  // (with respect to A().RowMatrixRowMap().
169  ExtMap_ = rcp (new map_type (global_invalid, ExtElements (),
170  Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
171  A_->getComm (), A_->getNode ()));
172  ExtMatrix_ = rcp (new crs_matrix_type (ExtMap_, ColMap_, 0));
173  ExtImporter_ = rcp (new import_type (A_->getRowMap (), ExtMap_));
174 
175  RCP<crs_matrix_type> ExtMatrixCRS =
176  rcp_dynamic_cast<crs_matrix_type, row_matrix_type> (ExtMatrix_);
177  ExtMatrixCRS->doImport (*ACRS, *ExtImporter_, Tpetra::INSERT);
178  ExtMatrixCRS->fillComplete (A_->getDomainMap (), RowMap_);
179 
180  Importer_ = rcp (new import_type (A_->getRowMap (), RowMap_));
181 
182  // fix indices for overlapping matrix
183  const size_t numMyRowsB = ExtMatrix_->getNodeNumRows ();
184 
185  GST NumMyNonzeros_tmp = A_->getNodeNumEntries () + ExtMatrix_->getNodeNumEntries ();
186  GST NumMyRows_tmp = numMyRowsA + numMyRowsB;
187  {
188  GST inArray[2], outArray[2];
189  inArray[0] = NumMyNonzeros_tmp;
190  inArray[1] = NumMyRows_tmp;
191  outArray[0] = 0;
192  outArray[1] = 0;
193  reduceAll<int, GST> (* (A_->getComm ()), REDUCE_SUM, 2, inArray, outArray);
194  NumGlobalNonzeros_ = outArray[0];
195  NumGlobalRows_ = outArray[1];
196  }
197  // reduceAll<int, GST> (* (A_->getComm ()), REDUCE_SUM, NumMyNonzeros_tmp,
198  // outArg (NumGlobalNonzeros_));
199  // reduceAll<int, GST> (* (A_->getComm ()), REDUCE_SUM, NumMyRows_tmp,
200  // outArg (NumGlobalRows_));
201 
202  MaxNumEntries_ = A_->getNodeMaxNumRowEntries ();
203  if (MaxNumEntries_ < ExtMatrix_->getNodeMaxNumRowEntries ()) {
204  MaxNumEntries_ = ExtMatrix_->getNodeMaxNumRowEntries ();
205  }
206 
207  // Create the graph (returned by getGraph()).
208  typedef Details::OverlappingRowGraph<row_graph_type> row_graph_impl_type;
209  RCP<row_graph_impl_type> graph =
210  rcp (new row_graph_impl_type (A_->getGraph (),
211  ExtMatrix_->getGraph (),
212  RowMap_,
213  ColMap_,
214  NumGlobalRows_,
215  NumGlobalRows_, // # global cols == # global rows
216  NumGlobalNonzeros_,
217  MaxNumEntries_,
218  rcp_const_cast<const import_type> (Importer_),
219  rcp_const_cast<const import_type> (ExtImporter_)));
220  graph_ = rcp_const_cast<const row_graph_type> (rcp_implicit_cast<row_graph_type> (graph));
221  // Resize temp arrays
222  Indices_.resize (MaxNumEntries_);
223  Values_.resize (MaxNumEntries_);
224 }
225 
226 
227 template<class MatrixType>
229 
230 
231 template<class MatrixType>
232 Teuchos::RCP<const Teuchos::Comm<int> >
234 {
235  return A_->getComm ();
236 }
237 
238 
239 template<class MatrixType>
240 Teuchos::RCP<typename MatrixType::node_type>
242 {
243  return A_->getNode();
244 }
245 
246 
247 template<class MatrixType>
248 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
250 {
251  // FIXME (mfh 12 July 2013) Is this really the right Map to return?
252  return RowMap_;
253 }
254 
255 
256 template<class MatrixType>
257 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
259 {
260  // FIXME (mfh 12 July 2013) Is this really the right Map to return?
261  return ColMap_;
262 }
263 
264 
265 template<class MatrixType>
266 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
268 {
269  // The original matrix's domain map is irrelevant; we want the map associated
270  // with the overlap. This can then be used by LocalFilter, for example, while
271  // letting LocalFilter still filter based on domain and range maps (instead of
272  // column and row maps).
273  // FIXME Ideally, this would be the same map but restricted to a local
274  // communicator. If replaceCommWithSubset were free, that would be the way to
275  // go. That would require a new Map ctor. For now, we'll stick with ColMap_'s
276  // global communicator.
277  return ColMap_;
278 }
279 
280 
281 template<class MatrixType>
282 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
284 {
285  return RowMap_;
286 }
287 
288 
289 template<class MatrixType>
290 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
292 {
293  return graph_;
294 }
295 
296 
297 template<class MatrixType>
299 {
300  return NumGlobalRows_;
301 }
302 
303 
304 template<class MatrixType>
306 {
307  return NumGlobalRows_;
308 }
309 
310 
311 template<class MatrixType>
313 {
314  return A_->getNodeNumRows () + ExtMatrix_->getNodeNumRows ();
315 }
316 
317 
318 template<class MatrixType>
320 {
321  return this->getNodeNumRows ();
322 }
323 
324 
325 template<class MatrixType>
326 typename MatrixType::global_ordinal_type
328 {
329  return A_->getIndexBase();
330 }
331 
332 
333 template<class MatrixType>
335 {
336  return NumGlobalNonzeros_;
337 }
338 
339 
340 template<class MatrixType>
342 {
343  return A_->getNodeNumEntries () + ExtMatrix_->getNodeNumEntries ();
344 }
345 
346 
347 template<class MatrixType>
348 size_t
350 getNumEntriesInGlobalRow (global_ordinal_type globalRow) const
351 {
352  const local_ordinal_type localRow = RowMap_->getLocalElement (globalRow);
353  if (localRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
354  return Teuchos::OrdinalTraits<size_t>::invalid();
355  } else {
356  return getNumEntriesInLocalRow (localRow);
357  }
358 }
359 
360 
361 template<class MatrixType>
362 size_t
364 getNumEntriesInLocalRow (local_ordinal_type localRow) const
365 {
366  using Teuchos::as;
367  const size_t numMyRowsA = A_->getNodeNumRows ();
368  if (as<size_t> (localRow) < numMyRowsA) {
369  return A_->getNumEntriesInLocalRow (localRow);
370  } else {
371  return ExtMatrix_->getNumEntriesInLocalRow (as<local_ordinal_type> (localRow - numMyRowsA));
372  }
373 }
374 
375 
376 template<class MatrixType>
378 {
379  throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getGlobalNumDiags() not supported.");
380 }
381 
382 
383 template<class MatrixType>
385 {
386  return A_->getNodeNumDiags();
387 }
388 
389 
390 template<class MatrixType>
392 {
393  throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getGlobalMaxNumRowEntries() not supported.");
394 }
395 
396 
397 template<class MatrixType>
399 {
400  return MaxNumEntries_;
401 }
402 
403 
404 template<class MatrixType>
406 {
407  return true;
408 }
409 
410 
411 template<class MatrixType>
413 {
414  return A_->isLowerTriangular();
415 }
416 
417 
418 template<class MatrixType>
420 {
421  return A_->isUpperTriangular();
422 }
423 
424 
425 template<class MatrixType>
427 {
428  return true;
429 }
430 
431 
432 template<class MatrixType>
434 {
435  return false;
436 }
437 
438 
439 template<class MatrixType>
441 {
442  return true;
443 }
444 
445 
446 template<class MatrixType>
447 void
449 getGlobalRowCopy (global_ordinal_type GlobalRow,
450  const Teuchos::ArrayView<global_ordinal_type> &Indices,
451  const Teuchos::ArrayView<scalar_type>& Values,
452  size_t& NumEntries) const
453 {
454  const local_ordinal_type LocalRow = RowMap_->getLocalElement (GlobalRow);
455  if (LocalRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
456  NumEntries = Teuchos::OrdinalTraits<size_t>::invalid ();
457  } else {
458  if (Teuchos::as<size_t> (LocalRow) < A_->getNodeNumRows ()) {
459  A_->getGlobalRowCopy (GlobalRow, Indices, Values, NumEntries);
460  } else {
461  ExtMatrix_->getGlobalRowCopy (GlobalRow, Indices, Values, NumEntries);
462  }
463  }
464 }
465 
466 
467 template<class MatrixType>
468 void
470 getLocalRowCopy (local_ordinal_type LocalRow,
471  const Teuchos::ArrayView<local_ordinal_type> &Indices,
472  const Teuchos::ArrayView<scalar_type> &Values,
473  size_t &NumEntries) const
474 {
475  using Teuchos::as;
476  const size_t numMyRowsA = A_->getNodeNumRows ();
477  if (as<size_t> (LocalRow) < numMyRowsA) {
478  A_->getLocalRowCopy (LocalRow, Indices, Values, NumEntries);
479  } else {
480  ExtMatrix_->getLocalRowCopy (LocalRow - as<local_ordinal_type> (numMyRowsA),
481  Indices, Values, NumEntries);
482  }
483 }
484 
485 
486 template<class MatrixType>
487 void
489 getGlobalRowView (global_ordinal_type GlobalRow,
490  Teuchos::ArrayView<const global_ordinal_type>& indices,
491  Teuchos::ArrayView<const scalar_type>& values) const
492 {
493  const local_ordinal_type LocalRow = RowMap_->getLocalElement (GlobalRow);
494  if (LocalRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid()) {
495  indices = Teuchos::null;
496  values = Teuchos::null;
497  } else {
498  if (Teuchos::as<size_t> (LocalRow) < A_->getNodeNumRows ()) {
499  A_->getGlobalRowView (GlobalRow, indices, values);
500  } else {
501  ExtMatrix_->getGlobalRowView (GlobalRow, indices, values);
502  }
503  }
504 }
505 
506 
507 template<class MatrixType>
508 void
510 getLocalRowView (local_ordinal_type LocalRow,
511  Teuchos::ArrayView<const local_ordinal_type>& indices,
512  Teuchos::ArrayView<const scalar_type>& values) const
513 {
514  using Teuchos::as;
515  const size_t numMyRowsA = A_->getNodeNumRows ();
516  if (as<size_t> (LocalRow) < numMyRowsA) {
517  A_->getLocalRowView (LocalRow, indices, values);
518  } else {
519  ExtMatrix_->getLocalRowView (LocalRow - as<local_ordinal_type> (numMyRowsA),
520  indices, values);
521  }
522 }
523 
524 
525 template<class MatrixType>
526 void
528 getLocalDiagCopy (Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& diag) const
529 {
530  throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getLocalDiagCopy not supported.");
531 }
532 
533 
534 template<class MatrixType>
535 void
537 leftScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
538 {
539  throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
540 }
541 
542 
543 template<class MatrixType>
544 void
546 rightScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
547 {
548  throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
549 }
550 
551 
552 template<class MatrixType>
553 typename OverlappingRowMatrix<MatrixType>::mag_type
555 {
556  throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support getFrobeniusNorm.");
557 }
558 
559 
560 template<class MatrixType>
561 void
563 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
564  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
565  Teuchos::ETransp mode,
566  scalar_type alpha,
567  scalar_type beta) const
568 {
569  using Teuchos::ArrayRCP;
570  using Teuchos::as;
571  typedef scalar_type RangeScalar;
572  typedef scalar_type DomainScalar;
573  typedef Teuchos::ScalarTraits<RangeScalar> STRS;
574 
575  TEUCHOS_TEST_FOR_EXCEPTION(
576  alpha != Teuchos::ScalarTraits<scalar_type>::one () ||
577  beta != Teuchos::ScalarTraits<scalar_type>::zero (), std::logic_error,
578  "Ifpack2::ReorderFilter::apply is only implemented for alpha = 1 and "
579  "beta = 0. You set alpha = " << alpha << " and beta = " << beta << ".");
580  TEUCHOS_TEST_FOR_EXCEPTION(
581  X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
582  "Ifpack2::OverlappingRowMatrix::apply: The input X and the output Y must "
583  "have the same number of columns. X.getNumVectors() = "
584  << X.getNumVectors() << " != Y.getNumVectors() = " << Y.getNumVectors()
585  << ".");
586 
587  // FIXME (mfh 13 July 2013) This would be a good candidate for a
588  // local parallel operator implementation. That would obviate the
589  // need for getting views of the data and make the code below a lot
590  // simpler.
591 
592  const RangeScalar zero = STRS::zero ();
593  ArrayRCP<ArrayRCP<const DomainScalar> > x_ptr = X.get2dView();
594  ArrayRCP<ArrayRCP<RangeScalar> > y_ptr = Y.get2dViewNonConst();
595  Y.putScalar(zero);
596  size_t NumVectors = Y.getNumVectors();
597 
598  const size_t numMyRowsA = A_->getNodeNumRows ();
599  for (size_t i = 0; i < numMyRowsA; ++i) {
600  size_t Nnz;
601  // Use this class's getrow to make the below code simpler
602  A_->getLocalRowCopy (i, Indices_ (),Values_ (), Nnz);
603  if (mode == Teuchos::NO_TRANS) {
604  for (size_t j = 0; j < Nnz; ++j)
605  for (size_t k = 0; k < NumVectors; ++k)
606  y_ptr[k][i] += as<RangeScalar> (Values_[j]) *
607  as<RangeScalar> (x_ptr[k][Indices_[j]]);
608  }
609  else if (mode == Teuchos::TRANS){
610  for (size_t j = 0; j < Nnz; ++j)
611  for (size_t k = 0; k < NumVectors; ++k)
612  y_ptr[k][Indices_[j]] += as<RangeScalar> (Values_[j]) *
613  as<RangeScalar> (x_ptr[k][i]);
614  }
615  else { // mode == Teuchos::CONJ_TRANS
616  for (size_t j = 0; j < Nnz; ++j)
617  for (size_t k = 0; k < NumVectors; ++k)
618  y_ptr[k][Indices_[j]] +=
619  STRS::conjugate (as<RangeScalar> (Values_[j])) *
620  as<RangeScalar> (x_ptr[k][i]);
621  }
622  }
623 
624  const size_t numMyRowsB = ExtMatrix_->getNodeNumRows ();
625  for (size_t i = 0 ; i < numMyRowsB ; ++i) {
626  size_t Nnz;
627  // Use this class's getrow to make the below code simpler
628  ExtMatrix_->getLocalRowCopy (i, Indices_ (), Values_ (), Nnz);
629  if (mode == Teuchos::NO_TRANS) {
630  for (size_t j = 0; j < Nnz; ++j)
631  for (size_t k = 0; k < NumVectors; ++k)
632  y_ptr[k][numMyRowsA+i] += as<RangeScalar> (Values_[j]) *
633  as<RangeScalar> (x_ptr[k][Indices_[j]]);
634  }
635  else if (mode == Teuchos::TRANS) {
636  for (size_t j = 0; j < Nnz; ++j)
637  for (size_t k = 0; k < NumVectors; ++k)
638  y_ptr[k][numMyRowsA+Indices_[j]] += as<RangeScalar> (Values_[j]) *
639  as<RangeScalar> (x_ptr[k][i]);
640  }
641  else { // mode == Teuchos::CONJ_TRANS
642  for (size_t j = 0; j < Nnz; ++j)
643  for (size_t k = 0; k < NumVectors; ++k)
644  y_ptr[k][numMyRowsA+Indices_[j]] +=
645  STRS::conjugate (as<RangeScalar> (Values_[j])) *
646  as<RangeScalar> (x_ptr[k][i]);
647  }
648  }
649 }
650 
651 
652 template<class MatrixType>
653 void
655 importMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
656  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
657  Tpetra::CombineMode CM)
658 {
659  OvX.doImport (X, *Importer_, CM);
660 }
661 
662 
663 template<class MatrixType>
664 void
665 OverlappingRowMatrix<MatrixType>::
666 exportMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
667  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
668  Tpetra::CombineMode CM)
669 {
670  X.doExport (OvX, *Importer_, CM);
671 }
672 
673 
674 template<class MatrixType>
676 {
677  return true;
678 }
679 
680 
681 template<class MatrixType>
683 {
684  return false;
685 }
686 
687 template<class MatrixType>
689 {
690  std::ostringstream oss;
691  if (isFillComplete()) {
692  oss << "{ isFillComplete: true"
693  << ", global rows: " << getGlobalNumRows()
694  << ", global columns: " << getGlobalNumCols()
695  << ", global entries: " << getGlobalNumEntries()
696  << " }";
697  }
698  else {
699  oss << "{ isFillComplete: false"
700  << ", global rows: " << getGlobalNumRows()
701  << " }";
702  }
703  return oss.str();
704 }
705 
706 template<class MatrixType>
707 void OverlappingRowMatrix<MatrixType>::describe(Teuchos::FancyOStream &out,
708  const Teuchos::EVerbosityLevel verbLevel) const
709 {
710  using std::endl;
711  using std::setw;
712  using Teuchos::as;
713  using Teuchos::VERB_DEFAULT;
714  using Teuchos::VERB_NONE;
715  using Teuchos::VERB_LOW;
716  using Teuchos::VERB_MEDIUM;
717  using Teuchos::VERB_HIGH;
718  using Teuchos::VERB_EXTREME;
719  using Teuchos::RCP;
720  using Teuchos::null;
721  using Teuchos::ArrayView;
722 
723  Teuchos::EVerbosityLevel vl = verbLevel;
724  if (vl == VERB_DEFAULT) {
725  vl = VERB_LOW;
726  }
727  RCP<const Teuchos::Comm<int> > comm = this->getComm();
728  const int myRank = comm->getRank();
729  const int numProcs = comm->getSize();
730  size_t width = 1;
731  for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
732  ++width;
733  }
734  width = std::max<size_t> (width, as<size_t> (11)) + 2;
735  Teuchos::OSTab tab(out);
736  // none: print nothing
737  // low: print O(1) info from node 0
738  // medium: print O(P) info, num entries per process
739  // high: print O(N) info, num entries per row
740  // extreme: print O(NNZ) info: print indices and values
741  //
742  // for medium and higher, print constituent objects at specified verbLevel
743  if (vl != VERB_NONE) {
744  if (myRank == 0) {
745  out << this->description() << std::endl;
746  }
747  // O(1) globals, minus what was already printed by description()
748  //if (isFillComplete() && myRank == 0) {
749  // out << "Global number of diagonal entries: " << getGlobalNumDiags() << std::endl;
750  // out << "Global max number of entries in a row: " << getGlobalMaxNumRowEntries() << std::endl;
751  //}
752  // constituent objects
753  if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
754  if (myRank == 0) {
755  out << endl << "Row map:" << endl;
756  }
757  getRowMap()->describe(out,vl);
758  //
759  if (getColMap() != null) {
760  if (getColMap() == getRowMap()) {
761  if (myRank == 0) {
762  out << endl << "Column map is row map.";
763  }
764  }
765  else {
766  if (myRank == 0) {
767  out << endl << "Column map:" << endl;
768  }
769  getColMap()->describe(out,vl);
770  }
771  }
772  if (getDomainMap() != null) {
773  if (getDomainMap() == getRowMap()) {
774  if (myRank == 0) {
775  out << endl << "Domain map is row map.";
776  }
777  }
778  else if (getDomainMap() == getColMap()) {
779  if (myRank == 0) {
780  out << endl << "Domain map is column map.";
781  }
782  }
783  else {
784  if (myRank == 0) {
785  out << endl << "Domain map:" << endl;
786  }
787  getDomainMap()->describe(out,vl);
788  }
789  }
790  if (getRangeMap() != null) {
791  if (getRangeMap() == getDomainMap()) {
792  if (myRank == 0) {
793  out << endl << "Range map is domain map." << endl;
794  }
795  }
796  else if (getRangeMap() == getRowMap()) {
797  if (myRank == 0) {
798  out << endl << "Range map is row map." << endl;
799  }
800  }
801  else {
802  if (myRank == 0) {
803  out << endl << "Range map: " << endl;
804  }
805  getRangeMap()->describe(out,vl);
806  }
807  }
808  if (myRank == 0) {
809  out << endl;
810  }
811  }
812  // O(P) data
813  if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
814  for (int curRank = 0; curRank < numProcs; ++curRank) {
815  if (myRank == curRank) {
816  out << "Process rank: " << curRank << std::endl;
817  out << " Number of entries: " << getNodeNumEntries() << std::endl;
818  if (isFillComplete()) {
819  out << " Number of diagonal entries: " << getNodeNumDiags() << std::endl;
820  }
821  out << " Max number of entries per row: " << getNodeMaxNumRowEntries() << std::endl;
822  }
823  comm->barrier();
824  comm->barrier();
825  comm->barrier();
826  }
827  }
828  // O(N) and O(NNZ) data
829  if (vl == VERB_HIGH || vl == VERB_EXTREME) {
830  for (int curRank = 0; curRank < numProcs; ++curRank) {
831  if (myRank == curRank) {
832  out << std::setw(width) << "Proc Rank"
833  << std::setw(width) << "Global Row"
834  << std::setw(width) << "Num Entries";
835  if (vl == VERB_EXTREME) {
836  out << std::setw(width) << "(Index,Value)";
837  }
838  out << endl;
839  for (size_t r = 0; r < getNodeNumRows (); ++r) {
840  const size_t nE = getNumEntriesInLocalRow(r);
841  typename MatrixType::global_ordinal_type gid = getRowMap()->getGlobalElement(r);
842  out << std::setw(width) << myRank
843  << std::setw(width) << gid
844  << std::setw(width) << nE;
845  if (vl == VERB_EXTREME) {
846  if (isGloballyIndexed()) {
847  ArrayView<const typename MatrixType::global_ordinal_type> rowinds;
848  ArrayView<const typename MatrixType::scalar_type> rowvals;
849  getGlobalRowView (gid, rowinds, rowvals);
850  for (size_t j = 0; j < nE; ++j) {
851  out << " (" << rowinds[j]
852  << ", " << rowvals[j]
853  << ") ";
854  }
855  }
856  else if (isLocallyIndexed()) {
857  ArrayView<const typename MatrixType::local_ordinal_type> rowinds;
858  ArrayView<const typename MatrixType::scalar_type> rowvals;
859  getLocalRowView (r, rowinds, rowvals);
860  for (size_t j=0; j < nE; ++j) {
861  out << " (" << getColMap()->getGlobalElement(rowinds[j])
862  << ", " << rowvals[j]
863  << ") ";
864  }
865  } // globally or locally indexed
866  } // vl == VERB_EXTREME
867  out << endl;
868  } // for each row r on this process
869 
870  } // if (myRank == curRank)
871  comm->barrier();
872  comm->barrier();
873  comm->barrier();
874  }
875 
876  out.setOutputToRootOnly(0);
877  out << "===========\nlocal matrix\n=================" << std::endl;
878  out.setOutputToRootOnly(-1);
879  A_->describe(out,Teuchos::VERB_EXTREME);
880  out.setOutputToRootOnly(0);
881  out << "===========\nend of local matrix\n=================" << std::endl;
882  comm->barrier();
883  out.setOutputToRootOnly(0);
884  out << "=================\nghost matrix\n=================" << std::endl;
885  out.setOutputToRootOnly(-1);
886  ExtMatrix_->describe(out,Teuchos::VERB_EXTREME);
887  out.setOutputToRootOnly(0);
888  out << "===========\nend of ghost matrix\n=================" << std::endl;
889  }
890  }
891 }
892 
893 template<class MatrixType>
894 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
895 OverlappingRowMatrix<MatrixType>::getUnderlyingMatrix() const
896 {
897  return A_;
898 }
899 
900 
901 } // namespace Ifpack2
902 
903 #define IFPACK2_OVERLAPPINGROWMATRIX_INSTANT(S,LO,GO,N) \
904  template class Ifpack2::OverlappingRowMatrix< Tpetra::RowMatrix<S, LO, GO, N> >;
905 
906 #endif // IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
virtual Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getDomainMap() const
The Map that describes the domain of this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:267
virtual global_size_t getGlobalNumDiags() const
The global number of diagonal entries.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:377
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:283
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:291
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:489
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:258
virtual bool hasTransposeApply() const
Whether this operator&#39;s apply() method can apply the adjoint (transpose).
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:675
virtual mag_type getFrobeniusNorm() const
Returns the Frobenius norm of the matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:554
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:350
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:233
virtual bool isUpperTriangular() const
Whether this matrix is upper triangular.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:419
virtual bool isLocallyIndexed() const
Whether this matrix is locally indexed.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:426
virtual global_size_t getGlobalNumCols() const
The global number of columns in this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:305
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:537
virtual bool isFillComplete() const
true if fillComplete() has been called, else false.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:440
virtual bool isGloballyIndexed() const
Whether this matrix is globally indexed.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:433
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:249
virtual bool hasColMap() const
Whether this matrix has a column Map.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:405
virtual Teuchos::RCP< node_type > getNode() const
The matrix&#39;s Node instance.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:241
virtual bool isLowerTriangular() const
Whether this matrix is lower triangular.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:412
virtual size_t getNodeNumEntries() const
The number of entries in this matrix owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:341
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:510
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:470
virtual bool supportsRowViews() const
true if row views are supported, else false.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:682
Sparse graph (Tpetra::RowGraph subclass) with ghost rows.
Definition: Ifpack2_Details_OverlappingRowGraph_decl.hpp:65
virtual size_t getNodeNumCols() const
The number of columns owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:319
virtual size_t getGlobalMaxNumRowEntries() const
The maximum number of entries in any row on any process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:391
~OverlappingRowMatrix()
Destructor.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:228
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:546
virtual size_t getNodeNumDiags() const
The number of diagonal entries owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:384
Sparse matrix (Tpetra::RowMatrix subclass) with ghost rows.
Definition: Ifpack2_OverlappingRowMatrix_decl.hpp:59
virtual size_t getNodeMaxNumRowEntries() const
The maximum number of entries in any row on the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:398
OverlappingRowMatrix(const Teuchos::RCP< const row_matrix_type > &A, const int overlapLevel)
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:55
virtual size_t getNodeNumRows() const
The number of rows owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:312
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:563
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
virtual global_size_t getGlobalNumEntries() const
The global number of entries in this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:334
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:449
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:364
virtual global_ordinal_type getIndexBase() const
The index base for global indices for this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:327
virtual global_size_t getGlobalNumRows() const
The global number of rows in this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:298
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:528