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,
ProfileType(StaticProfile+1) ));
222 C = rcp (
new crs_matrix_type (C_rowMap, 0,
ProfileType(StaticProfile+1) ,
223 constructorSublist));
227 TEUCHOS_TEST_FOR_EXCEPTION(
true,
228 std::invalid_argument,
229 "Tpetra::RowMatrix::add: The row maps must be the same for statically "
230 "allocated matrices in order to be sure that there is sufficient space "
231 "to do the addition");
235 #ifdef HAVE_TPETRA_DEBUG
236 TEUCHOS_TEST_FOR_EXCEPTION(C.is_null (), std::logic_error,
237 "Tpetra::RowMatrix::add: C should not be null at this point. "
238 "Please report this bug to the Tpetra developers.");
239 #endif // HAVE_TPETRA_DEBUG
246 if (alpha != STS::zero ()) {
247 const LO A_localNumRows =
static_cast<LO
> (A_rowMap->getNodeNumElements ());
248 for (LO localRow = 0; localRow < A_localNumRows; ++localRow) {
250 const GO globalRow = A_rowMap->getGlobalElement (localRow);
251 if (A_numEntries > static_cast<size_t> (ind.size ())) {
252 ind.resize (A_numEntries);
253 val.resize (A_numEntries);
255 ArrayView<GO> indView = ind (0, A_numEntries);
256 ArrayView<Scalar> valView = val (0, A_numEntries);
259 if (alpha != STS::one ()) {
260 for (
size_t k = 0; k < A_numEntries; ++k) {
264 C->insertGlobalValues (globalRow, indView, valView);
268 if (beta != STS::zero ()) {
269 const LO B_localNumRows =
static_cast<LO
> (B_rowMap->getNodeNumElements ());
270 for (LO localRow = 0; localRow < B_localNumRows; ++localRow) {
271 size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
272 const GO globalRow = B_rowMap->getGlobalElement (localRow);
273 if (B_numEntries > static_cast<size_t> (ind.size ())) {
274 ind.resize (B_numEntries);
275 val.resize (B_numEntries);
277 ArrayView<GO> indView = ind (0, B_numEntries);
278 ArrayView<Scalar> valView = val (0, B_numEntries);
279 B.getGlobalRowCopy (globalRow, indView, valView, B_numEntries);
281 if (beta != STS::one ()) {
282 for (
size_t k = 0; k < B_numEntries; ++k) {
286 C->insertGlobalValues (globalRow, indView, valView);
290 if (callFillComplete) {
291 if (fillCompleteSublist.is_null ()) {
292 C->fillComplete (theDomainMap, theRangeMap);
294 C->fillComplete (theDomainMap, theRangeMap, fillCompleteSublist);
298 return rcp_implicit_cast<this_type> (C);
302 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
305 pack (
const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
306 Teuchos::Array<char>& exports,
307 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
308 size_t& constantNumPackets,
311 #ifdef HAVE_TPETRA_DEBUG
312 const char tfecfFuncName[] =
"pack: ";
314 using Teuchos::reduceAll;
315 std::ostringstream msg;
318 this->packImpl (exportLIDs, exports, numPacketsPerLID,
319 constantNumPackets, distor);
320 }
catch (std::exception& e) {
325 const Teuchos::Comm<int>& comm = * (this->getComm ());
326 reduceAll<int, int> (comm, Teuchos::REDUCE_MAX,
327 lclBad, Teuchos::outArg (gblBad));
329 const int myRank = comm.getRank ();
330 const int numProcs = comm.getSize ();
331 for (
int r = 0; r < numProcs; ++r) {
332 if (r == myRank && lclBad != 0) {
333 std::ostringstream os;
334 os <<
"Proc " << myRank <<
": " << msg.str () << std::endl;
335 std::cerr << os.str ();
341 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
342 true, std::logic_error,
"packImpl() threw an exception on one or "
343 "more participating processes.");
347 this->packImpl (exportLIDs, exports, numPacketsPerLID,
348 constantNumPackets, distor);
349 #endif // HAVE_TPETRA_DEBUG
352 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
356 size_t& totalNumEntries,
357 const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs)
const
359 typedef LocalOrdinal LO;
360 typedef GlobalOrdinal GO;
361 typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
363 const size_type numExportLIDs = exportLIDs.size ();
367 for (size_type i = 0; i < numExportLIDs; ++i) {
368 const LO lclRow = exportLIDs[i];
369 size_t curNumEntries = this->getNumEntriesInLocalRow (lclRow);
372 if (curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid ()) {
375 totalNumEntries += curNumEntries;
386 const size_t allocSize =
387 static_cast<size_t> (numExportLIDs) *
sizeof (LO) +
388 totalNumEntries * (
sizeof (Scalar) +
sizeof (GO));
389 if (static_cast<size_t> (exports.size ()) < allocSize) {
390 exports.resize (allocSize);
394 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
396 RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
397 packRow (
char*
const numEntOut,
401 const LocalOrdinal lclRow)
const
403 using Teuchos::Array;
404 using Teuchos::ArrayView;
405 typedef LocalOrdinal LO;
406 typedef GlobalOrdinal GO;
409 const LO numEntLO =
static_cast<LO
> (numEnt);
410 memcpy (numEntOut, &numEntLO,
sizeof (LO));
412 if (this->supportsRowViews ()) {
413 if (this->isLocallyIndexed ()) {
417 ArrayView<const LO> indIn;
418 ArrayView<const Scalar> valIn;
419 this->getLocalRowView (lclRow, indIn, valIn);
420 const map_type& colMap = * (this->getColMap ());
423 for (
size_t k = 0; k < numEnt; ++k) {
424 const GO gblIndIn = colMap.getGlobalElement (indIn[k]);
425 memcpy (indOut + k *
sizeof (GO), &gblIndIn,
sizeof (GO));
427 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
429 else if (this->isGloballyIndexed ()) {
435 ArrayView<const GO> indIn;
436 ArrayView<const Scalar> valIn;
437 const map_type& rowMap = * (this->getRowMap ());
438 const GO gblRow = rowMap.getGlobalElement (lclRow);
439 this->getGlobalRowView (gblRow, indIn, valIn);
440 memcpy (indOut, indIn.getRawPtr (), numEnt *
sizeof (GO));
441 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
452 if (this->isLocallyIndexed ()) {
453 Array<LO> indIn (numEnt);
454 Array<Scalar> valIn (numEnt);
455 size_t theNumEnt = 0;
456 this->getLocalRowCopy (lclRow, indIn (), valIn (), theNumEnt);
457 if (theNumEnt != numEnt) {
460 const map_type& colMap = * (this->getColMap ());
463 for (
size_t k = 0; k < numEnt; ++k) {
464 const GO gblIndIn = colMap.getGlobalElement (indIn[k]);
465 memcpy (indOut + k *
sizeof (GO), &gblIndIn,
sizeof (GO));
467 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
469 else if (this->isGloballyIndexed ()) {
470 Array<GO> indIn (numEnt);
471 Array<Scalar> valIn (numEnt);
472 const map_type& rowMap = * (this->getRowMap ());
473 const GO gblRow = rowMap.getGlobalElement (lclRow);
474 size_t theNumEnt = 0;
475 this->getGlobalRowCopy (gblRow, indIn, valIn, theNumEnt);
476 if (theNumEnt != numEnt) {
479 memcpy (indOut, indIn.getRawPtr (), numEnt *
sizeof (GO));
480 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
491 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
493 RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
494 packImpl (
const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
495 Teuchos::Array<char>& exports,
496 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
497 size_t& constantNumPackets,
500 using Teuchos::Array;
501 using Teuchos::ArrayView;
503 using Teuchos::av_reinterpret_cast;
505 typedef LocalOrdinal LO;
506 typedef GlobalOrdinal GO;
507 typedef typename ArrayView<const LO>::size_type size_type;
508 const char tfecfFuncName[] =
"packImpl: ";
510 const size_type numExportLIDs = exportLIDs.size ();
511 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
512 numExportLIDs != numPacketsPerLID.size (), std::invalid_argument,
513 "exportLIDs.size() = " << numExportLIDs <<
" != numPacketsPerLID.size()"
514 " = " << numPacketsPerLID.size () <<
".");
519 constantNumPackets = 0;
524 size_t totalNumEntries = 0;
525 allocatePackSpace (exports, totalNumEntries, exportLIDs);
526 const size_t bufSize =
static_cast<size_t> (exports.size ());
538 size_type firstBadIndex = 0;
539 size_t firstBadOffset = 0;
540 size_t firstBadNumBytes = 0;
542 bool packErr =
false;
544 char*
const exportsRawPtr = exports.getRawPtr ();
546 for (size_type i = 0; i < numExportLIDs; ++i) {
547 const LO lclRow = exportLIDs[i];
548 const size_t numEnt = this->getNumEntriesInLocalRow (lclRow);
552 numPacketsPerLID[i] = 0;
555 char*
const numEntBeg = exportsRawPtr + offset;
556 char*
const numEntEnd = numEntBeg +
sizeof (LO);
557 char*
const valBeg = numEntEnd;
558 char*
const valEnd = valBeg + numEnt *
sizeof (Scalar);
559 char*
const indBeg = valEnd;
560 const size_t numBytes =
sizeof (LO) +
561 numEnt * (
sizeof (Scalar) +
sizeof (GO));
562 if (offset > bufSize || offset + numBytes > bufSize) {
564 firstBadOffset = offset;
565 firstBadNumBytes = numBytes;
569 packErr = ! packRow (numEntBeg, valBeg, indBeg, numEnt, lclRow);
572 firstBadOffset = offset;
573 firstBadNumBytes = numBytes;
579 numPacketsPerLID[i] = numBytes;
590 TEUCHOS_TEST_FOR_EXCEPTION(
591 outOfBounds, std::logic_error,
"First invalid offset into 'exports' "
592 "pack buffer at index i = " << firstBadIndex <<
". exportLIDs[i]: "
593 << exportLIDs[firstBadIndex] <<
", bufSize: " << bufSize <<
", offset: "
594 << firstBadOffset <<
", numBytes: " << firstBadNumBytes <<
".");
595 TEUCHOS_TEST_FOR_EXCEPTION(
596 packErr, std::logic_error,
"First error in packRow() at index i = "
597 << firstBadIndex <<
". exportLIDs[i]: " << exportLIDs[firstBadIndex]
598 <<
", bufSize: " << bufSize <<
", offset: " << firstBadOffset
599 <<
", numBytes: " << firstBadNumBytes <<
".");
602 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
606 LocalOrdinal& numEnt,
607 const LocalOrdinal*& lclColInds,
608 const Scalar*& vals)
const
613 Teuchos::ArrayView<const LocalOrdinal> lclColInds_av;
614 Teuchos::ArrayView<const Scalar> vals_av;
616 this->getLocalRowView (lclRow, lclColInds_av, vals_av);
617 numEnt =
static_cast<LocalOrdinal
> (lclColInds_av.size ());
623 lclColInds = lclColInds_av.getRawPtr ();
624 vals = vals_av.getRawPtr ();
627 return static_cast<LocalOrdinal
> (0);
638 #define TPETRA_ROWMATRIX_INSTANT(SCALAR,LO,GO,NODE) \
639 template class RowMatrix< SCALAR , LO , GO , NODE >;
642 #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.