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