Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Tpetra_RowMatrix_def.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Tpetra: Templated Linear Algebra Services Package
5 // Copyright (2008) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
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 #ifndef TPETRA_ROWMATRIX_DEF_HPP
43 #define TPETRA_ROWMATRIX_DEF_HPP
44 
45 #include "Tpetra_CrsMatrix.hpp"
46 #include "Tpetra_Map.hpp"
47 #include "Tpetra_RowGraph.hpp"
48 
49 namespace Tpetra {
50 
51  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
53 
54  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
55  Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
57  add (const Scalar& alpha,
59  const Scalar& beta,
60  const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& domainMap,
61  const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& rangeMap,
62  const Teuchos::RCP<Teuchos::ParameterList>& params) const
63  {
64  using Teuchos::Array;
65  using Teuchos::ArrayView;
66  using Teuchos::ParameterList;
67  using Teuchos::RCP;
68  using Teuchos::rcp;
69  using Teuchos::rcp_implicit_cast;
70  using Teuchos::sublist;
71  typedef LocalOrdinal LO;
72  typedef GlobalOrdinal GO;
73  typedef Teuchos::ScalarTraits<Scalar> STS;
74  typedef Map<LO, GO, Node> map_type;
77 
78  const this_type& B = *this; // a convenient abbreviation
79 
80  // If the user didn't supply a domain or range Map, then try to
81  // get one from B first (if it has them), then from A (if it has
82  // them). If we don't have any domain or range Maps, scold the
83  // user.
84  RCP<const map_type> A_domainMap = A.getDomainMap ();
85  RCP<const map_type> A_rangeMap = A.getRangeMap ();
86  RCP<const map_type> B_domainMap = B.getDomainMap ();
87  RCP<const map_type> B_rangeMap = B.getRangeMap ();
88 
89  RCP<const map_type> theDomainMap = domainMap;
90  RCP<const map_type> theRangeMap = rangeMap;
91 
92  if (domainMap.is_null ()) {
93  if (B_domainMap.is_null ()) {
94  TEUCHOS_TEST_FOR_EXCEPTION(
95  A_domainMap.is_null (), std::invalid_argument,
96  "Tpetra::RowMatrix::add: If neither A nor B have a domain Map, "
97  "then you must supply a nonnull domain Map to this method.");
98  theDomainMap = A_domainMap;
99  } else {
100  theDomainMap = B_domainMap;
101  }
102  }
103  if (rangeMap.is_null ()) {
104  if (B_rangeMap.is_null ()) {
105  TEUCHOS_TEST_FOR_EXCEPTION(
106  A_rangeMap.is_null (), std::invalid_argument,
107  "Tpetra::RowMatrix::add: If neither A nor B have a range Map, "
108  "then you must supply a nonnull range Map to this method.");
109  theRangeMap = A_rangeMap;
110  } else {
111  theRangeMap = B_rangeMap;
112  }
113  }
114 
115 #ifdef HAVE_TPETRA_DEBUG
116  // In a debug build, check that A and B have matching domain and
117  // range Maps, if they have domain and range Maps at all. (If
118  // they aren't fill complete, then they may not yet have them.)
119  if (! A_domainMap.is_null () && ! A_rangeMap.is_null ()) {
120  if (! B_domainMap.is_null () && ! B_rangeMap.is_null ()) {
121  TEUCHOS_TEST_FOR_EXCEPTION(
122  ! B_domainMap->isSameAs (*A_domainMap), std::invalid_argument,
123  "Tpetra::RowMatrix::add: The input RowMatrix A must have a domain Map "
124  "which is the same as (isSameAs) this RowMatrix's domain Map.");
125  TEUCHOS_TEST_FOR_EXCEPTION(
126  ! B_rangeMap->isSameAs (*A_rangeMap), std::invalid_argument,
127  "Tpetra::RowMatrix::add: The input RowMatrix A must have a range Map "
128  "which is the same as (isSameAs) this RowMatrix's range Map.");
129  TEUCHOS_TEST_FOR_EXCEPTION(
130  ! domainMap.is_null () && ! domainMap->isSameAs (*B_domainMap),
131  std::invalid_argument,
132  "Tpetra::RowMatrix::add: The input domain Map must be the same as "
133  "(isSameAs) this RowMatrix's domain Map.");
134  TEUCHOS_TEST_FOR_EXCEPTION(
135  ! rangeMap.is_null () && ! rangeMap->isSameAs (*B_rangeMap),
136  std::invalid_argument,
137  "Tpetra::RowMatrix::add: The input range Map must be the same as "
138  "(isSameAs) this RowMatrix's range Map.");
139  }
140  }
141  else if (! B_domainMap.is_null () && ! B_rangeMap.is_null ()) {
142  TEUCHOS_TEST_FOR_EXCEPTION(
143  ! domainMap.is_null () && ! domainMap->isSameAs (*B_domainMap),
144  std::invalid_argument,
145  "Tpetra::RowMatrix::add: The input domain Map must be the same as "
146  "(isSameAs) this RowMatrix's domain Map.");
147  TEUCHOS_TEST_FOR_EXCEPTION(
148  ! rangeMap.is_null () && ! rangeMap->isSameAs (*B_rangeMap),
149  std::invalid_argument,
150  "Tpetra::RowMatrix::add: The input range Map must be the same as "
151  "(isSameAs) this RowMatrix's range Map.");
152  }
153  else {
154  TEUCHOS_TEST_FOR_EXCEPTION(
155  domainMap.is_null () || rangeMap.is_null (), std::invalid_argument,
156  "Tpetra::RowMatrix::add: If neither A nor B have a domain and range "
157  "Map, then you must supply a nonnull domain and range Map to this "
158  "method.");
159  }
160 #endif // HAVE_TPETRA_DEBUG
161 
162  // What parameters do we pass to C's constructor? Do we call
163  // fillComplete on C after filling it? And if so, what parameters
164  // do we pass to C's fillComplete call?
165  bool callFillComplete = true;
166  RCP<ParameterList> constructorSublist;
167  RCP<ParameterList> fillCompleteSublist;
168  if (! params.is_null ()) {
169  callFillComplete = params->get ("Call fillComplete", callFillComplete);
170  constructorSublist = sublist (params, "Constructor parameters");
171  fillCompleteSublist = sublist (params, "fillComplete parameters");
172  }
173 
174  RCP<const map_type> A_rowMap = A.getRowMap ();
175  RCP<const map_type> B_rowMap = B.getRowMap ();
176  RCP<const map_type> C_rowMap = B_rowMap; // see discussion in documentation
177  RCP<crs_matrix_type> C; // The result matrix.
178 
179  // If A and B's row Maps are the same, we can compute an upper
180  // bound on the number of entries in each row of C, before
181  // actually computing the sum. A reasonable upper bound is the
182  // sum of the two entry counts in each row. If we choose this as
183  // the actual per-row upper bound, we can use static profile.
184  if (A_rowMap->isSameAs (*B_rowMap)) {
185  const LO localNumRows = static_cast<LO> (A_rowMap->getNodeNumElements ());
186  Array<size_t> C_maxNumEntriesPerRow (localNumRows, 0);
187 
188  // Get the number of entries in each row of A.
189  if (alpha != STS::zero ()) {
190  for (LO localRow = 0; localRow < localNumRows; ++localRow) {
191  const size_t A_numEntries = A.getNumEntriesInLocalRow (localRow);
192  C_maxNumEntriesPerRow[localRow] += A_numEntries;
193  }
194  }
195  // Get the number of entries in each row of B.
196  if (beta != STS::zero ()) {
197  for (LO localRow = 0; localRow < localNumRows; ++localRow) {
198  const size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
199  C_maxNumEntriesPerRow[localRow] += B_numEntries;
200  }
201  }
202  // Construct the result matrix C.
203  if (constructorSublist.is_null ()) {
204  C = rcp (new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow (),
205  StaticProfile));
206  } else {
207  C = rcp (new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow (),
208  StaticProfile, constructorSublist));
209  }
210  // Since A and B have the same row Maps, we could add them
211  // together all at once and merge values before we call
212  // insertGlobalValues. However, we don't really need to, since
213  // we've already allocated enough space in each row of C for C
214  // to do the merge itself.
215  }
216  else { // the row Maps of A and B are not the same
217  // Construct the result matrix C.
218  // true: !A_rowMap->isSameAs (*B_rowMap)
219  TEUCHOS_TEST_FOR_EXCEPTION(true,
220  std::invalid_argument,
221  "Tpetra::RowMatrix::add: The row maps must be the same for statically "
222  "allocated matrices in order to be sure that there is sufficient space "
223  "to do the addition");
224  }
225 
226 #ifdef HAVE_TPETRA_DEBUG
227  TEUCHOS_TEST_FOR_EXCEPTION(C.is_null (), std::logic_error,
228  "Tpetra::RowMatrix::add: C should not be null at this point. "
229  "Please report this bug to the Tpetra developers.");
230 #endif // HAVE_TPETRA_DEBUG
231  //
232  // Compute C = alpha*A + beta*B.
233  //
234  Array<GO> ind;
235  Array<Scalar> val;
236 
237  if (alpha != STS::zero ()) {
238  const LO A_localNumRows = static_cast<LO> (A_rowMap->getNodeNumElements ());
239  for (LO localRow = 0; localRow < A_localNumRows; ++localRow) {
240  size_t A_numEntries = A.getNumEntriesInLocalRow (localRow);
241  const GO globalRow = A_rowMap->getGlobalElement (localRow);
242  if (A_numEntries > static_cast<size_t> (ind.size ())) {
243  ind.resize (A_numEntries);
244  val.resize (A_numEntries);
245  }
246  ArrayView<GO> indView = ind (0, A_numEntries);
247  ArrayView<Scalar> valView = val (0, A_numEntries);
248  A.getGlobalRowCopy (globalRow, indView, valView, A_numEntries);
249 
250  if (alpha != STS::one ()) {
251  for (size_t k = 0; k < A_numEntries; ++k) {
252  valView[k] *= alpha;
253  }
254  }
255  C->insertGlobalValues (globalRow, indView, valView);
256  }
257  }
258 
259  if (beta != STS::zero ()) {
260  const LO B_localNumRows = static_cast<LO> (B_rowMap->getNodeNumElements ());
261  for (LO localRow = 0; localRow < B_localNumRows; ++localRow) {
262  size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
263  const GO globalRow = B_rowMap->getGlobalElement (localRow);
264  if (B_numEntries > static_cast<size_t> (ind.size ())) {
265  ind.resize (B_numEntries);
266  val.resize (B_numEntries);
267  }
268  ArrayView<GO> indView = ind (0, B_numEntries);
269  ArrayView<Scalar> valView = val (0, B_numEntries);
270  B.getGlobalRowCopy (globalRow, indView, valView, B_numEntries);
271 
272  if (beta != STS::one ()) {
273  for (size_t k = 0; k < B_numEntries; ++k) {
274  valView[k] *= beta;
275  }
276  }
277  C->insertGlobalValues (globalRow, indView, valView);
278  }
279  }
280 
281  if (callFillComplete) {
282  if (fillCompleteSublist.is_null ()) {
283  C->fillComplete (theDomainMap, theRangeMap);
284  } else {
285  C->fillComplete (theDomainMap, theRangeMap, fillCompleteSublist);
286  }
287  }
288 
289  return rcp_implicit_cast<this_type> (C);
290  }
291 
292 
293  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
294  void
296  pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
297  Teuchos::Array<char>& exports,
298  const Teuchos::ArrayView<size_t>& numPacketsPerLID,
299  size_t& constantNumPackets,
300  Distributor &distor) const
301  {
302 #ifdef HAVE_TPETRA_DEBUG
303  const char tfecfFuncName[] = "pack: ";
304  {
305  using Teuchos::reduceAll;
306  std::ostringstream msg;
307  int lclBad = 0;
308  try {
309  this->packImpl (exportLIDs, exports, numPacketsPerLID,
310  constantNumPackets, distor);
311  } catch (std::exception& e) {
312  lclBad = 1;
313  msg << e.what ();
314  }
315  int gblBad = 0;
316  const Teuchos::Comm<int>& comm = * (this->getComm ());
317  reduceAll<int, int> (comm, Teuchos::REDUCE_MAX,
318  lclBad, Teuchos::outArg (gblBad));
319  if (gblBad != 0) {
320  const int myRank = comm.getRank ();
321  const int numProcs = comm.getSize ();
322  for (int r = 0; r < numProcs; ++r) {
323  if (r == myRank && lclBad != 0) {
324  std::ostringstream os;
325  os << "Proc " << myRank << ": " << msg.str () << std::endl;
326  std::cerr << os.str ();
327  }
328  comm.barrier ();
329  comm.barrier ();
330  comm.barrier ();
331  }
332  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
333  true, std::logic_error, "packImpl() threw an exception on one or "
334  "more participating processes.");
335  }
336  }
337 #else
338  this->packImpl (exportLIDs, exports, numPacketsPerLID,
339  constantNumPackets, distor);
340 #endif // HAVE_TPETRA_DEBUG
341  }
342 
343  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
344  void
346  allocatePackSpace (Teuchos::Array<char>& exports,
347  size_t& totalNumEntries,
348  const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs) const
349  {
350  typedef LocalOrdinal LO;
351  typedef GlobalOrdinal GO;
352  typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
353  //const char tfecfFuncName[] = "allocatePackSpace: ";
354  const size_type numExportLIDs = exportLIDs.size ();
355 
356  // Count the total number of entries to send.
357  totalNumEntries = 0;
358  for (size_type i = 0; i < numExportLIDs; ++i) {
359  const LO lclRow = exportLIDs[i];
360  size_t curNumEntries = this->getNumEntriesInLocalRow (lclRow);
361  // FIXME (mfh 25 Jan 2015) We should actually report invalid row
362  // indices as an error. Just consider them nonowned for now.
363  if (curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid ()) {
364  curNumEntries = 0;
365  }
366  totalNumEntries += curNumEntries;
367  }
368 
369  // FIXME (mfh 24 Feb 2013) This code is only correct if
370  // sizeof(Scalar) is a meaningful representation of the amount of
371  // data in a Scalar instance. (LO and GO are always built-in
372  // integer types.)
373  //
374  // Allocate the exports array. It does NOT need padding for
375  // alignment, since we use memcpy to write to / read from send /
376  // receive buffers.
377  const size_t allocSize =
378  static_cast<size_t> (numExportLIDs) * sizeof (LO) +
379  totalNumEntries * (sizeof (Scalar) + sizeof (GO));
380  if (static_cast<size_t> (exports.size ()) < allocSize) {
381  exports.resize (allocSize);
382  }
383  }
384 
385  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
386  bool
387  RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
388  packRow (char* const numEntOut,
389  char* const valOut,
390  char* const indOut,
391  const size_t numEnt,
392  const LocalOrdinal lclRow) const
393  {
394  using Teuchos::Array;
395  using Teuchos::ArrayView;
396  typedef LocalOrdinal LO;
397  typedef GlobalOrdinal GO;
398  typedef Tpetra::Map<LO, GO, Node> map_type;
399 
400  const LO numEntLO = static_cast<LO> (numEnt);
401  memcpy (numEntOut, &numEntLO, sizeof (LO));
402 
403  if (this->supportsRowViews ()) {
404  if (this->isLocallyIndexed ()) {
405  // If the matrix is locally indexed on the calling process, we
406  // have to use its column Map (which it _must_ have in this
407  // case) to convert to global indices.
408  ArrayView<const LO> indIn;
409  ArrayView<const Scalar> valIn;
410  this->getLocalRowView (lclRow, indIn, valIn);
411  const map_type& colMap = * (this->getColMap ());
412  // Copy column indices one at a time, so that we don't need
413  // temporary storage.
414  for (size_t k = 0; k < numEnt; ++k) {
415  const GO gblIndIn = colMap.getGlobalElement (indIn[k]);
416  memcpy (indOut + k * sizeof (GO), &gblIndIn, sizeof (GO));
417  }
418  memcpy (valOut, valIn.getRawPtr (), numEnt * sizeof (Scalar));
419  }
420  else if (this->isGloballyIndexed ()) {
421  // If the matrix is globally indexed on the calling process,
422  // then we can use the column indices directly. However, we
423  // have to get the global row index. The calling process must
424  // have a row Map, since otherwise it shouldn't be participating
425  // in packing operations.
426  ArrayView<const GO> indIn;
427  ArrayView<const Scalar> valIn;
428  const map_type& rowMap = * (this->getRowMap ());
429  const GO gblRow = rowMap.getGlobalElement (lclRow);
430  this->getGlobalRowView (gblRow, indIn, valIn);
431  memcpy (indOut, indIn.getRawPtr (), numEnt * sizeof (GO));
432  memcpy (valOut, valIn.getRawPtr (), numEnt * sizeof (Scalar));
433  }
434  else {
435  if (numEnt != 0) {
436  return false;
437  }
438  }
439  }
440  else {
441  // FIXME (mfh 25 Jan 2015) Pass in valIn and indIn as scratch
442  // space, instead of allocating them on each call.
443  if (this->isLocallyIndexed ()) {
444  Array<LO> indIn (numEnt);
445  Array<Scalar> valIn (numEnt);
446  size_t theNumEnt = 0;
447  this->getLocalRowCopy (lclRow, indIn (), valIn (), theNumEnt);
448  if (theNumEnt != numEnt) {
449  return false;
450  }
451  const map_type& colMap = * (this->getColMap ());
452  // Copy column indices one at a time, so that we don't need
453  // temporary storage.
454  for (size_t k = 0; k < numEnt; ++k) {
455  const GO gblIndIn = colMap.getGlobalElement (indIn[k]);
456  memcpy (indOut + k * sizeof (GO), &gblIndIn, sizeof (GO));
457  }
458  memcpy (valOut, valIn.getRawPtr (), numEnt * sizeof (Scalar));
459  }
460  else if (this->isGloballyIndexed ()) {
461  Array<GO> indIn (numEnt);
462  Array<Scalar> valIn (numEnt);
463  const map_type& rowMap = * (this->getRowMap ());
464  const GO gblRow = rowMap.getGlobalElement (lclRow);
465  size_t theNumEnt = 0;
466  this->getGlobalRowCopy (gblRow, indIn, valIn, theNumEnt);
467  if (theNumEnt != numEnt) {
468  return false;
469  }
470  memcpy (indOut, indIn.getRawPtr (), numEnt * sizeof (GO));
471  memcpy (valOut, valIn.getRawPtr (), numEnt * sizeof (Scalar));
472  }
473  else {
474  if (numEnt != 0) {
475  return false;
476  }
477  }
478  }
479  return true;
480  }
481 
482  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
483  void
484  RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
485  packImpl (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
486  Teuchos::Array<char>& exports,
487  const Teuchos::ArrayView<size_t>& numPacketsPerLID,
488  size_t& constantNumPackets,
489  Distributor& /* distor */) const
490  {
491  using Teuchos::Array;
492  using Teuchos::ArrayView;
493  using Teuchos::as;
494  using Teuchos::av_reinterpret_cast;
495  using Teuchos::RCP;
496  typedef LocalOrdinal LO;
497  typedef GlobalOrdinal GO;
498  typedef typename ArrayView<const LO>::size_type size_type;
499  const char tfecfFuncName[] = "packImpl: ";
500 
501  const size_type numExportLIDs = exportLIDs.size ();
502  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
503  numExportLIDs != numPacketsPerLID.size (), std::invalid_argument,
504  "exportLIDs.size() = " << numExportLIDs << " != numPacketsPerLID.size()"
505  " = " << numPacketsPerLID.size () << ".");
506 
507  // Setting this to zero tells the caller to expect a possibly
508  // different ("nonconstant") number of packets per local index
509  // (i.e., a possibly different number of entries per row).
510  constantNumPackets = 0;
511 
512  // The pack buffer 'exports' enters this method possibly
513  // unallocated. Do the first two parts of "Count, allocate, fill,
514  // compute."
515  size_t totalNumEntries = 0;
516  allocatePackSpace (exports, totalNumEntries, exportLIDs);
517  const size_t bufSize = static_cast<size_t> (exports.size ());
518 
519  // Compute the number of "packets" (in this case, bytes) per
520  // export LID (in this case, local index of the row to send), and
521  // actually pack the data.
522  //
523  // FIXME (mfh 24 Feb 2013, 25 Jan 2015) This code is only correct
524  // if sizeof(Scalar) is a meaningful representation of the amount
525  // of data in a Scalar instance. (LO and GO are always built-in
526  // integer types.)
527 
528  // Variables for error reporting in the loop.
529  size_type firstBadIndex = 0; // only valid if outOfBounds == true.
530  size_t firstBadOffset = 0; // only valid if outOfBounds == true.
531  size_t firstBadNumBytes = 0; // only valid if outOfBounds == true.
532  bool outOfBounds = false;
533  bool packErr = false;
534 
535  char* const exportsRawPtr = exports.getRawPtr ();
536  size_t offset = 0; // current index into 'exports' array.
537  for (size_type i = 0; i < numExportLIDs; ++i) {
538  const LO lclRow = exportLIDs[i];
539  const size_t numEnt = this->getNumEntriesInLocalRow (lclRow);
540 
541  // Only pad this row if it has a nonzero number of entries.
542  if (numEnt == 0) {
543  numPacketsPerLID[i] = 0;
544  }
545  else {
546  char* const numEntBeg = exportsRawPtr + offset;
547  char* const numEntEnd = numEntBeg + sizeof (LO);
548  char* const valBeg = numEntEnd;
549  char* const valEnd = valBeg + numEnt * sizeof (Scalar);
550  char* const indBeg = valEnd;
551  const size_t numBytes = sizeof (LO) +
552  numEnt * (sizeof (Scalar) + sizeof (GO));
553  if (offset > bufSize || offset + numBytes > bufSize) {
554  firstBadIndex = i;
555  firstBadOffset = offset;
556  firstBadNumBytes = numBytes;
557  outOfBounds = true;
558  break;
559  }
560  packErr = ! packRow (numEntBeg, valBeg, indBeg, numEnt, lclRow);
561  if (packErr) {
562  firstBadIndex = i;
563  firstBadOffset = offset;
564  firstBadNumBytes = numBytes;
565  break;
566  }
567  // numPacketsPerLID[i] is the number of "packets" in the
568  // current local row i. Packet=char (really "byte") so use
569  // the number of bytes of the packed data for that row.
570  numPacketsPerLID[i] = numBytes;
571  offset += numBytes;
572  }
573  }
574 
575  // The point of these exceptions is to stop computation if the
576  // above checks found a bug. If HAVE_TPETRA_DEBUG is defined,
577  // Tpetra will do extra all-reduces to check for global
578  // consistency of the error state. Otherwise, throwing an
579  // exception here might cause deadlock, but we accept that as
580  // better than just continuing.
581  TEUCHOS_TEST_FOR_EXCEPTION(
582  outOfBounds, std::logic_error, "First invalid offset into 'exports' "
583  "pack buffer at index i = " << firstBadIndex << ". exportLIDs[i]: "
584  << exportLIDs[firstBadIndex] << ", bufSize: " << bufSize << ", offset: "
585  << firstBadOffset << ", numBytes: " << firstBadNumBytes << ".");
586  TEUCHOS_TEST_FOR_EXCEPTION(
587  packErr, std::logic_error, "First error in packRow() at index i = "
588  << firstBadIndex << ". exportLIDs[i]: " << exportLIDs[firstBadIndex]
589  << ", bufSize: " << bufSize << ", offset: " << firstBadOffset
590  << ", numBytes: " << firstBadNumBytes << ".");
591  }
592 
593  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
594  LocalOrdinal
596  getLocalRowViewRaw (const LocalOrdinal lclRow,
597  LocalOrdinal& numEnt,
598  const LocalOrdinal*& lclColInds,
599  const Scalar*& vals) const
600  {
601  // This is just the default implementation. Subclasses may want
602  // to implement this method in a more efficient way, e.g., to
603  // avoid creating Teuchos::ArrayView instances.
604  Teuchos::ArrayView<const LocalOrdinal> lclColInds_av;
605  Teuchos::ArrayView<const Scalar> vals_av;
606 
607  this->getLocalRowView (lclRow, lclColInds_av, vals_av);
608  numEnt = static_cast<LocalOrdinal> (lclColInds_av.size ());
609  if (numEnt == 0) {
610  lclColInds = NULL;
611  vals = NULL;
612  }
613  else {
614  lclColInds = lclColInds_av.getRawPtr ();
615  vals = vals_av.getRawPtr ();
616  }
617 
618  return static_cast<LocalOrdinal> (0);
619  }
620 
621 } // namespace Tpetra
622 
623 //
624 // Explicit instantiation macro
625 //
626 // Must be expanded from within the Tpetra namespace!
627 //
628 
629 #define TPETRA_ROWMATRIX_INSTANT(SCALAR,LO,GO,NODE) \
630  template class RowMatrix< SCALAR , LO , GO , NODE >;
631 
632 
633 #endif // TPETRA_ROWMATRIX_DEF_HPP
634 
KOKKOS_INLINE_FUNCTION bool outOfBounds(const IntegerType x, const IntegerType exclusiveUpperBound)
Is x out of bounds? That is, is x less than zero, or greater than or equal to the given exclusive upp...
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
The current number of entries on the calling process in the specified local row.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
virtual LocalOrdinal getLocalRowViewRaw(const LocalOrdinal lclRow, LocalOrdinal &numEnt, const LocalOrdinal *&lclColInds, const Scalar *&vals) const
Get a constant, nonpersisting, locally indexed view of the given row of the matrix, using &quot;raw&quot; pointers instead of Teuchos::ArrayView.
virtual ~RowMatrix()
Destructor (virtual for memory safety of derived classes).
virtual void getGlobalRowCopy(GlobalOrdinal GlobalRow, const Teuchos::ArrayView< GlobalOrdinal > &Indices, const Teuchos::ArrayView< Scalar > &Values, size_t &NumEntries) const =0
Get a copy of the given global row&#39;s entries.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
Sets up and executes a communication plan for a Tpetra DistObject.
virtual void pack(const Teuchos::ArrayView< const LocalOrdinal > &exportLIDs, Teuchos::Array< char > &exports, const Teuchos::ArrayView< size_t > &numPacketsPerLID, size_t &constantNumPackets, Distributor &distor) const
Pack this object&#39;s data for an Import or Export.
virtual Teuchos::RCP< RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > add(const Scalar &alpha, const RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta, const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &domainMap=Teuchos::null, const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rangeMap=Teuchos::null, const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null) const
Return a new RowMatrix which is the result of beta*this + alpha*A.
A parallel distribution of indices over processes.
A read-only, row-oriented interface to a sparse matrix.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
The Map that describes the distribution of rows over processes.