Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
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
49namespace 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->getLocalNumElements ());
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 } else {
206 C = rcp (new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow (),
207 constructorSublist));
208 }
209 // Since A and B have the same row Maps, we could add them
210 // together all at once and merge values before we call
211 // insertGlobalValues. However, we don't really need to, since
212 // we've already allocated enough space in each row of C for C
213 // to do the merge itself.
214 }
215 else { // the row Maps of A and B are not the same
216 // Construct the result matrix C.
217 // true: !A_rowMap->isSameAs (*B_rowMap)
218 TEUCHOS_TEST_FOR_EXCEPTION(true,
219 std::invalid_argument,
220 "Tpetra::RowMatrix::add: The row maps must be the same for statically "
221 "allocated matrices in order to be sure that there is sufficient space "
222 "to do the addition");
223 }
224
225#ifdef HAVE_TPETRA_DEBUG
226 TEUCHOS_TEST_FOR_EXCEPTION(C.is_null (), std::logic_error,
227 "Tpetra::RowMatrix::add: C should not be null at this point. "
228 "Please report this bug to the Tpetra developers.");
229#endif // HAVE_TPETRA_DEBUG
230 //
231 // Compute C = alpha*A + beta*B.
232 //
233 using gids_type = nonconst_global_inds_host_view_type;
234 using vals_type = nonconst_values_host_view_type;
235 gids_type ind;
236 vals_type val;
237
238 if (alpha != STS::zero ()) {
239 const LO A_localNumRows = static_cast<LO> (A_rowMap->getLocalNumElements ());
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 Kokkos::resize(ind,A_numEntries);
245 Kokkos::resize(val,A_numEntries);
246 }
247 gids_type indView = Kokkos::subview(ind, std::make_pair((size_t)0, A_numEntries));
248 vals_type valView = Kokkos::subview(val, std::make_pair((size_t)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, A_numEntries,
257 reinterpret_cast<const Scalar*>(valView.data()),
258 indView.data());
259 }
260 }
261
262 if (beta != STS::zero ()) {
263 const LO B_localNumRows = static_cast<LO> (B_rowMap->getLocalNumElements ());
264 for (LO localRow = 0; localRow < B_localNumRows; ++localRow) {
265 size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
266 const GO globalRow = B_rowMap->getGlobalElement (localRow);
267 if (B_numEntries > static_cast<size_t> (ind.size ())) {
268 Kokkos::resize(ind,B_numEntries);
269 Kokkos::resize(val,B_numEntries);
270 }
271 gids_type indView = Kokkos::subview(ind, std::make_pair((size_t)0, B_numEntries));
272 vals_type valView = Kokkos::subview(val, std::make_pair((size_t)0, B_numEntries));
273 B.getGlobalRowCopy (globalRow, indView, valView, B_numEntries);
274
275 if (beta != STS::one ()) {
276 for (size_t k = 0; k < B_numEntries; ++k) {
277 valView[k] *= beta;
278 }
279 }
280 C->insertGlobalValues (globalRow, B_numEntries,
281 reinterpret_cast<const Scalar*>(valView.data()),
282 indView.data());
283 }
284 }
285
286 if (callFillComplete) {
287 if (fillCompleteSublist.is_null ()) {
288 C->fillComplete (theDomainMap, theRangeMap);
289 } else {
290 C->fillComplete (theDomainMap, theRangeMap, fillCompleteSublist);
291 }
292 }
293
294 return rcp_implicit_cast<this_type> (C);
295 }
296
297
298 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
299 void
301 pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
302 Teuchos::Array<char>& exports,
303 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
304 size_t& constantNumPackets) const
305 {
306#ifdef HAVE_TPETRA_DEBUG
307 const char tfecfFuncName[] = "pack: ";
308 {
309 using Teuchos::reduceAll;
310 std::ostringstream msg;
311 int lclBad = 0;
312 try {
313 this->packImpl (exportLIDs, exports, numPacketsPerLID,
314 constantNumPackets);
315 } catch (std::exception& e) {
316 lclBad = 1;
317 msg << e.what ();
318 }
319 int gblBad = 0;
320 const Teuchos::Comm<int>& comm = * (this->getComm ());
321 reduceAll<int, int> (comm, Teuchos::REDUCE_MAX,
322 lclBad, Teuchos::outArg (gblBad));
323 if (gblBad != 0) {
324 const int myRank = comm.getRank ();
325 const int numProcs = comm.getSize ();
326 for (int r = 0; r < numProcs; ++r) {
327 if (r == myRank && lclBad != 0) {
328 std::ostringstream os;
329 os << "Proc " << myRank << ": " << msg.str () << std::endl;
330 std::cerr << os.str ();
331 }
332 comm.barrier ();
333 comm.barrier ();
334 comm.barrier ();
335 }
336 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
337 true, std::logic_error, "packImpl() threw an exception on one or "
338 "more participating processes.");
339 }
340 }
341#else
342 this->packImpl (exportLIDs, exports, numPacketsPerLID,
343 constantNumPackets);
344#endif // HAVE_TPETRA_DEBUG
345 }
346
347 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
348 void
350 allocatePackSpace (Teuchos::Array<char>& exports,
351 size_t& totalNumEntries,
352 const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs) const
353 {
354 typedef LocalOrdinal LO;
355 typedef GlobalOrdinal GO;
356 typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
357 //const char tfecfFuncName[] = "allocatePackSpace: ";
358 const size_type numExportLIDs = exportLIDs.size ();
359
360 // Count the total number of entries to send.
361 totalNumEntries = 0;
362 for (size_type i = 0; i < numExportLIDs; ++i) {
363 const LO lclRow = exportLIDs[i];
364 size_t curNumEntries = this->getNumEntriesInLocalRow (lclRow);
365 // FIXME (mfh 25 Jan 2015) We should actually report invalid row
366 // indices as an error. Just consider them nonowned for now.
367 if (curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid ()) {
368 curNumEntries = 0;
369 }
370 totalNumEntries += curNumEntries;
371 }
372
373 // FIXME (mfh 24 Feb 2013) This code is only correct if
374 // sizeof(Scalar) is a meaningful representation of the amount of
375 // data in a Scalar instance. (LO and GO are always built-in
376 // integer types.)
377 //
378 // Allocate the exports array. It does NOT need padding for
379 // alignment, since we use memcpy to write to / read from send /
380 // receive buffers.
381 const size_t allocSize =
382 static_cast<size_t> (numExportLIDs) * sizeof (LO) +
383 totalNumEntries * (sizeof (Scalar) + sizeof (GO));
384 if (static_cast<size_t> (exports.size ()) < allocSize) {
385 exports.resize (allocSize);
386 }
387 }
388
389 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
390 bool
391 RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
392 packRow (char* const numEntOut,
393 char* const valOut,
394 char* const indOut,
395 const size_t numEnt,
396 const LocalOrdinal lclRow) const
397 {
398 using Teuchos::Array;
399 using Teuchos::ArrayView;
400 typedef LocalOrdinal LO;
401 typedef GlobalOrdinal GO;
402 typedef Tpetra::Map<LO, GO, Node> map_type;
403
404 const LO numEntLO = static_cast<LO> (numEnt);
405 memcpy (numEntOut, &numEntLO, sizeof (LO));
406
407 if (this->supportsRowViews ()) {
408 if (this->isLocallyIndexed ()) {
409 // If the matrix is locally indexed on the calling process, we
410 // have to use its column Map (which it _must_ have in this
411 // case) to convert to global indices.
412 local_inds_host_view_type indIn;
413 values_host_view_type valIn;
414 this->getLocalRowView (lclRow, indIn, valIn);
415 const map_type& colMap = * (this->getColMap ());
416 // Copy column indices one at a time, so that we don't need
417 // temporary storage.
418 for (size_t k = 0; k < numEnt; ++k) {
419 const GO gblIndIn = colMap.getGlobalElement (indIn[k]);
420 memcpy (indOut + k * sizeof (GO), &gblIndIn, sizeof (GO));
421 }
422 memcpy (valOut, valIn.data (), numEnt * sizeof (Scalar));
423 }
424 else if (this->isGloballyIndexed ()) {
425 // If the matrix is globally indexed on the calling process,
426 // then we can use the column indices directly. However, we
427 // have to get the global row index. The calling process must
428 // have a row Map, since otherwise it shouldn't be participating
429 // in packing operations.
430 global_inds_host_view_type indIn;
431 values_host_view_type valIn;
432 const map_type& rowMap = * (this->getRowMap ());
433 const GO gblRow = rowMap.getGlobalElement (lclRow);
434 this->getGlobalRowView (gblRow, indIn, valIn);
435 memcpy (indOut, indIn.data (), numEnt * sizeof (GO));
436 memcpy (valOut, valIn.data (), numEnt * sizeof (Scalar));
437 }
438 else {
439 if (numEnt != 0) {
440 return false;
441 }
442 }
443 }
444 else {
445 // FIXME (mfh 25 Jan 2015) Pass in valIn and indIn as scratch
446 // space, instead of allocating them on each call.
447 if (this->isLocallyIndexed ()) {
448 nonconst_local_inds_host_view_type indIn("indIn",numEnt);
449 nonconst_values_host_view_type valIn("valIn",numEnt);
450 size_t theNumEnt = 0;
451 this->getLocalRowCopy (lclRow, indIn, valIn, theNumEnt);
452 if (theNumEnt != numEnt) {
453 return false;
454 }
455 const map_type& colMap = * (this->getColMap ());
456 // Copy column indices one at a time, so that we don't need
457 // temporary storage.
458 for (size_t k = 0; k < numEnt; ++k) {
459 const GO gblIndIn = colMap.getGlobalElement (indIn[k]);
460 memcpy (indOut + k * sizeof (GO), &gblIndIn, sizeof (GO));
461 }
462 memcpy (valOut, valIn.data(), numEnt * sizeof (Scalar));
463 }
464 else if (this->isGloballyIndexed ()) {
465 nonconst_global_inds_host_view_type indIn("indIn",numEnt);
466 nonconst_values_host_view_type valIn("valIn",numEnt);
467 const map_type& rowMap = * (this->getRowMap ());
468 const GO gblRow = rowMap.getGlobalElement (lclRow);
469 size_t theNumEnt = 0;
470 this->getGlobalRowCopy (gblRow, indIn, valIn, theNumEnt);
471 if (theNumEnt != numEnt) {
472 return false;
473 }
474 memcpy (indOut, indIn.data(), numEnt * sizeof (GO));
475 memcpy (valOut, valIn.data(), numEnt * sizeof (Scalar));
476 }
477 else {
478 if (numEnt != 0) {
479 return false;
481 }
482 }
483 return true;
484 }
485
486 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
487 void
488 RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
489 packImpl (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
490 Teuchos::Array<char>& exports,
491 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
492 size_t& constantNumPackets) const
493 {
494 using Teuchos::Array;
495 using Teuchos::ArrayView;
496 using Teuchos::as;
497 using Teuchos::av_reinterpret_cast;
498 using Teuchos::RCP;
499 typedef LocalOrdinal LO;
500 typedef GlobalOrdinal GO;
501 typedef typename ArrayView<const LO>::size_type size_type;
502 const char tfecfFuncName[] = "packImpl: ";
503
504 const size_type numExportLIDs = exportLIDs.size ();
505 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
506 numExportLIDs != numPacketsPerLID.size (), std::invalid_argument,
507 "exportLIDs.size() = " << numExportLIDs << " != numPacketsPerLID.size()"
508 " = " << numPacketsPerLID.size () << ".");
509
510 // Setting this to zero tells the caller to expect a possibly
511 // different ("nonconstant") number of packets per local index
512 // (i.e., a possibly different number of entries per row).
513 constantNumPackets = 0;
514
515 // The pack buffer 'exports' enters this method possibly
516 // unallocated. Do the first two parts of "Count, allocate, fill,
517 // compute."
518 size_t totalNumEntries = 0;
519 allocatePackSpace (exports, totalNumEntries, exportLIDs);
520 const size_t bufSize = static_cast<size_t> (exports.size ());
521
522 // Compute the number of "packets" (in this case, bytes) per
523 // export LID (in this case, local index of the row to send), and
524 // actually pack the data.
525 //
526 // FIXME (mfh 24 Feb 2013, 25 Jan 2015) This code is only correct
527 // if sizeof(Scalar) is a meaningful representation of the amount
528 // of data in a Scalar instance. (LO and GO are always built-in
529 // integer types.)
530
531 // Variables for error reporting in the loop.
532 size_type firstBadIndex = 0; // only valid if outOfBounds == true.
533 size_t firstBadOffset = 0; // only valid if outOfBounds == true.
534 size_t firstBadNumBytes = 0; // only valid if outOfBounds == true.
535 bool outOfBounds = false;
536 bool packErr = false;
537
538 char* const exportsRawPtr = exports.getRawPtr ();
539 size_t offset = 0; // current index into 'exports' array.
540 for (size_type i = 0; i < numExportLIDs; ++i) {
541 const LO lclRow = exportLIDs[i];
542 const size_t numEnt = this->getNumEntriesInLocalRow (lclRow);
543
544 // Only pad this row if it has a nonzero number of entries.
545 if (numEnt == 0) {
546 numPacketsPerLID[i] = 0;
547 }
548 else {
549 char* const numEntBeg = exportsRawPtr + offset;
550 char* const numEntEnd = numEntBeg + sizeof (LO);
551 char* const valBeg = numEntEnd;
552 char* const valEnd = valBeg + numEnt * sizeof (Scalar);
553 char* const indBeg = valEnd;
554 const size_t numBytes = sizeof (LO) +
555 numEnt * (sizeof (Scalar) + sizeof (GO));
556 if (offset > bufSize || offset + numBytes > bufSize) {
557 firstBadIndex = i;
558 firstBadOffset = offset;
559 firstBadNumBytes = numBytes;
560 outOfBounds = true;
561 break;
562 }
563 packErr = ! packRow (numEntBeg, valBeg, indBeg, numEnt, lclRow);
564 if (packErr) {
565 firstBadIndex = i;
566 firstBadOffset = offset;
567 firstBadNumBytes = numBytes;
568 break;
569 }
570 // numPacketsPerLID[i] is the number of "packets" in the
571 // current local row i. Packet=char (really "byte") so use
572 // the number of bytes of the packed data for that row.
573 numPacketsPerLID[i] = numBytes;
574 offset += numBytes;
575 }
576 }
577
578 // The point of these exceptions is to stop computation if the
579 // above checks found a bug. If HAVE_TPETRA_DEBUG is defined,
580 // Tpetra will do extra all-reduces to check for global
581 // consistency of the error state. Otherwise, throwing an
582 // exception here might cause deadlock, but we accept that as
583 // better than just continuing.
584 TEUCHOS_TEST_FOR_EXCEPTION(
585 outOfBounds, std::logic_error, "First invalid offset into 'exports' "
586 "pack buffer at index i = " << firstBadIndex << ". exportLIDs[i]: "
587 << exportLIDs[firstBadIndex] << ", bufSize: " << bufSize << ", offset: "
588 << firstBadOffset << ", numBytes: " << firstBadNumBytes << ".");
589 TEUCHOS_TEST_FOR_EXCEPTION(
590 packErr, std::logic_error, "First error in packRow() at index i = "
591 << firstBadIndex << ". exportLIDs[i]: " << exportLIDs[firstBadIndex]
592 << ", bufSize: " << bufSize << ", offset: " << firstBadOffset
593 << ", numBytes: " << firstBadNumBytes << ".");
594 }
595
596
597} // namespace Tpetra
598
599//
600// Explicit instantiation macro
601//
602// Must be expanded from within the Tpetra namespace!
603//
604
605#define TPETRA_ROWMATRIX_INSTANT(SCALAR,LO,GO,NODE) \
606 template class RowMatrix< SCALAR , LO , GO , NODE >;
607
608
609#endif // TPETRA_ROWMATRIX_DEF_HPP
610
Sparse matrix that presents a row-oriented interface that lets users read or modify entries.
A parallel distribution of indices over processes.
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 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....
A read-only, row-oriented interface to a sparse matrix.
virtual void pack(const Teuchos::ArrayView< const LocalOrdinal > &exportLIDs, Teuchos::Array< char > &exports, const Teuchos::ArrayView< size_t > &numPacketsPerLID, size_t &constantNumPackets) const
Pack this object's data for an Import or Export.
virtual ~RowMatrix()
Destructor (virtual for memory safety of derived classes).
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.
Namespace Tpetra contains the class and methods constituting the Tpetra library.