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  return A_->getDomainMap();
270 }
271 
272 
273 template<class MatrixType>
274 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
276 {
277  return A_->getRangeMap ();
278 }
279 
280 
281 template<class MatrixType>
282 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
284 {
285  return graph_;
286 }
287 
288 
289 template<class MatrixType>
291 {
292  return NumGlobalRows_;
293 }
294 
295 
296 template<class MatrixType>
298 {
299  return NumGlobalRows_;
300 }
301 
302 
303 template<class MatrixType>
305 {
306  return A_->getNodeNumRows () + ExtMatrix_->getNodeNumRows ();
307 }
308 
309 
310 template<class MatrixType>
312 {
313  return this->getNodeNumRows ();
314 }
315 
316 
317 template<class MatrixType>
318 typename MatrixType::global_ordinal_type
320 {
321  return A_->getIndexBase();
322 }
323 
324 
325 template<class MatrixType>
327 {
328  return NumGlobalNonzeros_;
329 }
330 
331 
332 template<class MatrixType>
334 {
335  return A_->getNodeNumEntries () + ExtMatrix_->getNodeNumEntries ();
336 }
337 
338 
339 template<class MatrixType>
340 size_t
342 getNumEntriesInGlobalRow (global_ordinal_type globalRow) const
343 {
344  const local_ordinal_type localRow = RowMap_->getLocalElement (globalRow);
345  if (localRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
346  return Teuchos::OrdinalTraits<size_t>::invalid();
347  } else {
348  return getNumEntriesInLocalRow (localRow);
349  }
350 }
351 
352 
353 template<class MatrixType>
354 size_t
356 getNumEntriesInLocalRow (local_ordinal_type localRow) const
357 {
358  using Teuchos::as;
359  const size_t numMyRowsA = A_->getNodeNumRows ();
360  if (as<size_t> (localRow) < numMyRowsA) {
361  return A_->getNumEntriesInLocalRow (localRow);
362  } else {
363  return ExtMatrix_->getNumEntriesInLocalRow (as<local_ordinal_type> (localRow - numMyRowsA));
364  }
365 }
366 
367 
368 template<class MatrixType>
370 {
371  throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getGlobalNumDiags() not supported.");
372 }
373 
374 
375 template<class MatrixType>
377 {
378  return A_->getNodeNumDiags();
379 }
380 
381 
382 template<class MatrixType>
384 {
385  throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getGlobalMaxNumRowEntries() not supported.");
386 }
387 
388 
389 template<class MatrixType>
391 {
392  return MaxNumEntries_;
393 }
394 
395 
396 template<class MatrixType>
398 {
399  return true;
400 }
401 
402 
403 template<class MatrixType>
405 {
406  return A_->isLowerTriangular();
407 }
408 
409 
410 template<class MatrixType>
412 {
413  return A_->isUpperTriangular();
414 }
415 
416 
417 template<class MatrixType>
419 {
420  return true;
421 }
422 
423 
424 template<class MatrixType>
426 {
427  return false;
428 }
429 
430 
431 template<class MatrixType>
433 {
434  return true;
435 }
436 
437 
438 template<class MatrixType>
439 void
441 getGlobalRowCopy (global_ordinal_type GlobalRow,
442  const Teuchos::ArrayView<global_ordinal_type> &Indices,
443  const Teuchos::ArrayView<scalar_type>& Values,
444  size_t& NumEntries) const
445 {
446  const local_ordinal_type LocalRow = RowMap_->getLocalElement (GlobalRow);
447  if (LocalRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
448  NumEntries = Teuchos::OrdinalTraits<size_t>::invalid ();
449  } else {
450  if (Teuchos::as<size_t> (LocalRow) < A_->getNodeNumRows ()) {
451  A_->getGlobalRowCopy (GlobalRow, Indices, Values, NumEntries);
452  } else {
453  ExtMatrix_->getGlobalRowCopy (GlobalRow, Indices, Values, NumEntries);
454  }
455  }
456 }
457 
458 
459 template<class MatrixType>
460 void
462 getLocalRowCopy (local_ordinal_type LocalRow,
463  const Teuchos::ArrayView<local_ordinal_type> &Indices,
464  const Teuchos::ArrayView<scalar_type> &Values,
465  size_t &NumEntries) const
466 {
467  using Teuchos::as;
468  const size_t numMyRowsA = A_->getNodeNumRows ();
469  if (as<size_t> (LocalRow) < numMyRowsA) {
470  A_->getLocalRowCopy (LocalRow, Indices, Values, NumEntries);
471  } else {
472  ExtMatrix_->getLocalRowCopy (LocalRow - as<local_ordinal_type> (numMyRowsA),
473  Indices, Values, NumEntries);
474  }
475 }
476 
477 
478 template<class MatrixType>
479 void
481 getGlobalRowView (global_ordinal_type GlobalRow,
482  Teuchos::ArrayView<const global_ordinal_type>& indices,
483  Teuchos::ArrayView<const scalar_type>& values) const
484 {
485  const local_ordinal_type LocalRow = RowMap_->getLocalElement (GlobalRow);
486  if (LocalRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid()) {
487  indices = Teuchos::null;
488  values = Teuchos::null;
489  } else {
490  if (Teuchos::as<size_t> (LocalRow) < A_->getNodeNumRows ()) {
491  A_->getGlobalRowView (GlobalRow, indices, values);
492  } else {
493  ExtMatrix_->getGlobalRowView (GlobalRow, indices, values);
494  }
495  }
496 }
497 
498 
499 template<class MatrixType>
500 void
502 getLocalRowView (local_ordinal_type LocalRow,
503  Teuchos::ArrayView<const local_ordinal_type>& indices,
504  Teuchos::ArrayView<const scalar_type>& values) const
505 {
506  using Teuchos::as;
507  const size_t numMyRowsA = A_->getNodeNumRows ();
508  if (as<size_t> (LocalRow) < numMyRowsA) {
509  A_->getLocalRowView (LocalRow, indices, values);
510  } else {
511  ExtMatrix_->getLocalRowView (LocalRow - as<local_ordinal_type> (numMyRowsA),
512  indices, values);
513  }
514 }
515 
516 
517 template<class MatrixType>
518 void
520 getLocalDiagCopy (Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& diag) const
521 {
522  throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getLocalDiagCopy not supported.");
523 }
524 
525 
526 template<class MatrixType>
527 void
529 leftScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
530 {
531  throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
532 }
533 
534 
535 template<class MatrixType>
536 void
538 rightScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
539 {
540  throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
541 }
542 
543 
544 template<class MatrixType>
545 typename OverlappingRowMatrix<MatrixType>::mag_type
547 {
548  throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support getFrobeniusNorm.");
549 }
550 
551 
552 template<class MatrixType>
553 void
555 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
556  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
557  Teuchos::ETransp mode,
558  scalar_type alpha,
559  scalar_type beta) const
560 {
561  using Teuchos::ArrayRCP;
562  using Teuchos::as;
563  typedef scalar_type RangeScalar;
564  typedef scalar_type DomainScalar;
565  typedef Teuchos::ScalarTraits<RangeScalar> STRS;
566 
567  TEUCHOS_TEST_FOR_EXCEPTION(
568  alpha != Teuchos::ScalarTraits<scalar_type>::one () ||
569  beta != Teuchos::ScalarTraits<scalar_type>::zero (), std::logic_error,
570  "Ifpack2::ReorderFilter::apply is only implemented for alpha = 1 and "
571  "beta = 0. You set alpha = " << alpha << " and beta = " << beta << ".");
572  TEUCHOS_TEST_FOR_EXCEPTION(
573  X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
574  "Ifpack2::OverlappingRowMatrix::apply: The input X and the output Y must "
575  "have the same number of columns. X.getNumVectors() = "
576  << X.getNumVectors() << " != Y.getNumVectors() = " << Y.getNumVectors()
577  << ".");
578 
579  // FIXME (mfh 13 July 2013) This would be a good candidate for a
580  // local parallel operator implementation. That would obviate the
581  // need for getting views of the data and make the code below a lot
582  // simpler.
583 
584  const RangeScalar zero = STRS::zero ();
585  ArrayRCP<ArrayRCP<const DomainScalar> > x_ptr = X.get2dView();
586  ArrayRCP<ArrayRCP<RangeScalar> > y_ptr = Y.get2dViewNonConst();
587  Y.putScalar(zero);
588  size_t NumVectors = Y.getNumVectors();
589 
590  const size_t numMyRowsA = A_->getNodeNumRows ();
591  for (size_t i = 0; i < numMyRowsA; ++i) {
592  size_t Nnz;
593  // Use this class's getrow to make the below code simpler
594  A_->getLocalRowCopy (i, Indices_ (),Values_ (), Nnz);
595  if (mode == Teuchos::NO_TRANS) {
596  for (size_t j = 0; j < Nnz; ++j)
597  for (size_t k = 0; k < NumVectors; ++k)
598  y_ptr[k][i] += as<RangeScalar> (Values_[j]) *
599  as<RangeScalar> (x_ptr[k][Indices_[j]]);
600  }
601  else if (mode == Teuchos::TRANS){
602  for (size_t j = 0; j < Nnz; ++j)
603  for (size_t k = 0; k < NumVectors; ++k)
604  y_ptr[k][Indices_[j]] += as<RangeScalar> (Values_[j]) *
605  as<RangeScalar> (x_ptr[k][i]);
606  }
607  else { // mode == Teuchos::CONJ_TRANS
608  for (size_t j = 0; j < Nnz; ++j)
609  for (size_t k = 0; k < NumVectors; ++k)
610  y_ptr[k][Indices_[j]] +=
611  STRS::conjugate (as<RangeScalar> (Values_[j])) *
612  as<RangeScalar> (x_ptr[k][i]);
613  }
614  }
615 
616  const size_t numMyRowsB = ExtMatrix_->getNodeNumRows ();
617  for (size_t i = 0 ; i < numMyRowsB ; ++i) {
618  size_t Nnz;
619  // Use this class's getrow to make the below code simpler
620  ExtMatrix_->getLocalRowCopy (i, Indices_ (), Values_ (), Nnz);
621  if (mode == Teuchos::NO_TRANS) {
622  for (size_t j = 0; j < Nnz; ++j)
623  for (size_t k = 0; k < NumVectors; ++k)
624  y_ptr[k][numMyRowsA+i] += as<RangeScalar> (Values_[j]) *
625  as<RangeScalar> (x_ptr[k][Indices_[j]]);
626  }
627  else if (mode == Teuchos::TRANS) {
628  for (size_t j = 0; j < Nnz; ++j)
629  for (size_t k = 0; k < NumVectors; ++k)
630  y_ptr[k][numMyRowsA+Indices_[j]] += as<RangeScalar> (Values_[j]) *
631  as<RangeScalar> (x_ptr[k][i]);
632  }
633  else { // mode == Teuchos::CONJ_TRANS
634  for (size_t j = 0; j < Nnz; ++j)
635  for (size_t k = 0; k < NumVectors; ++k)
636  y_ptr[k][numMyRowsA+Indices_[j]] +=
637  STRS::conjugate (as<RangeScalar> (Values_[j])) *
638  as<RangeScalar> (x_ptr[k][i]);
639  }
640  }
641 }
642 
643 
644 template<class MatrixType>
645 void
647 importMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
648  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
649  Tpetra::CombineMode CM)
650 {
651  OvX.doImport (X, *Importer_, CM);
652 }
653 
654 
655 template<class MatrixType>
656 void
658 exportMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
659  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
660  Tpetra::CombineMode CM)
661 {
662  X.doExport (OvX, *Importer_, CM);
663 }
664 
665 
666 template<class MatrixType>
668 {
669  return true;
670 }
671 
672 
673 template<class MatrixType>
675 {
676  return false;
677 }
678 
679 template<class MatrixType>
681 {
682  std::ostringstream oss;
683  if (isFillComplete()) {
684  oss << "{ isFillComplete: true"
685  << ", global rows: " << getGlobalNumRows()
686  << ", global columns: " << getGlobalNumCols()
687  << ", global entries: " << getGlobalNumEntries()
688  << " }";
689  }
690  else {
691  oss << "{ isFillComplete: false"
692  << ", global rows: " << getGlobalNumRows()
693  << " }";
694  }
695  return oss.str();
696 }
697 
698 template<class MatrixType>
699 void OverlappingRowMatrix<MatrixType>::describe(Teuchos::FancyOStream &out,
700  const Teuchos::EVerbosityLevel verbLevel) const
701 {
702  using std::endl;
703  using std::setw;
704  using Teuchos::as;
705  using Teuchos::VERB_DEFAULT;
706  using Teuchos::VERB_NONE;
707  using Teuchos::VERB_LOW;
708  using Teuchos::VERB_MEDIUM;
709  using Teuchos::VERB_HIGH;
710  using Teuchos::VERB_EXTREME;
711  using Teuchos::RCP;
712  using Teuchos::null;
713  using Teuchos::ArrayView;
714 
715  Teuchos::EVerbosityLevel vl = verbLevel;
716  if (vl == VERB_DEFAULT) {
717  vl = VERB_LOW;
718  }
719  RCP<const Teuchos::Comm<int> > comm = this->getComm();
720  const int myRank = comm->getRank();
721  const int numProcs = comm->getSize();
722  size_t width = 1;
723  for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
724  ++width;
725  }
726  width = std::max<size_t> (width, as<size_t> (11)) + 2;
727  Teuchos::OSTab tab(out);
728  // none: print nothing
729  // low: print O(1) info from node 0
730  // medium: print O(P) info, num entries per process
731  // high: print O(N) info, num entries per row
732  // extreme: print O(NNZ) info: print indices and values
733  //
734  // for medium and higher, print constituent objects at specified verbLevel
735  if (vl != VERB_NONE) {
736  if (myRank == 0) {
737  out << this->description() << std::endl;
738  }
739  // O(1) globals, minus what was already printed by description()
740  //if (isFillComplete() && myRank == 0) {
741  // out << "Global number of diagonal entries: " << getGlobalNumDiags() << std::endl;
742  // out << "Global max number of entries in a row: " << getGlobalMaxNumRowEntries() << std::endl;
743  //}
744  // constituent objects
745  if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
746  if (myRank == 0) {
747  out << endl << "Row map:" << endl;
748  }
749  getRowMap()->describe(out,vl);
750  //
751  if (getColMap() != null) {
752  if (getColMap() == getRowMap()) {
753  if (myRank == 0) {
754  out << endl << "Column map is row map.";
755  }
756  }
757  else {
758  if (myRank == 0) {
759  out << endl << "Column map:" << endl;
760  }
761  getColMap()->describe(out,vl);
762  }
763  }
764  if (getDomainMap() != null) {
765  if (getDomainMap() == getRowMap()) {
766  if (myRank == 0) {
767  out << endl << "Domain map is row map.";
768  }
769  }
770  else if (getDomainMap() == getColMap()) {
771  if (myRank == 0) {
772  out << endl << "Domain map is column map.";
773  }
774  }
775  else {
776  if (myRank == 0) {
777  out << endl << "Domain map:" << endl;
778  }
779  getDomainMap()->describe(out,vl);
780  }
781  }
782  if (getRangeMap() != null) {
783  if (getRangeMap() == getDomainMap()) {
784  if (myRank == 0) {
785  out << endl << "Range map is domain map." << endl;
786  }
787  }
788  else if (getRangeMap() == getRowMap()) {
789  if (myRank == 0) {
790  out << endl << "Range map is row map." << endl;
791  }
792  }
793  else {
794  if (myRank == 0) {
795  out << endl << "Range map: " << endl;
796  }
797  getRangeMap()->describe(out,vl);
798  }
799  }
800  if (myRank == 0) {
801  out << endl;
802  }
803  }
804  // O(P) data
805  if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
806  for (int curRank = 0; curRank < numProcs; ++curRank) {
807  if (myRank == curRank) {
808  out << "Process rank: " << curRank << std::endl;
809  out << " Number of entries: " << getNodeNumEntries() << std::endl;
810  if (isFillComplete()) {
811  out << " Number of diagonal entries: " << getNodeNumDiags() << std::endl;
812  }
813  out << " Max number of entries per row: " << getNodeMaxNumRowEntries() << std::endl;
814  }
815  comm->barrier();
816  comm->barrier();
817  comm->barrier();
818  }
819  }
820  // O(N) and O(NNZ) data
821  if (vl == VERB_HIGH || vl == VERB_EXTREME) {
822  for (int curRank = 0; curRank < numProcs; ++curRank) {
823  if (myRank == curRank) {
824  out << std::setw(width) << "Proc Rank"
825  << std::setw(width) << "Global Row"
826  << std::setw(width) << "Num Entries";
827  if (vl == VERB_EXTREME) {
828  out << std::setw(width) << "(Index,Value)";
829  }
830  out << endl;
831  for (size_t r = 0; r < getNodeNumRows (); ++r) {
832  const size_t nE = getNumEntriesInLocalRow(r);
833  typename MatrixType::global_ordinal_type gid = getRowMap()->getGlobalElement(r);
834  out << std::setw(width) << myRank
835  << std::setw(width) << gid
836  << std::setw(width) << nE;
837  if (vl == VERB_EXTREME) {
838  if (isGloballyIndexed()) {
839  ArrayView<const typename MatrixType::global_ordinal_type> rowinds;
840  ArrayView<const typename MatrixType::scalar_type> rowvals;
841  getGlobalRowView (gid, rowinds, rowvals);
842  for (size_t j = 0; j < nE; ++j) {
843  out << " (" << rowinds[j]
844  << ", " << rowvals[j]
845  << ") ";
846  }
847  }
848  else if (isLocallyIndexed()) {
849  ArrayView<const typename MatrixType::local_ordinal_type> rowinds;
850  ArrayView<const typename MatrixType::scalar_type> rowvals;
851  getLocalRowView (r, rowinds, rowvals);
852  for (size_t j=0; j < nE; ++j) {
853  out << " (" << getColMap()->getGlobalElement(rowinds[j])
854  << ", " << rowvals[j]
855  << ") ";
856  }
857  } // globally or locally indexed
858  } // vl == VERB_EXTREME
859  out << endl;
860  } // for each row r on this process
861 
862  } // if (myRank == curRank)
863  comm->barrier();
864  comm->barrier();
865  comm->barrier();
866  }
867 
868  out.setOutputToRootOnly(0);
869  out << "===========\nlocal matrix\n=================" << std::endl;
870  out.setOutputToRootOnly(-1);
871  A_->describe(out,Teuchos::VERB_EXTREME);
872  out.setOutputToRootOnly(0);
873  out << "===========\nend of local matrix\n=================" << std::endl;
874  comm->barrier();
875  out.setOutputToRootOnly(0);
876  out << "=================\nghost matrix\n=================" << std::endl;
877  out.setOutputToRootOnly(-1);
878  ExtMatrix_->describe(out,Teuchos::VERB_EXTREME);
879  out.setOutputToRootOnly(0);
880  out << "===========\nend of ghost matrix\n=================" << std::endl;
881  }
882  }
883 }
884 
885 template<class MatrixType>
886 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
888 {
889  return A_;
890 }
891 
892 
893 } // namespace Ifpack2
894 
895 #define IFPACK2_OVERLAPPINGROWMATRIX_INSTANT(S,LO,GO,N) \
896  template class Ifpack2::OverlappingRowMatrix< Tpetra::RowMatrix<S, LO, GO, N> >;
897 
898 #endif // IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
virtual bool isGloballyIndexed() const
Whether this matrix is globally indexed.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:425
virtual mag_type getFrobeniusNorm() const
Returns the Frobenius norm of the matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:546
virtual size_t getNodeNumRows() const
The number of rows owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:304
virtual bool isFillComplete() const
true if fillComplete() has been called, else false.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:432
virtual bool isUpperTriangular() const
Whether this matrix is upper triangular.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:411
virtual Teuchos::RCP< node_type > getNode() const
The matrix&#39;s Node instance.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:241
virtual size_t getNodeNumDiags() const
The number of diagonal entries owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:376
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:529
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:233
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 isLowerTriangular() const
Whether this matrix is lower triangular.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:404
virtual global_size_t getGlobalNumRows() const
The global number of rows in this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:290
virtual global_size_t getGlobalNumCols() const
The global number of columns in this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:297
virtual bool isLocallyIndexed() const
Whether this matrix is locally indexed.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:418
virtual global_size_t getGlobalNumDiags() const
The global number of diagonal entries.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:369
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:275
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
Sparse graph (Tpetra::RowGraph subclass) with ghost rows.
Definition: Ifpack2_Details_OverlappingRowGraph_decl.hpp:65
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:342
virtual size_t getNodeNumEntries() const
The number of entries in this matrix owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:333
virtual global_size_t getGlobalNumEntries() const
The global number of entries in this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:326
~OverlappingRowMatrix()
Destructor.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:228
virtual bool hasTransposeApply() const
Whether this operator&#39;s apply() method can apply the adjoint (transpose).
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:667
virtual global_ordinal_type getIndexBase() const
The index base for global indices for this matrix.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:319
virtual bool hasColMap() const
Whether this matrix has a column Map.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:397
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:356
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:538
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:283
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:481
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:520
Sparse matrix (Tpetra::RowMatrix subclass) with ghost rows.
Definition: Ifpack2_OverlappingRowMatrix_decl.hpp:59
OverlappingRowMatrix(const Teuchos::RCP< const row_matrix_type > &A, const int overlapLevel)
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:55
virtual size_t getNodeMaxNumRowEntries() const
The maximum number of entries in any row on the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:390
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:555
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:462
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:441
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
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 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:502
virtual bool supportsRowViews() const
true if row views are supported, else false.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:674
virtual size_t getGlobalMaxNumRowEntries() const
The maximum number of entries in any row on any process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:383
virtual size_t getNodeNumCols() const
The number of columns owned by the calling process.
Definition: Ifpack2_OverlappingRowMatrix_def.hpp:311