10 #ifndef TPETRA_DETAILS_COOMATRIX_HPP
11 #define TPETRA_DETAILS_COOMATRIX_HPP
18 #include "TpetraCore_config.h"
19 #include "Tpetra_Details_PackTriples.hpp"
20 #include "Tpetra_DistObject.hpp"
23 #include "Teuchos_TypeNameTraits.hpp"
25 #include <initializer_list>
29 #include <type_traits>
41 template <
class IndexType>
51 template <
class IndexType>
56 return (x.row < y.row) || (x.row == y.row && x.col < y.col);
62 template <
class SC,
class GO>
72 typedef std::map<CooGraphEntry<GO>, SC,
78 entries_type entries_;
101 auto result = this->entries_.insert({ij, val});
102 if (!result.second) {
103 result.first->second += val;
117 const GO gblColInds[],
119 const std::size_t numEnt) {
120 for (std::size_t k = 0; k < numEnt; ++k) {
125 const SC val = vals[k];
126 auto result = this->entries_.insert({ij, val});
127 if (!result.second) {
128 result.first->second += val;
136 return this->entries_.size();
144 for (
auto iter = this->entries_.begin();
145 iter != this->entries_.end(); ++iter) {
146 f(iter->first.row, iter->first.col, iter->second);
158 const GO srcGblRow) {
162 entries_type& tgtEntries = this->entries_;
163 const entries_type& srcEntries = src.entries_;
169 auto srcBeg = srcEntries.lower_bound({srcGblRow, std::numeric_limits<GO>::min()});
170 auto srcEnd = srcEntries.upper_bound({srcGblRow + 1, std::numeric_limits<GO>::min()});
171 auto tgtBeg = tgtEntries.lower_bound({tgtGblRow, std::numeric_limits<GO>::min()});
172 auto tgtEnd = tgtEntries.upper_bound({tgtGblRow + 1, std::numeric_limits<GO>::min()});
176 if (srcBeg != srcEnd) {
177 for (
auto tgtCur = tgtBeg; tgtCur != tgtEnd; ++tgtCur) {
178 auto srcCur = srcBeg;
179 while (srcCur != srcEnd && srcCur->first.col < tgtCur->first.col) {
191 if (srcCur != srcEnd) {
192 if (srcCur->first.col == tgtCur->first.col) {
193 tgtCur->second += srcCur->second;
199 (void)tgtEntries.insert(tgtCur, *srcCur);
217 const ::Teuchos::Comm<int>& comm,
218 std::ostream* errStrm = NULL)
const {
220 using ::Tpetra::Details::countPackTriples;
221 using ::Tpetra::Details::countPackTriplesCount;
222 const char prefix[] =
"Tpetra::Details::Impl::CooMatrixImpl::countPackRow: ";
223 #ifdef HAVE_TPETRACORE_MPI
224 int errCode = MPI_SUCCESS;
227 #endif // HAVE_TPETRACORE_MPI
230 const GO minGO = std::numeric_limits<GO>::min();
231 auto beg = this->entries_.lower_bound({gblRow, minGO});
232 auto end = this->entries_.upper_bound({gblRow + 1, minGO});
234 for (
auto iter = beg; iter != end; ++iter) {
236 if (numEnt == std::numeric_limits<int>::max()) {
237 if (errStrm != NULL) {
238 *errStrm << prefix <<
"In (global) row " << gblRow <<
", the number "
239 "of entries thus far, numEnt = "
240 << numEnt <<
", has reached the "
241 "maximum int value. We need to store this count as int for MPI's "
249 int numPacketsForCount = 0;
254 if (errStrm != NULL) {
255 *errStrm << prefix <<
"countPackTriplesCount "
256 "returned errCode = "
257 << errCode <<
" != 0." << endl;
261 if (numPacketsForCount < 0) {
262 if (errStrm != NULL) {
263 *errStrm << prefix <<
"countPackTriplesCount returned "
264 "numPacketsForCount = "
265 << numPacketsForCount <<
" < 0." << endl;
271 int numPacketsForTriples = 0;
273 errCode = countPackTriples<SC, GO>(numEnt, comm, numPacketsForTriples);
274 TEUCHOS_TEST_FOR_EXCEPTION(errCode != 0, std::runtime_error, prefix <<
"countPackTriples "
275 "returned errCode = "
276 << errCode <<
" != 0.");
277 TEUCHOS_TEST_FOR_EXCEPTION(numPacketsForTriples < 0, std::logic_error, prefix <<
"countPackTriples "
278 "returned numPacketsForTriples = "
279 << numPacketsForTriples <<
" < 0.");
282 numPackets = numPacketsForCount + numPacketsForTriples;
302 const int outBufSize,
304 const ::Teuchos::Comm<int>& comm,
305 std::vector<GO>& gblRowInds,
306 std::vector<GO>& gblColInds,
307 std::vector<SC>& vals,
308 const GO gblRow)
const {
309 using ::Tpetra::Details::packTriples;
310 using ::Tpetra::Details::packTriplesCount;
311 const char prefix[] =
"Tpetra::Details::Impl::CooMatrixImpl::packRow: ";
313 const GO minGO = std::numeric_limits<GO>::min();
314 auto beg = this->entries_.lower_bound({gblRow, minGO});
315 auto end = this->entries_.upper_bound({gblRow + 1, minGO});
319 gblRowInds.resize(0);
320 gblColInds.resize(0);
324 for (
auto iter = beg; iter != end; ++iter) {
325 gblRowInds.push_back(iter->first.row);
326 gblColInds.push_back(iter->first.col);
327 vals.push_back(iter->second);
329 TEUCHOS_TEST_FOR_EXCEPTION(numEnt == std::numeric_limits<int>::max(), std::runtime_error, prefix <<
"In (global) row " << gblRow <<
", the number of entries thus far, "
331 << numEnt <<
", has reached the maximum int value. "
332 "We need to store this count as int for MPI's sake.");
338 TEUCHOS_TEST_FOR_EXCEPTION(errCode != 0, std::runtime_error, prefix <<
"In (global) row " << gblRow <<
", packTriplesCount returned "
340 << errCode <<
" != 0.");
352 TEUCHOS_TEST_FOR_EXCEPTION(errCode != 0, std::runtime_error, prefix <<
"In (global) row " << gblRow <<
", packTriples returned errCode = " << errCode <<
" != 0.");
366 GO lclMinRowInd = std::numeric_limits<GO>::max();
367 for (
typename entries_type::const_iterator iter = this->entries_.begin();
368 iter != this->entries_.end(); ++iter) {
369 const GO rowInd = iter->first.row;
370 if (rowInd < lclMinRowInd) {
371 lclMinRowInd = rowInd;
373 if (rowInds.empty() || rowInds.back() != rowInd) {
374 rowInds.push_back(rowInd);
380 template <
class OffsetType>
382 buildCrs(std::vector<OffsetType>& rowOffsets,
385 static_assert(std::is_integral<OffsetType>::value,
386 "OffsetType must be a built-in integer type.");
394 rowOffsets.push_back(0);
396 typename entries_type::const_iterator iter = this->entries_.begin();
397 GO prevGblRowInd = iter->first.row;
400 for (; iter != this->entries_.end(); ++iter, ++k) {
401 const GO gblRowInd = iter->first.row;
402 if (k == 0 || gblRowInd != prevGblRowInd) {
406 rowOffsets.push_back(k);
407 prevGblRowInd = gblRowInd;
409 gblColInds[k] = iter->first.col;
411 static_assert(std::is_same<
typename std::decay<decltype(iter->second)>::type, SC>::value,
412 "The type of iter->second != SC.");
413 vals[k] = iter->second;
415 rowOffsets.push_back(static_cast<OffsetType>(numEnt));
431 template <
class OffsetType,
class LO>
436 std::function<LO(
const GO)> gblToLcl)
const {
437 static_assert(std::is_integral<OffsetType>::value,
438 "OffsetType must be a built-in integer type.");
439 static_assert(std::is_integral<LO>::value,
440 "LO must be a built-in integer type.");
448 rowOffsets.push_back(0);
450 typename entries_type::const_iterator iter = this->entries_.begin();
451 GO prevGblRowInd = iter->first.row;
454 for (; iter != this->entries_.end(); ++iter, ++k) {
455 const GO gblRowInd = iter->first.row;
456 if (k == 0 || gblRowInd != prevGblRowInd) {
460 rowOffsets.push_back(k);
461 prevGblRowInd = gblRowInd;
463 lclColInds[k] = gblToLcl(iter->first.col);
464 vals[k] = iter->second;
466 rowOffsets.push_back(static_cast<OffsetType>(numEnt));
531 typedef LO local_ordinal_type;
532 typedef GO global_ordinal_type;
533 typedef NT node_type;
534 typedef typename NT::device_type device_type;
536 typedef ::Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type>
map_type;
540 typedef ::Tpetra::DistObject<
packet_type, local_ordinal_type,
541 global_ordinal_type, node_type>
555 , localError_(new bool(false))
556 , errs_(new std::shared_ptr<std::ostringstream>())
568 , localError_(new bool(false))
569 , errs_(new std::shared_ptr<std::ostringstream>())
585 this->impl_.sumIntoGlobalValue(gblRowInd, gblColInd, val);
598 const GO gblColInds[],
600 const std::size_t numEnt) {
601 this->impl_.sumIntoGlobalValues(gblRowInds, gblColInds, vals, numEnt);
607 std::initializer_list<GO> gblColInds,
608 std::initializer_list<SC> vals,
609 const std::size_t numEnt) {
610 this->impl_.sumIntoGlobalValues(gblRowInds.begin(), gblColInds.begin(),
611 vals.begin(), numEnt);
617 return this->impl_.getLclNumEntries();
620 template <
class OffsetType>
622 buildCrs(::Kokkos::View<OffsetType*, device_type>& rowOffsets,
623 ::Kokkos::View<GO*, device_type>& gblColInds,
624 #
if KOKKOS_VERSION >= 40799
625 ::Kokkos::View<typename ::KokkosKernels::ArithTraits<SC>::val_type*, device_type>& vals) {
627 ::Kokkos::View<typename ::Kokkos::ArithTraits<SC>::val_type*, device_type>& vals) {
629 static_assert(std::is_integral<OffsetType>::value,
630 "OffsetType must be a built-in integer type.");
631 using ::Kokkos::create_mirror_view;
633 using ::Kokkos::View;
634 #if KOKKOS_VERSION >= 40799
635 typedef typename ::KokkosKernels::ArithTraits<SC>::val_type ISC;
637 typedef typename ::Kokkos::ArithTraits<SC>::val_type ISC;
642 gblColInds = View<GO*, device_type>(
"gblColInds", numEnt);
643 vals = View<ISC*, device_type>(
"vals", numEnt);
644 auto gblColInds_h = create_mirror_view(gblColInds);
645 auto vals_h = create_mirror_view(vals);
647 std::vector<std::size_t> rowOffsetsSV;
648 this->impl_.buildCrs(rowOffsetsSV,
652 View<OffsetType*, device_type>(
"rowOffsets", rowOffsetsSV.size());
653 typename View<OffsetType*, device_type>::host_mirror_type
654 rowOffsets_h(rowOffsetsSV.data(), rowOffsetsSV.size());
672 fillComplete(const ::Teuchos::RCP<const ::Teuchos::Comm<int> >& comm) {
673 if (comm.is_null()) {
674 this->
map_ = ::Teuchos::null;
676 this->
map_ = this->buildMap(comm);
690 TEUCHOS_TEST_FOR_EXCEPTION(this->
getMap().is_null(), std::runtime_error,
692 "CooMatrix::fillComplete: This object does not yet have a Map. "
693 "You must call the version of fillComplete "
694 "that takes an input communicator.");
734 return ((*errs_).get() == NULL) ? std::string(
"") : (*errs_)->str();
750 std::shared_ptr<bool> localError_;
759 std::shared_ptr<std::shared_ptr<std::ostringstream> > errs_;
763 markLocalErrorAndGetStream() {
764 *(this->localError_) =
true;
765 if ((*errs_).get() == NULL) {
766 *errs_ = std::shared_ptr<std::ostringstream>(
new std::ostringstream());
775 using Teuchos::TypeNameTraits;
777 std::ostringstream os;
778 os <<
"\"Tpetra::Details::CooMatrix\": {"
779 <<
"SC: " << TypeNameTraits<SC>::name()
780 <<
", LO: " << TypeNameTraits<LO>::name()
781 <<
", GO: " << TypeNameTraits<GO>::name()
782 <<
", NT: " << TypeNameTraits<NT>::name();
783 if (this->getObjectLabel() !=
"") {
784 os <<
", Label: \"" << this->getObjectLabel() <<
"\"";
786 os <<
", Has Map: " << (this->
map_.is_null() ?
"false" :
"true")
795 const Teuchos::EVerbosityLevel verbLevel =
796 Teuchos::Describable::verbLevel_default)
const {
798 using ::Teuchos::EVerbosityLevel;
799 using ::Teuchos::OSTab;
800 using ::Teuchos::TypeNameTraits;
801 using ::Teuchos::VERB_DEFAULT;
802 using ::Teuchos::VERB_LOW;
803 using ::Teuchos::VERB_MEDIUM;
804 using ::Tpetra::Details::gathervPrint;
806 const auto vl = (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
808 auto comm = this->
getMap().is_null() ? Teuchos::null : this->
getMap()->getComm();
809 const int myRank = comm.is_null() ? 0 : comm->getRank();
812 if (vl != Teuchos::VERB_NONE) {
816 out <<
"\"Tpetra::Details::CooMatrix\":" << endl;
820 out <<
"Template parameters:" << endl;
823 out <<
"SC: " << TypeNameTraits<SC>::name() << endl
824 <<
"LO: " << TypeNameTraits<LO>::name() << endl
825 <<
"GO: " << TypeNameTraits<GO>::name() << endl
826 <<
"NT: " << TypeNameTraits<NT>::name() << endl;
828 if (this->getObjectLabel() !=
"") {
829 out <<
"Label: \"" << this->getObjectLabel() <<
"\"" << endl;
831 out <<
"Has Map: " << (this->
map_.is_null() ?
"false" :
"true") << endl;
835 if (!this->
map_.is_null()) {
837 out <<
"Map:" << endl;
840 this->
map_->describe(out, vl);
845 std::ostringstream os;
846 os <<
"Process " << myRank <<
":" << endl;
849 const std::size_t numEnt = this->impl_.getLclNumEntries();
850 os <<
" Local number of entries: " << numEnt << endl;
851 if (vl > VERB_MEDIUM) {
852 os <<
" Local entries:" << endl;
854 this->impl_.forAllEntries([&os](
const GO row,
const GO col,
const SC& val) {
876 Teuchos::RCP<const map_type>
877 buildMap(const ::Teuchos::RCP<const ::Teuchos::Comm<int> >& comm) {
878 using ::Teuchos::outArg;
879 using ::Teuchos::rcp;
880 using ::Teuchos::REDUCE_MIN;
881 using ::Teuchos::reduceAll;
882 typedef ::Tpetra::global_size_t GST;
886 if (comm.is_null()) {
887 return ::Teuchos::null;
898 std::vector<GO> rowInds;
899 const GO lclMinGblRowInd = this->impl_.getMyGlobalRowIndices(rowInds);
902 GO gblMinGblRowInd = 0;
903 reduceAll<int, GO>(*comm, REDUCE_MIN, lclMinGblRowInd,
904 outArg(gblMinGblRowInd));
905 const GO indexBase = gblMinGblRowInd;
906 const GST INV = Tpetra::Details::OrdinalTraits<GST>::invalid();
907 return rcp(
new map_type(INV, rowInds.data(), rowInds.size(),
916 return static_cast<size_t>(0);
926 const char prefix[] =
"Tpetra::Details::CooMatrix::checkSizes: ";
928 const this_COO_type* src =
dynamic_cast<const this_COO_type*
>(&source);
931 std::ostream& err = markLocalErrorAndGetStream();
932 err << prefix <<
"The source object of the Import or Export "
933 "must be a CooMatrix with the same template parameters as the "
936 }
else if (this->
map_.is_null()) {
937 std::ostream& err = markLocalErrorAndGetStream();
938 err << prefix <<
"The target object of the Import or Export "
939 "must be a CooMatrix with a nonnull Map."
942 return !(*(this->localError_));
947 typename ::Tpetra::DistObject<char, LO, GO, NT>::buffer_device_type;
950 copyAndPermute(const ::Tpetra::SrcDistObject& sourceObject,
951 const size_t numSameIDs,
952 const Kokkos::DualView<
const LO*,
954 const Kokkos::DualView<
const LO*,
959 const char prefix[] =
"Tpetra::Details::CooMatrix::copyAndPermute: ";
964 if (*(this->localError_)) {
965 std::ostream& err = this->markLocalErrorAndGetStream();
966 err << prefix <<
"The target object of the Import or Export is "
967 "already in an error state."
972 const this_COO_type* src =
dynamic_cast<const this_COO_type*
>(&sourceObject);
973 if (src ==
nullptr) {
974 std::ostream& err = this->markLocalErrorAndGetStream();
975 err << prefix <<
"Input argument 'sourceObject' is not a CooMatrix."
980 const size_t numPermuteIDs =
981 static_cast<size_t>(permuteToLIDs.extent(0));
982 if (numPermuteIDs != static_cast<size_t>(permuteFromLIDs.extent(0))) {
983 std::ostream& err = this->markLocalErrorAndGetStream();
984 err << prefix <<
"permuteToLIDs.extent(0) = "
985 << numPermuteIDs <<
" != permuteFromLIDs.extent(0) = "
986 << permuteFromLIDs.extent(0) <<
"." << endl;
989 if (
sizeof(
int) <=
sizeof(
size_t) &&
990 numPermuteIDs > static_cast<size_t>(std::numeric_limits<int>::max())) {
991 std::ostream& err = this->markLocalErrorAndGetStream();
992 err << prefix <<
"numPermuteIDs = " << numPermuteIDs
993 <<
", a size_t, overflows int." << endl;
1000 if (
sizeof(
int) <=
sizeof(
size_t) &&
1001 numSameIDs > static_cast<size_t>(std::numeric_limits<int>::max())) {
1002 std::ostream& err = this->markLocalErrorAndGetStream();
1003 err << prefix <<
"numSameIDs = " << numSameIDs
1004 <<
", a size_t, overflows int." << endl;
1011 const LO numSame =
static_cast<int>(numSameIDs);
1016 LO numInvalidSameRows = 0;
1017 for (LO lclRow = 0; lclRow < numSame; ++lclRow) {
1021 const GO gblRow = this->
map_->getGlobalElement(lclRow);
1022 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1023 ++numInvalidSameRows;
1026 this->impl_.mergeIntoRow(gblRow, src->impl_, gblRow);
1035 const LO numPermute =
static_cast<int>(numPermuteIDs);
1040 LO numInvalidRowsFrom = 0;
1045 LO numInvalidRowsTo = 0;
1047 TEUCHOS_ASSERT(!permuteFromLIDs.need_sync_host());
1048 TEUCHOS_ASSERT(!permuteToLIDs.need_sync_host());
1049 auto permuteFromLIDs_h = permuteFromLIDs.view_host();
1050 auto permuteToLIDs_h = permuteToLIDs.view_host();
1052 for (LO k = 0; k < numPermute; ++k) {
1053 const LO lclRowFrom = permuteFromLIDs_h[k];
1054 const LO lclRowTo = permuteToLIDs_h[k];
1055 const GO gblRowFrom = src->map_->getGlobalElement(lclRowFrom);
1056 const GO gblRowTo = this->
map_->getGlobalElement(lclRowTo);
1058 bool bothConversionsValid =
true;
1059 if (gblRowFrom == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1060 ++numInvalidRowsFrom;
1061 bothConversionsValid =
false;
1063 if (gblRowTo == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1065 bothConversionsValid =
false;
1067 if (bothConversionsValid) {
1068 this->impl_.mergeIntoRow(gblRowTo, src->impl_, gblRowFrom);
1073 if (numInvalidSameRows != 0 || numInvalidRowsFrom != 0 ||
1074 numInvalidRowsTo != 0) {
1077 std::vector<std::pair<LO, GO> > invalidSameRows;
1078 invalidSameRows.reserve(numInvalidSameRows);
1079 std::vector<std::pair<LO, GO> > invalidRowsFrom;
1080 invalidRowsFrom.reserve(numInvalidRowsFrom);
1081 std::vector<std::pair<LO, GO> > invalidRowsTo;
1082 invalidRowsTo.reserve(numInvalidRowsTo);
1084 for (LO lclRow = 0; lclRow < numSame; ++lclRow) {
1088 const GO gblRow = this->
map_->getGlobalElement(lclRow);
1089 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1090 invalidSameRows.push_back({lclRow, gblRow});
1094 for (LO k = 0; k < numPermute; ++k) {
1095 const LO lclRowFrom = permuteFromLIDs_h[k];
1096 const LO lclRowTo = permuteToLIDs_h[k];
1097 const GO gblRowFrom = src->map_->getGlobalElement(lclRowFrom);
1098 const GO gblRowTo = this->
map_->getGlobalElement(lclRowTo);
1100 if (gblRowFrom == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1101 invalidRowsFrom.push_back({lclRowFrom, gblRowFrom});
1103 if (gblRowTo == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1104 invalidRowsTo.push_back({lclRowTo, gblRowTo});
1108 std::ostringstream os;
1109 if (numInvalidSameRows != 0) {
1110 os <<
"Invalid permute \"same\" (local, global) index pairs: [";
1111 for (std::size_t k = 0; k < invalidSameRows.size(); ++k) {
1112 const auto& p = invalidSameRows[k];
1113 os <<
"(" << p.first <<
"," << p.second <<
")";
1114 if (k + 1 < invalidSameRows.size()) {
1118 os <<
"]" << std::endl;
1120 if (numInvalidRowsFrom != 0) {
1121 os <<
"Invalid permute \"from\" (local, global) index pairs: [";
1122 for (std::size_t k = 0; k < invalidRowsFrom.size(); ++k) {
1123 const auto& p = invalidRowsFrom[k];
1124 os <<
"(" << p.first <<
"," << p.second <<
")";
1125 if (k + 1 < invalidRowsFrom.size()) {
1129 os <<
"]" << std::endl;
1131 if (numInvalidRowsTo != 0) {
1132 os <<
"Invalid permute \"to\" (local, global) index pairs: [";
1133 for (std::size_t k = 0; k < invalidRowsTo.size(); ++k) {
1134 const auto& p = invalidRowsTo[k];
1135 os <<
"(" << p.first <<
"," << p.second <<
")";
1136 if (k + 1 < invalidRowsTo.size()) {
1140 os <<
"]" << std::endl;
1143 std::ostream& err = this->markLocalErrorAndGetStream();
1144 err << prefix << os.str();
1150 packAndPrepare(const ::Tpetra::SrcDistObject& sourceObject,
1151 const Kokkos::DualView<
const local_ordinal_type*,
1155 Kokkos::DualView<
size_t*,
1158 size_t& constantNumPackets) {
1160 using Teuchos::Comm;
1163 const char prefix[] =
"Tpetra::Details::CooMatrix::packAndPrepare: ";
1164 const char suffix[] =
1165 " This should never happen. "
1166 "Please report this bug to the Tpetra developers.";
1170 constantNumPackets = 0;
1172 const this_COO_type* src =
dynamic_cast<const this_COO_type*
>(&sourceObject);
1173 if (src ==
nullptr) {
1174 std::ostream& err = this->markLocalErrorAndGetStream();
1175 err << prefix <<
"Input argument 'sourceObject' is not a CooMatrix."
1177 }
else if (*(src->localError_)) {
1178 std::ostream& err = this->markLocalErrorAndGetStream();
1179 err << prefix <<
"The source (input) object of the Import or Export "
1180 "is already in an error state on this process."
1182 }
else if (*(this->localError_)) {
1183 std::ostream& err = this->markLocalErrorAndGetStream();
1184 err << prefix <<
"The target (output, \"this\") object of the Import "
1185 "or Export is already in an error state on this process."
1191 if (*(this->localError_)) {
1196 auto numPacketsPerLID_tmp = numPacketsPerLID;
1197 numPacketsPerLID_tmp.sync_host();
1198 numPacketsPerLID_tmp.modify_host();
1205 const size_t numExports = exportLIDs.extent(0);
1206 if (numExports == 0) {
1210 RCP<const Comm<int> > comm = src->getMap().is_null() ? Teuchos::null : src->getMap()->getComm();
1211 if (comm.is_null() || comm->getSize() == 1) {
1212 if (numExports != static_cast<size_t>(0)) {
1213 std::ostream& err = this->markLocalErrorAndGetStream();
1214 err << prefix <<
"The input communicator is either null or "
1215 "has only one process, but numExports = "
1216 << numExports <<
" != 0."
1227 numPacketsPerLID.sync_host();
1228 numPacketsPerLID.modify_host();
1230 TEUCHOS_ASSERT(!exportLIDs.need_sync_host());
1231 auto exportLIDs_h = exportLIDs.view_host();
1233 int totalNumPackets = 0;
1234 size_t numInvalidRowInds = 0;
1235 std::ostringstream errStrm;
1236 for (
size_t k = 0; k < numExports; ++k) {
1237 const LO lclRow = exportLIDs_h[k];
1240 const GO gblRow = src->map_->getGlobalElement(lclRow);
1241 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1243 ++numInvalidRowInds;
1244 numPacketsPerLID.view_host()[k] = 0;
1252 src->impl_.countPackRow(numPackets, gblRow, *comm, &errStrm);
1254 std::ostream& err = this->markLocalErrorAndGetStream();
1255 err << prefix << errStrm.str() << endl;
1256 numPacketsPerLID.view_host()[k] = 0;
1262 const long long newTotalNumPackets =
1263 static_cast<long long>(totalNumPackets) +
1264 static_cast<long long>(numPackets);
1265 if (newTotalNumPackets >
1266 static_cast<long long>(std::numeric_limits<int>::max())) {
1267 std::ostream& err = this->markLocalErrorAndGetStream();
1268 err << prefix <<
"The new total number of packets "
1269 << newTotalNumPackets <<
" does not fit in int." << endl;
1273 for (
size_t k2 = k; k2 < numExports; ++k2) {
1274 numPacketsPerLID.view_host()[k2] = 0;
1278 numPacketsPerLID.view_host()[k] =
static_cast<size_t>(numPackets);
1279 totalNumPackets =
static_cast<int>(newTotalNumPackets);
1284 if (numInvalidRowInds != 0) {
1285 std::vector<std::pair<LO, GO> > invalidRowInds;
1286 for (
size_t k = 0; k < numExports; ++k) {
1287 const LO lclRow = exportLIDs_h[k];
1291 const GO gblRow = src->map_->getGlobalElement(lclRow);
1292 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1293 invalidRowInds.push_back({lclRow, gblRow});
1296 std::ostringstream os;
1297 os << prefix <<
"We found " << numInvalidRowInds <<
" invalid row ind"
1298 << (numInvalidRowInds !=
static_cast<size_t>(1) ?
"ices" :
"ex")
1299 <<
" out of " << numExports <<
" in exportLIDs. Here is the list "
1300 <<
" of invalid row indices: [";
1301 for (
size_t k = 0; k < invalidRowInds.size(); ++k) {
1302 os <<
"(LID: " << invalidRowInds[k].first <<
", GID: "
1303 << invalidRowInds[k].second <<
")";
1304 if (k + 1 < invalidRowInds.size()) {
1310 std::ostream& err = this->markLocalErrorAndGetStream();
1311 err << prefix << os.str() << std::endl;
1316 const bool reallocated =
1318 "CooMatrix exports");
1320 exports.sync_host();
1323 exports.modify_host();
1327 std::vector<GO> gblRowInds;
1328 std::vector<GO> gblColInds;
1329 std::vector<SC> vals;
1331 int outBufCurPos = 0;
1333 for (
size_t k = 0; k < numExports; ++k) {
1334 const LO lclRow = exportLIDs.view_host()[k];
1337 const GO gblRow = src->map_->getGlobalElement(lclRow);
1339 src->impl_.packRow(outBuf, totalNumPackets, outBufCurPos, *comm,
1340 gblRowInds, gblColInds, vals, gblRow);
1345 unpackAndCombine(
const Kokkos::DualView<
const local_ordinal_type*,
1350 Kokkos::DualView<
size_t*,
1356 using Teuchos::Comm;
1358 const char prefix[] =
"Tpetra::Details::CooMatrix::unpackAndCombine: ";
1359 const char suffix[] =
1360 " This should never happen. "
1361 "Please report this bug to the Tpetra developers.";
1363 TEUCHOS_ASSERT(!importLIDs.need_sync_host());
1364 auto importLIDs_h = importLIDs.view_host();
1366 const std::size_t numImports = importLIDs.extent(0);
1367 if (numImports == 0) {
1369 }
else if (imports.extent(0) == 0) {
1370 std::ostream& err = this->markLocalErrorAndGetStream();
1371 err << prefix <<
"importLIDs.extent(0) = " << numImports <<
" != 0, "
1372 <<
"but imports.extent(0) = 0. This doesn't make sense, because "
1373 <<
"for every incoming LID, CooMatrix packs at least the count of "
1374 <<
"triples associated with that LID, even if the count is zero. "
1375 <<
"importLIDs = [";
1376 for (std::size_t k = 0; k < numImports; ++k) {
1377 err << importLIDs_h[k];
1378 if (k + 1 < numImports) {
1382 err <<
"]. " << suffix << endl;
1386 RCP<const Comm<int> > comm = this->
getMap().is_null() ? Teuchos::null : this->
getMap()->getComm();
1387 if (comm.is_null() || comm->getSize() == 1) {
1388 if (numImports != static_cast<size_t>(0)) {
1389 std::ostream& err = this->markLocalErrorAndGetStream();
1390 err << prefix <<
"The input communicator is either null or "
1391 "has only one process, but numImports = "
1392 << numImports <<
" != 0."
1405 if (static_cast<size_t>(imports.extent(0)) >
1406 static_cast<size_t>(std::numeric_limits<int>::max())) {
1407 std::ostream& err = this->markLocalErrorAndGetStream();
1408 err << prefix <<
"imports.extent(0) = "
1409 << imports.extent(0) <<
" does not fit in int." << endl;
1412 const int inBufSize =
static_cast<int>(imports.extent(0));
1414 if (imports.need_sync_host()) {
1415 imports.sync_host();
1417 if (numPacketsPerLID.need_sync_host()) {
1418 numPacketsPerLID.sync_host();
1420 auto imports_h = imports.view_host();
1421 auto numPacketsPerLID_h = numPacketsPerLID.view_host();
1425 std::vector<GO> gblRowInds;
1426 std::vector<GO> gblColInds;
1427 std::vector<SC> vals;
1430 int inBufCurPos = 0;
1431 size_t numInvalidRowInds = 0;
1433 std::ostringstream errStrm;
1434 for (
size_t k = 0; k < numImports; ++k) {
1435 const LO lclRow = importLIDs_h(k);
1436 const GO gblRow = this->
map_->getGlobalElement(lclRow);
1437 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1438 ++numInvalidRowInds;
1444 const int origInBufCurPos = inBufCurPos;
1448 numEnt, *comm, &errStrm);
1449 if (errCode != 0 || numEnt < 0 || inBufCurPos < origInBufCurPos) {
1450 std::ostream& err = this->markLocalErrorAndGetStream();
1452 err << prefix <<
"In unpack loop, k=" << k <<
": ";
1454 err <<
" unpackTriplesCount returned errCode = " << errCode
1455 <<
" != 0." << endl;
1458 err <<
" unpackTriplesCount returned errCode = 0, but numEnt = "
1459 << numEnt <<
" < 0." << endl;
1461 if (inBufCurPos < origInBufCurPos) {
1462 err <<
" After unpackTriplesCount, inBufCurPos = " << inBufCurPos
1463 <<
" < origInBufCurPos = " << origInBufCurPos <<
"." << endl;
1465 err <<
" unpackTriplesCount report: " << errStrm.str() << endl;
1466 err << suffix << endl;
1473 #ifdef HAVE_TPETRA_DEBUG
1480 #endif // HAVE_TPETRA_DEBUG
1485 gblRowInds.resize(numEnt);
1486 gblColInds.resize(numEnt);
1487 vals.resize(numEnt);
1490 gblRowInds.data(), gblColInds.data(),
1491 vals.data(), numEnt, *comm, &errStrm);
1493 std::ostream& err = this->markLocalErrorAndGetStream();
1494 err << prefix <<
"unpackTriples returned errCode = "
1495 << errCode <<
" != 0. It reports: " << errStrm.str()
1502 #ifdef HAVE_TPETRA_DEBUG
1509 #endif // HAVE_TPETRA_DEBUG
1512 vals.data(), numEnt);
1517 if (numInvalidRowInds != 0) {
1521 std::ostream& err = this->markLocalErrorAndGetStream();
1523 std::vector<std::pair<LO, GO> > invalidRowInds;
1524 for (
size_t k = 0; k < numImports; ++k) {
1525 const LO lclRow = importLIDs_h(k);
1526 const GO gblRow = this->
map_->getGlobalElement(lclRow);
1527 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1528 invalidRowInds.push_back({lclRow, gblRow});
1532 err << prefix <<
"We found " << numInvalidRowInds <<
" invalid row ind"
1533 << (numInvalidRowInds !=
static_cast<size_t>(1) ?
"ices" :
"ex")
1534 <<
" out of " << numImports <<
" in importLIDs. Here is the list "
1535 <<
" of invalid row indices: [";
1536 for (
size_t k = 0; k < invalidRowInds.size(); ++k) {
1537 err <<
"(LID: " << invalidRowInds[k].first <<
", GID: "
1538 << invalidRowInds[k].second <<
")";
1539 if (k + 1 < invalidRowInds.size()) {
1552 #endif // TPETRA_DETAILS_COOMATRIX_HPP
char packet_type
Type for packing and unpacking data.
int packTriplesCount(const int, char[], const int, int &, const ::Teuchos::Comm< int > &, std::ostream *errStrm)
Pack the count (number) of matrix triples.
void sumIntoGlobalValues(std::initializer_list< GO > gblRowInds, std::initializer_list< GO > gblColInds, std::initializer_list< SC > vals, const std::size_t numEnt)
Initializer-list overload of the above method (which see).
Node node_type
The Node type. If you don't know what this is, don't use it.
CooMatrixImpl()=default
Default constructor.
std::string errorMessages() const
The current stream of error messages.
Declaration of a function that prints strings from each process.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print a descriptiion of this object to the given output stream; overrides Teuchos::Describable method...
void packRow(packet_type outBuf[], const int outBufSize, int &outBufCurPos, const ::Teuchos::Comm< int > &comm, std::vector< GO > &gblRowInds, std::vector< GO > &gblColInds, std::vector< SC > &vals, const GO gblRow) const
Pack the given row of the matrix.
Function comparing two CooGraphEntry structs, lexicographically, first by row index, then by column index.
void sumIntoGlobalValue(const GO gblRowInd, const GO gblColInd, const SC &val)
Insert one entry locally into the sparse matrix, if it does not exist there yet. If it does exist...
Sparse matrix used only for file input / output.
void buildLocallyIndexedCrs(std::vector< OffsetType > &rowOffsets, LO lclColInds[], SC vals[], std::function< LO(const GO)> gblToLcl) const
Build a locally indexed version of CRS storage.
void fillComplete()
Special version of fillComplete that assumes that the matrix already has a Map, and reuses its commun...
GlobalOrdinal global_ordinal_type
The type of global indices.
virtual ~CooMatrix()
Destructor (virtual for memory safety of derived classes).
void gathervPrint(std::ostream &out, const std::string &s, const Teuchos::Comm< int > &comm)
On Process 0 in the given communicator, print strings from each process in that communicator, in rank order.
void fillComplete(const ::Teuchos::RCP< const ::Teuchos::Comm< int > > &comm)
Tell the matrix that you are done inserting entries locally, and that the matrix should build its Map...
void deep_copy(MultiVector< DS, DL, DG, DN > &dst, const MultiVector< SS, SL, SG, SN > &src)
Copy the contents of the MultiVector src into dst.
std::size_t getLclNumEntries() const
Number of entries in the sparse matrix on the calling process.
void sumIntoGlobalValue(const GO gblRowInd, const GO gblColInd, const SC &val)
Insert one entry locally into the sparse matrix, if it does not exist there yet. If it does exist...
Type of each (row index, column index) pair in the Tpetra::Details::CooMatrix (see below)...
bool localError() const
Whether this object had an error on the calling process.
int unpackTriplesCount(const char[], const int, int &, int &, const ::Teuchos::Comm< int > &, std::ostream *errStrm)
Unpack just the count of triples from the given input buffer.
SC scalar_type
Type of each entry (value) in the sparse matrix.
CooMatrix()
Default constructor.
virtual size_t constantNumberOfPackets() const
By returning 0, tell DistObject that this class may not necessarily have a constant number of "packet...
void sumIntoGlobalValues(const GO gblRowInds[], const GO gblColInds[], const SC vals[], const std::size_t numEnt)
Insert multiple entries locally into the sparse matrix.
CooMatrix(const ::Teuchos::RCP< const map_type > &map)
Constructor that takes a Map.
virtual bool checkSizes(const ::Tpetra::SrcDistObject &source)
Compare the source and target (this) objects for compatibility.
void sumIntoGlobalValues(const GO gblRowInds[], const GO gblColInds[], const SC vals[], const std::size_t numEnt)
Insert multiple entries locally into the sparse matrix.
CombineMode
Rule for combining data in an Import or Export.
Teuchos::RCP< const map_type > map_
The Map over which this object is distributed.
bool reallocDualViewIfNeeded(Kokkos::DualView< ValueType *, DeviceType > &dv, const size_t newSize, const char newLabel[], const size_t tooBigFactor=2, const bool needFenceBeforeRealloc=true)
Reallocate the DualView in/out argument, if needed.
int countPackTriplesCount(const ::Teuchos::Comm< int > &, int &size, std::ostream *errStrm)
Compute the buffer size required by packTriples for packing the number of matrix entries ("triples")...
::Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Type of the Map specialization to give to the constructor.
Declaration and definition of Tpetra::Details::reallocDualViewIfNeeded, an implementation detail of T...
void forAllEntries(std::function< void(const GO, const GO, const SC &)> f) const
Execute the given function for all entries of the sparse matrix, sequentially (no thread parallelism)...
LocalOrdinal local_ordinal_type
The type of local indices.
int countPackRow(int &numPackets, const GO gblRow, const ::Teuchos::Comm< int > &comm, std::ostream *errStrm=NULL) const
Count the number of packets (bytes, in this case) needed to pack the given row of the matrix...
std::size_t getLclNumEntries() const
Number of entries in the sparse matrix on the calling process.
GO getMyGlobalRowIndices(std::vector< GO > &rowInds) const
Get the global row indices on this process, sorted and made unique, and return the minimum global row...
typename::Tpetra::DistObject< char, LO, GO, NT >::buffer_device_type buffer_device_type
Kokkos::Device specialization for DistObject communication buffers.
void mergeIntoRow(const GO tgtGblRow, const CooMatrixImpl< SC, GO > &src, const GO srcGblRow)
Into global row tgtGblRow of *this, merge global row srcGblRow of src.
char packet_type
This class transfers data as bytes (MPI_BYTE).
virtual std::string description() const
One-line descriptiion of this object; overrides Teuchos::Describable method.
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
int unpackTriples(const char[], const int, int &, OrdinalType[], OrdinalType[], ScalarType[], const int, const ::Teuchos::Comm< int > &, std::ostream *errStrm=NULL)
Unpack matrix entries ("triples" (i, j, A(i,j))) from the given input buffer.
Base class for distributed Tpetra objects that support data redistribution.
Implementation detail of Tpetra::Details::CooMatrix (which see below).
int packTriples(const OrdinalType[], const OrdinalType[], const ScalarType[], const int, char[], const int, int &, const ::Teuchos::Comm< int > &, std::ostream *errStrm=NULL)
Pack matrix entries ("triples" (i, j, A(i,j))) into the given output buffer.