1 #ifndef TPETRA_DETAILS_MAKECOLMAP_DEF_HPP
2 #define TPETRA_DETAILS_MAKECOLMAP_DEF_HPP
24 #include "Tpetra_RowGraph.hpp"
26 #include "Teuchos_Array.hpp"
27 #include "Kokkos_Bitset.hpp"
34 template <
class LO,
class GO,
class NT>
36 Teuchos::Array<int>& remotePIDs,
38 size_t numLocalColGIDs,
39 size_t numRemoteColGIDs,
40 std::set<GO>& RemoteGIDSet,
41 std::vector<GO>& RemoteGIDUnorderedVector,
42 std::vector<bool>& GIDisLocal,
43 const bool sortEachProcsGids,
44 std::ostream* errStrm) {
47 using Teuchos::ArrayView;
50 const char prefix[] =
"Tpetra::Details::makeColMapImpl: ";
51 typedef ::Tpetra::Map<LO, GO, NT> map_type;
81 if (domMap->getComm()->getSize() == 1) {
82 if (numRemoteColGIDs != 0) {
84 if (errStrm != NULL) {
85 *errStrm << prefix <<
"The domain Map only has one process, but "
86 << numRemoteColGIDs <<
" column "
87 << (numRemoteColGIDs != 1 ?
"indices are" :
"index is")
88 <<
" not in the domain Map. Either these indices are "
89 "invalid or the domain Map is invalid. Remember that nonsquare "
90 "matrices, or matrices where the row and range Maps differ, "
91 "require calling the version of fillComplete that takes the "
92 "domain and range Maps as input."
96 if (numLocalColGIDs == domMap->getLocalNumElements()) {
106 Array<GO> myColumns(numLocalColGIDs + numRemoteColGIDs);
108 ArrayView<GO> LocalColGIDs = myColumns(0, numLocalColGIDs);
109 ArrayView<GO> remoteColGIDs = myColumns(numLocalColGIDs, numRemoteColGIDs);
112 if (sortEachProcsGids) {
114 std::copy(RemoteGIDSet.begin(), RemoteGIDSet.end(),
115 remoteColGIDs.begin());
118 std::copy(RemoteGIDUnorderedVector.begin(),
119 RemoteGIDUnorderedVector.end(), remoteColGIDs.begin());
125 if (static_cast<size_t>(remotePIDs.size()) != numRemoteColGIDs) {
126 remotePIDs.resize(numRemoteColGIDs);
131 domMap->getRemoteIndexList(remoteColGIDs, remotePIDs());
141 if (errStrm != NULL) {
142 *errStrm << prefix <<
"Some column indices are not in the domain Map."
143 "Either these column indices are invalid or the domain Map is "
144 "invalid. Likely cause: For a nonsquare matrix, you must give the "
145 "domain and range Maps as input to fillComplete."
165 sort2(remotePIDs.begin(), remotePIDs.end(), remoteColGIDs.begin(),
true);
177 const size_t numDomainElts = domMap->getLocalNumElements();
178 if (numLocalColGIDs == numDomainElts) {
182 if (domMap->isContiguous()) {
187 GO curColMapGid = domMap->getMinGlobalIndex();
188 for (
size_t k = 0; k < numLocalColGIDs; ++k, ++curColMapGid) {
189 LocalColGIDs[k] = curColMapGid;
192 ArrayView<const GO> domainElts = domMap->getLocalElementList();
193 std::copy(domainElts.begin(), domainElts.end(), LocalColGIDs.begin());
198 size_t numLocalCount = 0;
199 if (domMap->isContiguous()) {
204 GO curColMapGid = domMap->getMinGlobalIndex();
205 for (
size_t i = 0; i < numDomainElts; ++i, ++curColMapGid) {
207 LocalColGIDs[numLocalCount++] = curColMapGid;
211 ArrayView<const GO> domainElts = domMap->getLocalElementList();
212 for (
size_t i = 0; i < numDomainElts; ++i) {
214 LocalColGIDs[numLocalCount++] = domainElts[i];
218 if (numLocalCount != numLocalColGIDs) {
219 if (errStrm != NULL) {
220 *errStrm << prefix <<
"numLocalCount = " << numLocalCount
221 <<
" != numLocalColGIDs = " << numLocalColGIDs
222 <<
". This should never happen. "
223 "Please report this bug to the Tpetra developers."
280 Tpetra::Details::OrdinalTraits<global_size_t>::invalid();
285 const GO indexBase = domMap->getIndexBase();
286 colMap = rcp(
new map_type(INV, myColumns, indexBase, domMap->getComm()));
290 template <
class LO,
class GO,
class NT>
292 Teuchos::Array<int>& remotePIDs,
295 const bool sortEachProcsGids,
296 std::ostream* errStrm) {
298 using Teuchos::Array;
299 using Teuchos::ArrayView;
301 typedef ::Tpetra::Map<LO, GO, NT> map_type;
302 const char prefix[] =
"Tpetra::Details::makeColMap: ";
312 if (domMap.is_null() || domMap->getComm().is_null()) {
313 colMap = Teuchos::null;
322 if (colMap.is_null() || !graph.
hasColMap()) {
324 if (errStrm != NULL) {
325 *errStrm << prefix <<
"The graph is locally indexed on the calling "
326 "process, but has no column Map (either getColMap() returns null, "
327 "or hasColMap() returns false)."
339 if (colMap->isContiguous()) {
341 const LO numCurGids =
static_cast<LO
>(colMap->getLocalNumElements());
342 myColumns.resize(numCurGids);
343 const GO myFirstGblInd = colMap->getMinGlobalIndex();
344 for (LO k = 0; k < numCurGids; ++k) {
345 myColumns[k] = myFirstGblInd +
static_cast<GO
>(k);
348 ArrayView<const GO> curGids = graph.
getColMap()->getLocalElementList();
350 const LO numCurGids =
static_cast<LO
>(curGids.size());
351 myColumns.resize(numCurGids);
352 for (LO k = 0; k < numCurGids; ++k) {
353 myColumns[k] = curGids[k];
381 const LO LINV = Tpetra::Details::OrdinalTraits<LO>::invalid();
382 size_t numLocalColGIDs = 0;
383 size_t numRemoteColGIDs = 0;
387 std::vector<bool> GIDisLocal(domMap->getLocalNumElements(),
false);
388 std::set<GO> RemoteGIDSet;
391 std::vector<GO> RemoteGIDUnorderedVector;
396 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
398 typename RowGraph<LO, GO, NT>::global_inds_host_view_type rowGids;
401 const LO numEnt =
static_cast<LO
>(rowGids.size());
403 for (LO k = 0; k < numEnt; ++k) {
404 const GO gid = rowGids[k];
405 const LO lid = domMap->getLocalElement(gid);
407 const bool alreadyFound = GIDisLocal[lid];
409 GIDisLocal[lid] =
true;
413 const bool notAlreadyFound = RemoteGIDSet.insert(gid).second;
414 if (notAlreadyFound) {
415 if (!sortEachProcsGids) {
422 RemoteGIDUnorderedVector.push_back(gid);
432 return makeColMapImpl<LO, GO, NT>(
435 numLocalColGIDs, numRemoteColGIDs,
436 RemoteGIDSet, RemoteGIDUnorderedVector, GIDisLocal,
437 sortEachProcsGids, errStrm);
449 Tpetra::Details::OrdinalTraits<global_size_t>::invalid();
454 const GO indexBase = domMap->getIndexBase();
455 colMap = rcp(
new map_type(INV, myColumns, indexBase, domMap->getComm()));
459 template <
typename GOView,
typename bitset_t>
460 struct GatherPresentEntries {
461 using GO =
typename GOView::non_const_value_type;
463 GatherPresentEntries(GO minGID_,
const GOView& gids_,
const bitset_t& present_)
466 , present(present_) {}
468 KOKKOS_INLINE_FUNCTION
void operator()(
const GO i)
const {
469 present.set(gids(i) - minGID);
477 template <
typename LO,
typename GO,
typename device_t,
typename LocalMapType,
typename const_bitset_t,
bool doingRemotes>
479 using mem_space =
typename device_t::memory_space;
480 using GOView = Kokkos::View<GO*, mem_space>;
481 using SingleView = Kokkos::View<GO, mem_space>;
483 ListGIDs(GO minGID_, GOView& gidList_, SingleView& numElems_, const_bitset_t& present_,
const LocalMapType& localDomainMap_)
486 , numElems(numElems_)
488 , localDomainMap(localDomainMap_) {}
490 KOKKOS_INLINE_FUNCTION
void operator()(
const GO i, GO& lcount,
const bool finalPass)
const {
491 bool isRemote = localDomainMap.getLocalElement(i + minGID) == ::Tpetra::Details::OrdinalTraits<LO>::invalid();
492 if (present.test(i) && doingRemotes == isRemote) {
495 gidList(lcount) = minGID + i;
499 if ((i == static_cast<GO>(present.size() - 1)) && finalPass) {
508 const_bitset_t present;
509 const LocalMapType localDomainMap;
512 template <
typename GO,
typename mem_space>
513 struct MinMaxReduceFunctor {
514 using MinMaxValue =
typename Kokkos::MinMax<GO>::value_type;
515 using GOView = Kokkos::View<GO*, mem_space>;
517 MinMaxReduceFunctor(
const GOView& gids_)
520 KOKKOS_INLINE_FUNCTION
void operator()(
const GO i, MinMaxValue& lminmax)
const {
522 if (gid < lminmax.min_val)
523 lminmax.min_val = gid;
524 if (gid > lminmax.max_val)
525 lminmax.max_val = gid;
531 template <
class LO,
class GO,
class NT>
534 Kokkos::View<GO*, typename NT::memory_space> gids,
535 std::ostream* errStrm) {
536 using Kokkos::RangePolicy;
537 using Teuchos::Array;
539 using device_t =
typename NT::device_type;
540 using exec_space =
typename device_t::execution_space;
541 using memory_space =
typename device_t::memory_space;
547 using bitset_t = Kokkos::Bitset<typename exec_space::memory_space>;
548 using const_bitset_t = Kokkos::ConstBitset<typename exec_space::memory_space>;
549 using GOView = Kokkos::View<GO*, memory_space>;
550 using SingleView = Kokkos::View<GO, memory_space>;
552 using LocalMap =
typename map_type::local_map_type;
553 GO nentries = gids.extent(0);
554 GO minGID = Teuchos::OrdinalTraits<GO>::max();
556 using MinMaxValue =
typename Kokkos::MinMax<GO>::value_type;
557 MinMaxValue minMaxGID = {minGID, maxGID};
558 Kokkos::parallel_reduce(RangePolicy<exec_space>(0, nentries),
559 MinMaxReduceFunctor<GO, memory_space>(gids),
560 Kokkos::MinMax<GO>(minMaxGID));
561 minGID = minMaxGID.min_val;
562 maxGID = minMaxGID.max_val;
565 bitset_t presentGIDs(maxGID - minGID + 1);
566 Kokkos::parallel_for(RangePolicy<exec_space>(0, nentries), GatherPresentEntries<GOView, bitset_t>(minGID, gids, presentGIDs));
567 const_bitset_t constPresentGIDs(presentGIDs);
569 SingleView numLocals(
"Num local GIDs");
570 SingleView numRemotes(
"Num remote GIDs");
571 GOView localGIDView(Kokkos::ViewAllocateWithoutInitializing(
"Local GIDs"), constPresentGIDs.count());
572 GOView remoteGIDView(Kokkos::ViewAllocateWithoutInitializing(
"Remote GIDs"), constPresentGIDs.count());
573 LocalMap localDomMap = domMap->getLocalMap();
575 Kokkos::parallel_scan(RangePolicy<exec_space>(0, constPresentGIDs.size()),
576 ListGIDs<LO, GO, device_t, LocalMap, const_bitset_t, false>(minGID, localGIDView, numLocals, constPresentGIDs, localDomMap));
578 Kokkos::parallel_scan(RangePolicy<exec_space>(0, constPresentGIDs.size()),
579 ListGIDs<LO, GO, device_t, LocalMap, const_bitset_t, true>(minGID, remoteGIDView, numRemotes, constPresentGIDs, localDomMap));
581 GO numLocalColGIDs = 0;
582 GO numRemoteColGIDs = 0;
588 auto localsHost = Kokkos::create_mirror_view(localGIDView);
589 auto remotesHost = Kokkos::create_mirror_view(remoteGIDView);
595 Kokkos::fence(
"Tpetra::makeColMap");
597 std::set<GO> RemoteGIDSet;
598 std::vector<GO> RemoteGIDUnorderedVector;
599 std::vector<bool> GIDisLocal(domMap->getLocalNumElements(),
false);
600 for (GO i = 0; i < numLocalColGIDs; i++) {
601 GO gid = localsHost(i);
604 GIDisLocal[domMap->getLocalElement(gid)] =
true;
606 RemoteGIDUnorderedVector.reserve(numRemoteColGIDs);
607 for (GO i = 0; i < numRemoteColGIDs; i++) {
608 GO gid = remotesHost(i);
609 RemoteGIDSet.insert(gid);
610 RemoteGIDUnorderedVector.push_back(gid);
613 Array<int> remotePIDs;
615 return makeColMapImpl<LO, GO, NT>(
619 static_cast<size_t>(numLocalColGIDs),
620 static_cast<size_t>(numRemoteColGIDs),
622 RemoteGIDUnorderedVector,
636 #define TPETRA_DETAILS_MAKECOLMAP_INSTANT(LO, GO, NT) \
637 namespace Details { \
639 makeColMap(Teuchos::RCP<const Tpetra::Map<LO, GO, NT>>&, \
640 Teuchos::Array<int>&, \
641 const Teuchos::RCP<const Tpetra::Map<LO, GO, NT>>&, \
642 const RowGraph<LO, GO, NT>&, \
646 makeColMap(Teuchos::RCP<const Tpetra::Map<LO, GO, NT>>&, \
647 const Teuchos::RCP<const Tpetra::Map<LO, GO, NT>>&, \
648 Kokkos::View<GO*, typename NT::memory_space>, \
652 #endif // TPETRA_DETAILS_MAKECOLMAP_DEF_HPP
virtual bool hasColMap() const =0
Whether the graph has a well-defined column Map.
An abstract interface for graphs accessed by rows.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
The Map that describes this graph's distribution of rows over processes.
LookupStatus
Return status of Map remote index lookup (getRemoteIndexList()).
size_t getLocalNumElements() const
The number of elements belonging to the calling process.
int makeColMap(Teuchos::RCP< const Tpetra::Map< LO, GO, NT >> &colMap, Teuchos::Array< int > &remotePIDs, const Teuchos::RCP< const Tpetra::Map< LO, GO, NT >> &domMap, const RowGraph< LO, GO, NT > &graph, const bool sortEachProcsGids=true, std::ostream *errStrm=NULL)
Make the graph's column Map.
"Local" part of Map suitable for Kokkos kernels.
size_t global_size_t
Global size_t object.
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.
global_ordinal_type getGlobalElement(local_ordinal_type localIndex) const
The global index corresponding to the given local index.
void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2, const bool stableSort=false)
Sort the first array, and apply the resulting permutation to the second array.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getColMap() const =0
The Map that describes this graph's distribution of columns over processes.
virtual bool isGloballyIndexed() const =0
If graph indices are in the global range, this function returns true. Otherwise, this function return...
virtual void getGlobalRowView(const GlobalOrdinal gblRow, global_inds_host_view_type &gblColInds) const =0
Get a const, non-persisting view of the given global row's global column indices, as a Teuchos::Array...
A parallel distribution of indices over processes.
Stand-alone utility functions and macros.
virtual bool isLocallyIndexed() const =0
If graph indices are in the local range, this function returns true. Otherwise, this function returns...