42 #ifndef TPETRA_ROWMATRIX_DEF_HPP 43 #define TPETRA_ROWMATRIX_DEF_HPP 45 #include <Tpetra_ConfigDefs.hpp> 46 #include <Tpetra_CrsMatrix.hpp> 47 #include <Tpetra_Map.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::ArrayRCP;
66 using Teuchos::ArrayView;
67 using Teuchos::ParameterList;
70 using Teuchos::rcp_implicit_cast;
71 using Teuchos::sublist;
72 typedef LocalOrdinal LO;
73 typedef GlobalOrdinal GO;
74 typedef Teuchos::ScalarTraits<Scalar> STS;
79 const this_type& B = *
this;
87 RCP<const map_type> B_domainMap = B.getDomainMap ();
88 RCP<const map_type> B_rangeMap = B.getRangeMap ();
90 RCP<const map_type> theDomainMap = domainMap;
91 RCP<const map_type> theRangeMap = rangeMap;
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;
101 theDomainMap = B_domainMap;
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;
112 theRangeMap = B_rangeMap;
116 #ifdef HAVE_TPETRA_DEBUG 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.");
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.");
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 " 161 #endif // HAVE_TPETRA_DEBUG 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");
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;
178 RCP<crs_matrix_type> C;
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);
190 if (alpha != STS::zero ()) {
191 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
193 C_maxNumEntriesPerRow[localRow] += A_numEntries;
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;
204 if (constructorSublist.is_null ()) {
205 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow,
208 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow,
219 if (constructorSublist.is_null ()) {
223 constructorSublist));
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 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) {
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);
247 ArrayView<GO> indView = ind (0, A_numEntries);
248 ArrayView<Scalar> valView = val (0, A_numEntries);
251 if (alpha != STS::one ()) {
252 for (
size_t k = 0; k < A_numEntries; ++k) {
256 C->insertGlobalValues (globalRow, indView, valView);
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);
269 ArrayView<GO> indView = ind (0, B_numEntries);
270 ArrayView<Scalar> valView = val (0, B_numEntries);
271 B.getGlobalRowCopy (globalRow, indView, valView, B_numEntries);
273 if (beta != STS::one ()) {
274 for (
size_t k = 0; k < B_numEntries; ++k) {
278 C->insertGlobalValues (globalRow, indView, valView);
282 if (callFillComplete) {
283 if (fillCompleteSublist.is_null ()) {
284 C->fillComplete (theDomainMap, theRangeMap);
286 C->fillComplete (theDomainMap, theRangeMap, fillCompleteSublist);
290 return rcp_implicit_cast<this_type> (C);
294 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
297 pack (
const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
298 Teuchos::Array<char>& exports,
299 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
300 size_t& constantNumPackets,
303 #ifdef HAVE_TPETRA_DEBUG 304 const char tfecfFuncName[] =
"pack: ";
306 using Teuchos::reduceAll;
307 std::ostringstream msg;
310 this->packImpl (exportLIDs, exports, numPacketsPerLID,
311 constantNumPackets, distor);
312 }
catch (std::exception& e) {
317 const Teuchos::Comm<int>& comm = * (this->getComm ());
318 reduceAll<int, int> (comm, Teuchos::REDUCE_MAX,
319 lclBad, Teuchos::outArg (gblBad));
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 ();
333 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
334 true, std::logic_error,
"packImpl() threw an exception on one or " 335 "more participating processes.");
339 this->packImpl (exportLIDs, exports, numPacketsPerLID,
340 constantNumPackets, distor);
341 #endif // HAVE_TPETRA_DEBUG 344 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
348 size_t& totalNumEntries,
349 const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs)
const 351 typedef LocalOrdinal LO;
352 typedef GlobalOrdinal GO;
353 typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
355 const size_type numExportLIDs = exportLIDs.size ();
359 for (size_type i = 0; i < numExportLIDs; ++i) {
360 const LO lclRow = exportLIDs[i];
361 size_t curNumEntries = this->getNumEntriesInLocalRow (lclRow);
364 if (curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid ()) {
367 totalNumEntries += curNumEntries;
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);
386 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
388 RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
389 packRow (
char*
const numEntOut,
393 const LocalOrdinal lclRow)
const 395 using Teuchos::Array;
396 using Teuchos::ArrayView;
397 typedef LocalOrdinal LO;
398 typedef GlobalOrdinal GO;
401 const LO numEntLO =
static_cast<LO
> (numEnt);
402 memcpy (numEntOut, &numEntLO,
sizeof (LO));
404 if (this->supportsRowViews ()) {
405 if (this->isLocallyIndexed ()) {
409 ArrayView<const LO> indIn;
410 ArrayView<const Scalar> valIn;
411 this->getLocalRowView (lclRow, indIn, valIn);
412 const map_type& colMap = * (this->getColMap ());
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));
419 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
421 else if (this->isGloballyIndexed ()) {
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));
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) {
452 const map_type& colMap = * (this->getColMap ());
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));
459 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
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) {
471 memcpy (indOut, indIn.getRawPtr (), numEnt *
sizeof (GO));
472 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
483 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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 492 using Teuchos::Array;
493 using Teuchos::ArrayView;
495 using Teuchos::av_reinterpret_cast;
497 typedef LocalOrdinal LO;
498 typedef GlobalOrdinal GO;
499 typedef typename ArrayView<const LO>::size_type size_type;
500 const char tfecfFuncName[] =
"packImpl: ";
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 () <<
".");
511 constantNumPackets = 0;
516 size_t totalNumEntries = 0;
517 allocatePackSpace (exports, totalNumEntries, exportLIDs);
518 const size_t bufSize =
static_cast<size_t> (exports.size ());
530 size_type firstBadIndex = 0;
531 size_t firstBadOffset = 0;
532 size_t firstBadNumBytes = 0;
534 bool packErr =
false;
536 char*
const exportsRawPtr = exports.getRawPtr ();
538 for (size_type i = 0; i < numExportLIDs; ++i) {
539 const LO lclRow = exportLIDs[i];
540 const size_t numEnt = this->getNumEntriesInLocalRow (lclRow);
544 numPacketsPerLID[i] = 0;
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) {
556 firstBadOffset = offset;
557 firstBadNumBytes = numBytes;
561 packErr = ! packRow (numEntBeg, valBeg, indBeg, numEnt, lclRow);
564 firstBadOffset = offset;
565 firstBadNumBytes = numBytes;
571 numPacketsPerLID[i] = numBytes;
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 <<
".");
594 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
598 LocalOrdinal& numEnt,
599 const LocalOrdinal*& lclColInds,
600 const Scalar*& vals)
const 605 Teuchos::ArrayView<const LocalOrdinal> lclColInds_av;
606 Teuchos::ArrayView<const Scalar> vals_av;
608 this->getLocalRowView (lclRow, lclColInds_av, vals_av);
609 numEnt =
static_cast<LocalOrdinal
> (lclColInds_av.size ());
615 lclColInds = lclColInds_av.getRawPtr ();
616 vals = vals_av.getRawPtr ();
619 return static_cast<LocalOrdinal
> (0);
630 #define TPETRA_ROWMATRIX_INSTANT(SCALAR,LO,GO,NODE) \ 632 template class RowMatrix< SCALAR , LO , GO , NODE >; 635 #endif // TPETRA_ROWMATRIX_DEF_HPP 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 > ¶ms=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'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'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.