42#ifndef TPETRA_IMPORT_UTIL2_HPP
43#define TPETRA_IMPORT_UTIL2_HPP
50#include "Tpetra_ConfigDefs.hpp"
51#include "Tpetra_Import.hpp"
52#include "Tpetra_HashTable.hpp"
53#include "Tpetra_Map.hpp"
55#include "Tpetra_Distributor.hpp"
58#include "Tpetra_Vector.hpp"
59#include "Kokkos_DualView.hpp"
60#include <Teuchos_Array.hpp>
67namespace Import_Util {
71template<
typename Scalar,
typename Ordinal>
74 const Teuchos::ArrayView<Ordinal>& CRS_colind,
75 const Teuchos::ArrayView<Scalar>&CRS_vals);
77template<
typename Ordinal>
80 const Teuchos::ArrayView<Ordinal>& CRS_colind);
82template<
typename rowptr_array_type,
typename colind_array_type,
typename vals_array_type>
85 const colind_array_type& CRS_colind,
86 const vals_array_type& CRS_vals);
88template<
typename rowptr_array_type,
typename colind_array_type>
91 const colind_array_type& CRS_colind);
97template<
typename Scalar,
typename Ordinal>
100 const Teuchos::ArrayView<Ordinal>& CRS_colind,
101 const Teuchos::ArrayView<Scalar>& CRS_vals);
103template<
typename Ordinal>
106 const Teuchos::ArrayView<Ordinal>& CRS_colind);
123template <
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
126 const Teuchos::ArrayView<LocalOrdinal> &columnIndices_LID,
127 const Teuchos::ArrayView<GlobalOrdinal> &columnIndices_GID,
129 const Teuchos::ArrayView<const int> &owningPids,
130 Teuchos::Array<int> &remotePids,
149 template <
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
151 bool useReverseModeForOwnership,
152 const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& transferForMigratingData,
153 bool useReverseModeForMigration,
165namespace Import_Util {
168template<
typename PID,
typename GlobalOrdinal>
169bool sort_PID_then_GID(
const std::pair<PID,GlobalOrdinal> &a,
170 const std::pair<PID,GlobalOrdinal> &b)
173 return (a.first < b.first);
174 return (a.second < b.second);
177template<
typename PID,
178 typename GlobalOrdinal,
179 typename LocalOrdinal>
180bool sort_PID_then_pair_GID_LID(
const std::pair<PID, std::pair< GlobalOrdinal, LocalOrdinal > > &a,
181 const std::pair<PID, std::pair< GlobalOrdinal, LocalOrdinal > > &b)
184 return a.first < b.first;
186 return (a.second.first < b.second.first);
189template<
typename Scalar,
190 typename LocalOrdinal,
191 typename GlobalOrdinal,
194reverseNeighborDiscovery(
const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & SourceMatrix,
195 const typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::row_ptrs_host_view_type & rowptr,
196 const typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_inds_host_view_type & colind,
200 Teuchos::ArrayRCP<int>& type3PIDs,
201 Teuchos::ArrayRCP<LocalOrdinal>& type3LIDs,
202 Teuchos::RCP<
const Teuchos::Comm<int> >& rcomm)
204#ifdef HAVE_TPETRACORE_MPI
205 using Teuchos::TimeMonitor;
206 using ::Tpetra::Details::Behavior;
207 using Kokkos::AllowPadding;
208 using Kokkos::view_alloc;
209 using Kokkos::WithoutInitializing;
210 typedef LocalOrdinal LO;
211 typedef GlobalOrdinal GO;
212 typedef std::pair<GO,GO> pidgidpair_t;
214 const std::string prefix {
" Import_Util2::ReverseND:: "};
215 const std::string label (
"IU2::Neighbor");
218 if(MyImporter.is_null())
return;
220 std::ostringstream errstr;
222 auto const comm = MyDomainMap->getComm();
224 MPI_Comm rawComm = getRawMpiComm(*comm);
225 const int MyPID = rcomm->getRank ();
233 Distributor & Distor = MyImporter->getDistributor();
234 const size_t NumRecvs = Distor.getNumReceives();
235 const size_t NumSends = Distor.getNumSends();
236 auto RemoteLIDs = MyImporter->getRemoteLIDs();
237 auto const ProcsFrom = Distor.getProcsFrom();
238 auto const ProcsTo = Distor.getProcsTo();
240 auto LengthsFrom = Distor.getLengthsFrom();
241 auto MyColMap = SourceMatrix.getColMap();
242 const size_t numCols = MyColMap->getLocalNumElements ();
243 RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > target = MyImporter->getTargetMap();
247 Teuchos::Array<int> RemotePIDOrder(numCols,-1);
250 for (
size_t i = 0, j = 0; i < NumRecvs; ++i) {
251 for (
size_t k = 0; k < LengthsFrom[i]; ++k) {
252 const int pid = ProcsFrom[i];
254 RemotePIDOrder[RemoteLIDs[j]] = i;
266 Teuchos::Array<int> ReverseSendSizes(NumRecvs,0);
268 Teuchos::Array< Teuchos::ArrayRCP<pidgidpair_t > > RSB(NumRecvs);
271#ifdef HAVE_TPETRA_MMM_TIMINGS
272 TimeMonitor set_all(*TimeMonitor::getNewTimer(prefix + std::string(
"isMMallSetRSB")));
285 Teuchos::Array<std::set<pidgidpair_t>> pidsets(NumRecvs);
287#ifdef HAVE_TPETRA_MMM_TIMINGS
288 TimeMonitor set_insert(*TimeMonitor::getNewTimer(prefix + std::string(
"isMMallSetRSBinsert")));
290 for(
size_t i=0; i < NumExportLIDs; i++) {
291 LO lid = ExportLIDs[i];
292 GO exp_pid = ExportPIDs[i];
293 for(
auto j=rowptr[lid]; j<rowptr[lid+1]; j++){
294 int pid_order = RemotePIDOrder[colind[j]];
296 GO gid = MyColMap->getGlobalElement(colind[j]);
297 auto tpair = pidgidpair_t(exp_pid,gid);
298 pidsets[pid_order].insert(pidsets[pid_order].end(),tpair);
305#ifdef HAVE_TPETRA_MMM_TIMINGS
306 TimeMonitor set_cpy(*TimeMonitor::getNewTimer(prefix + std::string(
"isMMallSetRSBcpy")));
309 for(
auto &&ps : pidsets) {
311 RSB[jj] = Teuchos::arcp(
new pidgidpair_t[ s ],0, s ,
true);
312 std::copy(ps.begin(),ps.end(),RSB[jj]);
313 ReverseSendSizes[jj]=s;
319 Teuchos::Array<int> ReverseRecvSizes(NumSends,-1);
320 Teuchos::Array<MPI_Request> rawBreq(ProcsFrom.size()+ProcsTo.size(), MPI_REQUEST_NULL);
322 const int mpi_tag_base_ = 3;
325 for(
int i=0;i<ProcsTo.size();++i) {
326 int Rec_Tag = mpi_tag_base_ + ProcsTo[i];
327 int * thisrecv = (
int *) (&ReverseRecvSizes[i]);
328 MPI_Request rawRequest = MPI_REQUEST_NULL;
329 MPI_Irecv(
const_cast<int*
>(thisrecv),
336 rawBreq[mpireq_idx++]=rawRequest;
338 for(
int i=0;i<ProcsFrom.size();++i) {
339 int Send_Tag = mpi_tag_base_ + MyPID;
340 int * mysend = (
int *)(&ReverseSendSizes[i]);
341 MPI_Request rawRequest = MPI_REQUEST_NULL;
349 rawBreq[mpireq_idx++]=rawRequest;
351 Teuchos::Array<MPI_Status> rawBstatus(rawBreq.size());
352#ifdef HAVE_TPETRA_DEBUG
355 MPI_Waitall (rawBreq.size(), rawBreq.getRawPtr(),
356 rawBstatus.getRawPtr());
359#ifdef HAVE_TPETRA_DEBUG
361 errstr <<MyPID<<
"sE1 reverseNeighborDiscovery Mpi_Waitall error on send ";
363 std::cerr<<errstr.str()<<std::flush;
367 int totalexportpairrecsize = 0;
368 for (
size_t i = 0; i < NumSends; ++i) {
369 totalexportpairrecsize += ReverseRecvSizes[i];
370#ifdef HAVE_TPETRA_DEBUG
371 if(ReverseRecvSizes[i]<0) {
372 errstr << MyPID <<
"E4 reverseNeighborDiscovery: got < 0 for receive size "<<ReverseRecvSizes[i]<<std::endl;
377 Teuchos::ArrayRCP<pidgidpair_t >AllReverseRecv= Teuchos::arcp(
new pidgidpair_t[totalexportpairrecsize],0,totalexportpairrecsize,
true);
380 for(
int i=0;i<ProcsTo.size();++i) {
381 int recv_data_size = ReverseRecvSizes[i]*2;
382 int recvData_MPI_Tag = mpi_tag_base_*2 + ProcsTo[i];
383 MPI_Request rawRequest = MPI_REQUEST_NULL;
384 GO * rec_bptr = (GO*) (&AllReverseRecv[offset]);
385 offset+=ReverseRecvSizes[i];
388 ::Tpetra::Details::MpiTypeTraits<GO>::getType(rec_bptr[0]),
393 rawBreq[mpireq_idx++]=rawRequest;
395 for(
int ii=0;ii<ProcsFrom.size();++ii) {
396 GO * send_bptr = (GO*) (RSB[ii].getRawPtr());
397 MPI_Request rawSequest = MPI_REQUEST_NULL;
398 int send_data_size = ReverseSendSizes[ii]*2;
399 int sendData_MPI_Tag = mpi_tag_base_*2+MyPID;
402 ::Tpetra::Details::MpiTypeTraits<GO>::getType(send_bptr[0]),
408 rawBreq[mpireq_idx++]=rawSequest;
410#ifdef HAVE_TPETRA_DEBUG
413 MPI_Waitall (rawBreq.size(),
415 rawBstatus.getRawPtr());
416#ifdef HAVE_TPETRA_DEBUG
418 errstr <<MyPID<<
"E3.r reverseNeighborDiscovery Mpi_Waitall error on receive ";
420 std::cerr<<errstr.str()<<std::flush;
423 std::sort(AllReverseRecv.begin(), AllReverseRecv.end(), Tpetra::Import_Util::sort_PID_then_GID<GlobalOrdinal, GlobalOrdinal>);
425 auto newEndOfPairs = std::unique(AllReverseRecv.begin(), AllReverseRecv.end());
427 if(AllReverseRecv.begin() == newEndOfPairs)
return;
428 int ARRsize = std::distance(AllReverseRecv.begin(),newEndOfPairs);
429 auto rPIDs = Teuchos::arcp(
new int[ARRsize],0,ARRsize,
true);
430 auto rLIDs = Teuchos::arcp(
new LocalOrdinal[ARRsize],0,ARRsize,
true);
433 for(
auto itr = AllReverseRecv.begin(); itr!=newEndOfPairs; ++itr ) {
434 if((
int)(itr->first) != MyPID) {
435 rPIDs[tsize]=(int)itr->first;
436 LocalOrdinal lid = MyDomainMap->getLocalElement(itr->second);
442 type3PIDs=rPIDs.persistingView(0,tsize);
443 type3LIDs=rLIDs.persistingView(0,tsize);
446 std::cerr<<errstr.str()<<std::flush;
450 MPI_Abort (MPI_COMM_WORLD, -1);
456template<
typename Scalar,
typename Ordinal>
459 const Teuchos::ArrayView<Ordinal> &
CRS_colind,
460 const Teuchos::ArrayView<Scalar> &
CRS_vals)
472 if(start >=
nnz)
continue;
482 while (
m<
n)
m =
m*3+1;
506template<
typename Ordinal>
508sortCrsEntries (
const Teuchos::ArrayView<size_t> &
CRS_rowptr,
509 const Teuchos::ArrayView<Ordinal> &
CRS_colind)
512 Teuchos::ArrayView<Tpetra::Details::DefaultTypes::scalar_type>
CRS_vals;
518template<
class RowOffsetsType,
class ColumnIndicesType,
class ValuesType>
519class SortCrsEntries {
521 typedef typename ColumnIndicesType::non_const_value_type ordinal_type;
522 typedef typename ValuesType::non_const_value_type scalar_type;
525 SortCrsEntries (
const RowOffsetsType& ptr,
526 const ColumnIndicesType& ind,
527 const ValuesType& val) :
532 static_assert (std::is_signed<ordinal_type>::value,
"The type of each "
533 "column index -- that is, the type of each entry of ind "
534 "-- must be signed in order for this functor to work.");
537 KOKKOS_FUNCTION
void operator() (
const size_t i)
const
539 const size_t nnz = ind_.extent (0);
540 const size_t start = ptr_(i);
541 const bool permute_values_array = val_.extent(0) > 0;
544 const size_t NumEntries = ptr_(i+1) - start;
546 const ordinal_type n =
static_cast<ordinal_type
> (NumEntries);
548 while (m<n) m = m*3+1;
552 ordinal_type max = n - m;
553 for (ordinal_type j = 0; j < max; j++) {
554 for (ordinal_type k = j; k >= 0; k -= m) {
555 const size_t sk = start+k;
556 if (ind_(sk+m) >= ind_(sk)) {
559 if (permute_values_array) {
560 const scalar_type dtemp = val_(sk+m);
561 val_(sk+m) = val_(sk);
564 const ordinal_type itemp = ind_(sk+m);
565 ind_(sk+m) = ind_(sk);
575 sortCrsEntries (
const RowOffsetsType& ptr,
576 const ColumnIndicesType& ind,
577 const ValuesType& val)
584 if (ptr.extent (0) == 0) {
587 const size_t NumRows = ptr.extent (0) - 1;
589 typedef size_t index_type;
590 typedef typename ValuesType::execution_space execution_space;
593 typedef Kokkos::RangePolicy<execution_space, index_type> range_type;
594 Kokkos::parallel_for (
"sortCrsEntries", range_type (0, NumRows),
595 SortCrsEntries (ptr, ind, val));
600 ColumnIndicesType ind_;
606template<
typename rowptr_array_type,
typename colind_array_type,
typename vals_array_type>
609 const colind_array_type& CRS_colind,
610 const vals_array_type& CRS_vals)
612 Impl::SortCrsEntries<rowptr_array_type, colind_array_type,
613 vals_array_type>::sortCrsEntries (CRS_rowptr, CRS_colind, CRS_vals);
616template<
typename rowptr_array_type,
typename colind_array_type>
619 const colind_array_type& CRS_colind)
622 typedef typename colind_array_type::execution_space execution_space;
623 typedef Tpetra::Details::DefaultTypes::scalar_type scalar_type;
624 typedef typename Kokkos::View<scalar_type*, execution_space> scalar_view_type;
625 scalar_view_type CRS_vals;
627 scalar_view_type>(CRS_rowptr, CRS_colind, CRS_vals);
631template<
typename Scalar,
typename Ordinal>
634 const Teuchos::ArrayView<Ordinal> &
CRS_colind,
635 const Teuchos::ArrayView<Scalar> &
CRS_vals)
706template<
typename Ordinal>
708sortAndMergeCrsEntries (
const Teuchos::ArrayView<size_t> &
CRS_rowptr,
709 const Teuchos::ArrayView<Ordinal> &
CRS_colind)
711 Teuchos::ArrayView<Tpetra::Details::DefaultTypes::scalar_type>
CRS_vals;
716template <
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
719 const Teuchos::ArrayView<LocalOrdinal> &
colind_LID,
720 const Teuchos::ArrayView<GlobalOrdinal> &
colind_GID,
722 const Teuchos::ArrayView<const int> &
owningPIDs,
723 Teuchos::Array<int> &remotePIDs,
731 const char prefix[] =
"lowCommunicationMakeColMapAndReindex: ";
778 const LO LID = domainMap.getLocalElement (GID);
792 PID == -1, std::invalid_argument,
prefix <<
"Cannot figure out if "
809 if (domainMap.getComm ()->getSize () == 1) {
814 "process in the domain Map's communicator, which means that there are no "
815 "\"remote\" indices. Nevertheless, some column indices are not in the "
924 std::runtime_error,
prefix <<
"Local ID count test failed.");
928 const GST minus_one = Teuchos::OrdinalTraits<GST>::invalid ();
930 domainMap.getComm ()));
967template <
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
991#ifdef HAVE_TPETRA_DEBUG
1002 Teuchos::ArrayRCP<int>
pids =
temp.getDataNonConst();
1007 else {
TEUCHOS_TEST_FOR_EXCEPTION(1,std::runtime_error,
"Tpetra::Import_Util::getTwoTransferOwnershipVector owner must be 1-to-1");}
Declaration of the Tpetra::CrsMatrix class.
Add specializations of Teuchos::Details::MpiTypeTraits for Kokkos::complex<float> and Kokkos::complex...
Declaration and definition of Tpetra::Details::reallocDualViewIfNeeded, an implementation detail of T...
void lowCommunicationMakeColMapAndReindex(const Teuchos::ArrayView< const size_t > &rowPointers, const Teuchos::ArrayView< LocalOrdinal > &columnIndices_LID, const Teuchos::ArrayView< GlobalOrdinal > &columnIndices_GID, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &domainMap, const Teuchos::ArrayView< const int > &owningPids, Teuchos::Array< int > &remotePids, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &colMap)
lowCommunicationMakeColMapAndReindex
void getTwoTransferOwnershipVector(const ::Tpetra::Details::Transfer< LocalOrdinal, GlobalOrdinal, Node > &transferThatDefinesOwnership, bool useReverseModeForOwnership, const ::Tpetra::Details::Transfer< LocalOrdinal, GlobalOrdinal, Node > &transferForMigratingData, bool useReverseModeForMigration, Tpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node > &owningPIDs)
Generates an list of owning PIDs based on two transfer (aka import/export objects) Let: OwningMap = u...
void sortAndMergeCrsEntries(const Teuchos::ArrayView< size_t > &CRS_rowptr, const Teuchos::ArrayView< Ordinal > &CRS_colind, const Teuchos::ArrayView< Scalar > &CRS_vals)
Sort and merge the entries of the (raw CSR) matrix by column index within each row.
void sortCrsEntries(const Teuchos::ArrayView< size_t > &CRS_rowptr, const Teuchos::ArrayView< Ordinal > &CRS_colind, const Teuchos::ArrayView< Scalar > &CRS_vals)
Sort the entries of the (raw CSR) matrix by column index within each row.
void getPids(const Tpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer, Teuchos::Array< int > &pids, bool use_minus_one_for_local)
Like getPidGidPairs, but just gets the PIDs, ordered by the column Map.
Stand-alone utility functions and macros.
Struct that holds views of the contents of a CrsMatrix.
Common base class of Import and Export.
Teuchos::ArrayView< const LO > getExportLIDs() const
List of entries in the source Map that will be sent to other processes.
size_t getNumExportIDs() const
Number of entries that must be sent by the calling process to other processes.
Teuchos::ArrayView< const int > getExportPIDs() const
List of processes to which entries will be sent.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
void sort3(const IT1 &first1, const IT1 &last1, const IT2 &first2, const IT3 &first3)
Sort the first array, and apply the same permutation to the second and third arrays.
void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2)
Sort the first array, and apply the resulting permutation to the second array.
@ REPLACE
Replace existing values with new values.