42 #ifndef TPETRA_ROWMATRIX_DEF_HPP
43 #define TPETRA_ROWMATRIX_DEF_HPP
45 #include "Tpetra_CrsMatrix.hpp"
46 #include "Tpetra_Map.hpp"
47 #include "Tpetra_RowGraph.hpp"
51 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
54 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
55 Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
57 add (
const Scalar& alpha,
62 const Teuchos::RCP<Teuchos::ParameterList>& params)
const
65 using Teuchos::ArrayView;
66 using Teuchos::ParameterList;
69 using Teuchos::rcp_implicit_cast;
70 using Teuchos::sublist;
71 typedef LocalOrdinal LO;
72 typedef GlobalOrdinal GO;
73 typedef Teuchos::ScalarTraits<Scalar> STS;
78 const this_type& B = *
this;
86 RCP<const map_type> B_domainMap = B.getDomainMap ();
87 RCP<const map_type> B_rangeMap = B.getRangeMap ();
89 RCP<const map_type> theDomainMap = domainMap;
90 RCP<const map_type> theRangeMap = rangeMap;
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;
100 theDomainMap = B_domainMap;
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;
111 theRangeMap = B_rangeMap;
115 #ifdef HAVE_TPETRA_DEBUG
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.");
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.");
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 "
160 #endif // HAVE_TPETRA_DEBUG
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");
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->getNodeNumElements ());
186 Array<size_t> C_maxNumEntriesPerRow (localNumRows, 0);
189 if (alpha != STS::zero ()) {
190 for (LO localRow = 0; localRow < localNumRows; ++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 (),
207 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow (),
208 StaticProfile, constructorSublist));
218 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
219 if (constructorSublist.is_null ()) {
220 C = rcp (
new crs_matrix_type (C_rowMap, 0, DynamicProfile));
222 C = rcp (
new crs_matrix_type (C_rowMap, 0, DynamicProfile,
223 constructorSublist));
228 #ifdef HAVE_TPETRA_DEBUG
229 TEUCHOS_TEST_FOR_EXCEPTION(C.is_null (), std::logic_error,
230 "Tpetra::RowMatrix::add: C should not be null at this point. "
231 "Please report this bug to the Tpetra developers.");
232 #endif // HAVE_TPETRA_DEBUG
239 if (alpha != STS::zero ()) {
240 const LO A_localNumRows =
static_cast<LO
> (A_rowMap->getNodeNumElements ());
241 for (LO localRow = 0; localRow < A_localNumRows; ++localRow) {
243 const GO globalRow = A_rowMap->getGlobalElement (localRow);
244 if (A_numEntries > static_cast<size_t> (ind.size ())) {
245 ind.resize (A_numEntries);
246 val.resize (A_numEntries);
248 ArrayView<GO> indView = ind (0, A_numEntries);
249 ArrayView<Scalar> valView = val (0, A_numEntries);
252 if (alpha != STS::one ()) {
253 for (
size_t k = 0; k < A_numEntries; ++k) {
257 C->insertGlobalValues (globalRow, indView, valView);
261 if (beta != STS::zero ()) {
262 const LO B_localNumRows =
static_cast<LO
> (B_rowMap->getNodeNumElements ());
263 for (LO localRow = 0; localRow < B_localNumRows; ++localRow) {
264 size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
265 const GO globalRow = B_rowMap->getGlobalElement (localRow);
266 if (B_numEntries > static_cast<size_t> (ind.size ())) {
267 ind.resize (B_numEntries);
268 val.resize (B_numEntries);
270 ArrayView<GO> indView = ind (0, B_numEntries);
271 ArrayView<Scalar> valView = val (0, B_numEntries);
272 B.getGlobalRowCopy (globalRow, indView, valView, B_numEntries);
274 if (beta != STS::one ()) {
275 for (
size_t k = 0; k < B_numEntries; ++k) {
279 C->insertGlobalValues (globalRow, indView, valView);
283 if (callFillComplete) {
284 if (fillCompleteSublist.is_null ()) {
285 C->fillComplete (theDomainMap, theRangeMap);
287 C->fillComplete (theDomainMap, theRangeMap, fillCompleteSublist);
291 return rcp_implicit_cast<this_type> (C);
295 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
298 pack (
const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
299 Teuchos::Array<char>& exports,
300 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
301 size_t& constantNumPackets,
304 #ifdef HAVE_TPETRA_DEBUG
305 const char tfecfFuncName[] =
"pack: ";
307 using Teuchos::reduceAll;
308 std::ostringstream msg;
311 this->packImpl (exportLIDs, exports, numPacketsPerLID,
312 constantNumPackets, distor);
313 }
catch (std::exception& e) {
318 const Teuchos::Comm<int>& comm = * (this->getComm ());
319 reduceAll<int, int> (comm, Teuchos::REDUCE_MAX,
320 lclBad, Teuchos::outArg (gblBad));
322 const int myRank = comm.getRank ();
323 const int numProcs = comm.getSize ();
324 for (
int r = 0; r < numProcs; ++r) {
325 if (r == myRank && lclBad != 0) {
326 std::ostringstream os;
327 os <<
"Proc " << myRank <<
": " << msg.str () << std::endl;
328 std::cerr << os.str ();
334 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
335 true, std::logic_error,
"packImpl() threw an exception on one or "
336 "more participating processes.");
340 this->packImpl (exportLIDs, exports, numPacketsPerLID,
341 constantNumPackets, distor);
342 #endif // HAVE_TPETRA_DEBUG
345 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
349 size_t& totalNumEntries,
350 const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs)
const
352 typedef LocalOrdinal LO;
353 typedef GlobalOrdinal GO;
354 typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
356 const size_type numExportLIDs = exportLIDs.size ();
360 for (size_type i = 0; i < numExportLIDs; ++i) {
361 const LO lclRow = exportLIDs[i];
362 size_t curNumEntries = this->getNumEntriesInLocalRow (lclRow);
365 if (curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid ()) {
368 totalNumEntries += curNumEntries;
379 const size_t allocSize =
380 static_cast<size_t> (numExportLIDs) *
sizeof (LO) +
381 totalNumEntries * (
sizeof (Scalar) +
sizeof (GO));
382 if (static_cast<size_t> (exports.size ()) < allocSize) {
383 exports.resize (allocSize);
387 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
389 RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
390 packRow (
char*
const numEntOut,
394 const LocalOrdinal lclRow)
const
396 using Teuchos::Array;
397 using Teuchos::ArrayView;
398 typedef LocalOrdinal LO;
399 typedef GlobalOrdinal GO;
402 const LO numEntLO =
static_cast<LO
> (numEnt);
403 memcpy (numEntOut, &numEntLO,
sizeof (LO));
405 if (this->supportsRowViews ()) {
406 if (this->isLocallyIndexed ()) {
410 ArrayView<const LO> indIn;
411 ArrayView<const Scalar> valIn;
412 this->getLocalRowView (lclRow, indIn, valIn);
413 const map_type& colMap = * (this->getColMap ());
416 for (
size_t k = 0; k < numEnt; ++k) {
417 const GO gblIndIn = colMap.getGlobalElement (indIn[k]);
418 memcpy (indOut + k *
sizeof (GO), &gblIndIn,
sizeof (GO));
420 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
422 else if (this->isGloballyIndexed ()) {
428 ArrayView<const GO> indIn;
429 ArrayView<const Scalar> valIn;
430 const map_type& rowMap = * (this->getRowMap ());
431 const GO gblRow = rowMap.getGlobalElement (lclRow);
432 this->getGlobalRowView (gblRow, indIn, valIn);
433 memcpy (indOut, indIn.getRawPtr (), numEnt *
sizeof (GO));
434 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
445 if (this->isLocallyIndexed ()) {
446 Array<LO> indIn (numEnt);
447 Array<Scalar> valIn (numEnt);
448 size_t theNumEnt = 0;
449 this->getLocalRowCopy (lclRow, indIn (), valIn (), theNumEnt);
450 if (theNumEnt != numEnt) {
453 const map_type& colMap = * (this->getColMap ());
456 for (
size_t k = 0; k < numEnt; ++k) {
457 const GO gblIndIn = colMap.getGlobalElement (indIn[k]);
458 memcpy (indOut + k *
sizeof (GO), &gblIndIn,
sizeof (GO));
460 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
462 else if (this->isGloballyIndexed ()) {
463 Array<GO> indIn (numEnt);
464 Array<Scalar> valIn (numEnt);
465 const map_type& rowMap = * (this->getRowMap ());
466 const GO gblRow = rowMap.getGlobalElement (lclRow);
467 size_t theNumEnt = 0;
468 this->getGlobalRowCopy (gblRow, indIn, valIn, theNumEnt);
469 if (theNumEnt != numEnt) {
472 memcpy (indOut, indIn.getRawPtr (), numEnt *
sizeof (GO));
473 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
484 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
486 RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
487 packImpl (
const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
488 Teuchos::Array<char>& exports,
489 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
490 size_t& constantNumPackets,
493 using Teuchos::Array;
494 using Teuchos::ArrayView;
496 using Teuchos::av_reinterpret_cast;
498 typedef LocalOrdinal LO;
499 typedef GlobalOrdinal GO;
500 typedef typename ArrayView<const LO>::size_type size_type;
501 const char tfecfFuncName[] =
"packImpl: ";
503 const size_type numExportLIDs = exportLIDs.size ();
504 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
505 numExportLIDs != numPacketsPerLID.size (), std::invalid_argument,
506 "exportLIDs.size() = " << numExportLIDs <<
" != numPacketsPerLID.size()"
507 " = " << numPacketsPerLID.size () <<
".");
512 constantNumPackets = 0;
517 size_t totalNumEntries = 0;
518 allocatePackSpace (exports, totalNumEntries, exportLIDs);
519 const size_t bufSize =
static_cast<size_t> (exports.size ());
531 size_type firstBadIndex = 0;
532 size_t firstBadOffset = 0;
533 size_t firstBadNumBytes = 0;
535 bool packErr =
false;
537 char*
const exportsRawPtr = exports.getRawPtr ();
539 for (size_type i = 0; i < numExportLIDs; ++i) {
540 const LO lclRow = exportLIDs[i];
541 const size_t numEnt = this->getNumEntriesInLocalRow (lclRow);
545 numPacketsPerLID[i] = 0;
548 char*
const numEntBeg = exportsRawPtr + offset;
549 char*
const numEntEnd = numEntBeg +
sizeof (LO);
550 char*
const valBeg = numEntEnd;
551 char*
const valEnd = valBeg + numEnt *
sizeof (Scalar);
552 char*
const indBeg = valEnd;
553 const size_t numBytes =
sizeof (LO) +
554 numEnt * (
sizeof (Scalar) +
sizeof (GO));
555 if (offset > bufSize || offset + numBytes > bufSize) {
557 firstBadOffset = offset;
558 firstBadNumBytes = numBytes;
562 packErr = ! packRow (numEntBeg, valBeg, indBeg, numEnt, lclRow);
565 firstBadOffset = offset;
566 firstBadNumBytes = numBytes;
572 numPacketsPerLID[i] = numBytes;
583 TEUCHOS_TEST_FOR_EXCEPTION(
584 outOfBounds, std::logic_error,
"First invalid offset into 'exports' "
585 "pack buffer at index i = " << firstBadIndex <<
". exportLIDs[i]: "
586 << exportLIDs[firstBadIndex] <<
", bufSize: " << bufSize <<
", offset: "
587 << firstBadOffset <<
", numBytes: " << firstBadNumBytes <<
".");
588 TEUCHOS_TEST_FOR_EXCEPTION(
589 packErr, std::logic_error,
"First error in packRow() at index i = "
590 << firstBadIndex <<
". exportLIDs[i]: " << exportLIDs[firstBadIndex]
591 <<
", bufSize: " << bufSize <<
", offset: " << firstBadOffset
592 <<
", numBytes: " << firstBadNumBytes <<
".");
595 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
599 LocalOrdinal& numEnt,
600 const LocalOrdinal*& lclColInds,
601 const Scalar*& vals)
const
606 Teuchos::ArrayView<const LocalOrdinal> lclColInds_av;
607 Teuchos::ArrayView<const Scalar> vals_av;
609 this->getLocalRowView (lclRow, lclColInds_av, vals_av);
610 numEnt =
static_cast<LocalOrdinal
> (lclColInds_av.size ());
616 lclColInds = lclColInds_av.getRawPtr ();
617 vals = vals_av.getRawPtr ();
620 return static_cast<LocalOrdinal
> (0);
631 #define TPETRA_ROWMATRIX_INSTANT(SCALAR,LO,GO,NODE) \
632 template class RowMatrix< SCALAR , LO , GO , NODE >;
635 #endif // TPETRA_ROWMATRIX_DEF_HPP
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 "raw" 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'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'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 > ¶ms=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.