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