157 "Map, then you must supply a nonnull domain and range Map to this "
168 if (!
params.is_null ()) {
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;
177 RCP<crs_matrix_type> C;
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);
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;
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;
203 if (constructorSublist.is_null ()) {
204 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow ()));
206 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow (),
207 constructorSublist));
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");
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.");
233 using gids_type = nonconst_global_inds_host_view_type;
234 using vals_type = nonconst_values_host_view_type;
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);
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);
251 if (alpha != STS::one ()) {
252 for (
size_t k = 0; k < A_numEntries; ++k) {
256 C->insertGlobalValues (globalRow, A_numEntries,
257 reinterpret_cast<const Scalar*
>(valView.data()),
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);
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);
275 if (beta != STS::one ()) {
276 for (
size_t k = 0; k < B_numEntries; ++k) {
280 C->insertGlobalValues (globalRow, B_numEntries,
281 reinterpret_cast<const Scalar*
>(valView.data()),
286 if (callFillComplete) {
287 if (fillCompleteSublist.is_null ()) {
288 C->fillComplete (theDomainMap, theRangeMap);
290 C->fillComplete (theDomainMap, theRangeMap, fillCompleteSublist);
294 return rcp_implicit_cast<this_type> (C);
298 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
302 Teuchos::Array<char>& exports,
306#ifdef HAVE_TPETRA_DEBUG
309 using Teuchos::reduceAll;
310 std::ostringstream
msg;
315 }
catch (std::exception&
e) {
320 const Teuchos::Comm<int>& comm = * (this->getComm ());
324 const int myRank = comm.getRank ();
325 const int numProcs = comm.getSize ();
328 std::ostringstream
os;
329 os <<
"Proc " <<
myRank <<
": " <<
msg.str () << std::endl;
330 std::cerr <<
os.str ();
337 true, std::logic_error,
"packImpl() threw an exception on one or "
338 "more participating processes.");
347 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
351 size_t& totalNumEntries,
352 const Teuchos::ArrayView<const LocalOrdinal>&
exportLIDs)
const
356 typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
367 if (
curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid ()) {
370 totalNumEntries += curNumEntries;
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);
389 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
391 RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
392 packRow (
char*
const numEntOut,
396 const LocalOrdinal lclRow)
const
398 using Teuchos::Array;
399 using Teuchos::ArrayView;
400 typedef LocalOrdinal LO;
401 typedef GlobalOrdinal GO;
404 const LO numEntLO =
static_cast<LO
> (numEnt);
405 memcpy (numEntOut, &numEntLO,
sizeof (LO));
407 if (this->supportsRowViews ()) {
408 if (this->isLocallyIndexed ()) {
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 ());
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));
422 memcpy (valOut, valIn.data (), numEnt * sizeof (Scalar));
424 else if (this->isGloballyIndexed ()) {
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));
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) {
455 const map_type&
colMap = * (this->getColMap ());
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));
462 memcpy (valOut, valIn.data(), numEnt * sizeof (Scalar));
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) {
474 memcpy (indOut, indIn.data(), numEnt * sizeof (GO));
475 memcpy (valOut, valIn.data(), numEnt * sizeof (Scalar));
486 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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
494 using Teuchos::Array;
495 using Teuchos::ArrayView;
497 using Teuchos::av_reinterpret_cast;
499 typedef LocalOrdinal LO;
500 typedef GlobalOrdinal GO;
501 typedef typename ArrayView<const LO>::size_type size_type;
502 const char tfecfFuncName[] =
"packImpl: ";
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 () <<
".");
513 constantNumPackets = 0;
518 size_t totalNumEntries = 0;
519 allocatePackSpace (exports, totalNumEntries, exportLIDs);
520 const size_t bufSize =
static_cast<size_t> (exports.size ());
535 bool outOfBounds =
false;
542 const size_t numEnt = this->getNumEntriesInLocalRow (
lclRow);
554 const size_t numBytes =
sizeof (LO) +
585 outOfBounds, std::logic_error,
"First invalid offset into 'exports' "
586 "pack buffer at index i = " <<
firstBadIndex <<
". exportLIDs[i]: "
590 packErr, std::logic_error,
"First error in packRow() at index i = "
605#define TPETRA_ROWMATRIX_INSTANT(SCALAR,LO,GO,NODE) \
606 template class RowMatrix< SCALAR , LO , GO , NODE >;