Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Tpetra_Import_def.hpp
1 // @HEADER
2 // *****************************************************************************
3 // Tpetra: Templated Linear Algebra Services Package
4 //
5 // Copyright 2008 NTESS and the Tpetra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef TPETRA_IMPORT_DEF_HPP
11 #define TPETRA_IMPORT_DEF_HPP
12 
13 #include "Tpetra_Distributor.hpp"
14 #include "Tpetra_Map.hpp"
15 #include "Tpetra_ImportExportData.hpp"
16 #include "Tpetra_Util.hpp"
17 #include "Tpetra_Import_Util.hpp"
18 #include "Tpetra_Export.hpp"
20 #include "Tpetra_Details_DualViewUtil.hpp"
23 #include "Teuchos_as.hpp"
24 #ifdef HAVE_TPETRA_MMM_TIMINGS
25 #include "Teuchos_TimeMonitor.hpp"
26 #endif
27 #include <array>
28 #include <memory>
29 
30 namespace Teuchos {
31 template <class T>
32 std::string toString(const std::vector<T>& x) {
33  std::ostringstream os;
34  os << "[";
35  const std::size_t N = x.size();
36  for (std::size_t k = 0; k < N; ++k) {
37  os << x[k];
38  if (k + std::size_t(1) < N) {
39  os << ",";
40  }
41  }
42  os << "]";
43  return os.str();
44 }
45 
46 template <class ElementType, class DeviceType>
47 std::string toString(const Kokkos::View<const ElementType*, DeviceType>& x) {
48  std::ostringstream os;
49  os << "[";
50  const std::size_t N = std::size_t(x.extent(0));
51  for (std::size_t k = 0; k < N; ++k) {
52  os << x[k];
53  if (k + std::size_t(1) < N) {
54  os << ",";
55  }
56  }
57  os << "]";
58  return os.str();
59 }
60 } // namespace Teuchos
61 
62 namespace Tpetra {
63 
64 // head: init(source, target, true, remotePIDs, Teuchos::null);
65 
66 template <class LocalOrdinal, class GlobalOrdinal, class Node>
67 void Import<LocalOrdinal, GlobalOrdinal, Node>::
68  init(const Teuchos::RCP<const map_type>& source,
69  const Teuchos::RCP<const map_type>& /* target */,
70  bool useRemotePIDs,
71  Teuchos::Array<int>& remotePIDs,
72  const Teuchos::RCP<Teuchos::ParameterList>& plist) {
73  using std::endl;
74  using Teuchos::Array;
75  using Teuchos::null;
76  using Teuchos::Ptr;
77  using Teuchos::rcp;
78  using ::Tpetra::Details::ProfilingRegion;
79  ProfilingRegion regionImportInit("Tpetra::Import::init");
80 
81  std::unique_ptr<std::string> verbPrefix;
82  if (this->verbose()) {
83  std::ostringstream os;
84  const int myRank = source->getComm()->getRank();
85  os << "Proc " << myRank << ": Tpetra::Import::init: ";
86  verbPrefix = std::unique_ptr<std::string>(new std::string(os.str()));
87  os << endl;
88  this->verboseOutputStream() << os.str();
89  }
90 
91  Array<GlobalOrdinal> remoteGIDs;
92 
93 #ifdef HAVE_TPETRA_MMM_TIMINGS
94  using Teuchos::TimeMonitor;
95  std::string label;
96  if (!plist.is_null())
97  label = plist->get("Timer Label", label);
98  std::string prefix = std::string("Tpetra ") + label + std::string(":iport_ctor:preIData: ");
99 #else
100  (void)plist;
101 #endif
102  {
103 #ifdef HAVE_TPETRA_MMM_TIMINGS
104  auto MM(*TimeMonitor::getNewTimer(prefix));
105 #endif
106  if (this->verbose()) {
107  std::ostringstream os;
108  os << *verbPrefix << "Call setupSamePermuteRemote" << endl;
109  this->verboseOutputStream() << os.str();
110  }
111  setupSamePermuteRemote(remoteGIDs);
112  }
113  {
114 #ifdef HAVE_TPETRA_MMM_TIMINGS
115  prefix = std::string("Tpetra ") + label + std::string(":iport_ctor:preSetupExport: ");
116  auto MM2(*TimeMonitor::getNewTimer(prefix));
117 #endif
118  if (source->isDistributed()) {
119  if (this->verbose()) {
120  std::ostringstream os;
121  os << *verbPrefix << "Call setupExport" << endl;
122  this->verboseOutputStream() << os.str();
123  }
124  setupExport(remoteGIDs, useRemotePIDs, remotePIDs);
125  } else if (this->verbose()) {
126  std::ostringstream os;
127  os << *verbPrefix << "Source Map not distributed; skip setupExport"
128  << endl;
129  this->verboseOutputStream() << os.str();
130  }
131  }
132 
133  TEUCHOS_ASSERT(!this->TransferData_->permuteFromLIDs_.need_sync_device());
134  TEUCHOS_ASSERT(!this->TransferData_->permuteFromLIDs_.need_sync_host());
135  TEUCHOS_ASSERT(!this->TransferData_->permuteToLIDs_.need_sync_device());
136  TEUCHOS_ASSERT(!this->TransferData_->permuteToLIDs_.need_sync_host());
137  TEUCHOS_ASSERT(!this->TransferData_->remoteLIDs_.need_sync_device());
138  TEUCHOS_ASSERT(!this->TransferData_->remoteLIDs_.need_sync_host());
139  TEUCHOS_ASSERT(!this->TransferData_->exportLIDs_.need_sync_device());
140  TEUCHOS_ASSERT(!this->TransferData_->exportLIDs_.need_sync_host());
141 
142  this->detectRemoteExportLIDsContiguous();
143 
144  if (this->verbose()) {
145  std::ostringstream os;
146  os << *verbPrefix << "Done!" << endl;
147  this->verboseOutputStream() << os.str();
148  }
149 }
150 
151 template <class LocalOrdinal, class GlobalOrdinal, class Node>
153  Import(const Teuchos::RCP<const map_type>& source,
154  const Teuchos::RCP<const map_type>& target)
155  : base_type(source, target, Teuchos::null, Teuchos::null, "Import") {
156  Teuchos::Array<int> dummy;
157 #ifdef HAVE_TPETRA_MMM_TIMINGS
158  Teuchos::RCP<Teuchos::ParameterList> mypars = rcp(new Teuchos::ParameterList);
159  mypars->set("Timer Label", "Naive_tAFC");
160  init(source, target, false, dummy, mypars);
161 #else
162  init(source, target, false, dummy, Teuchos::null);
163 #endif
164 }
165 
166 template <class LocalOrdinal, class GlobalOrdinal, class Node>
168  Import(const Teuchos::RCP<const map_type>& source,
169  const Teuchos::RCP<const map_type>& target,
170  const Teuchos::RCP<Teuchos::FancyOStream>& out)
171  : base_type(source, target, out, Teuchos::null, "Import") {
172  Teuchos::Array<int> dummy;
173  init(source, target, false, dummy, Teuchos::null);
174 }
175 
176 template <class LocalOrdinal, class GlobalOrdinal, class Node>
178  Import(const Teuchos::RCP<const map_type>& source,
179  const Teuchos::RCP<const map_type>& target,
180  const Teuchos::RCP<Teuchos::ParameterList>& plist)
181  : base_type(source, target, Teuchos::null, plist, "Import") {
182  Teuchos::Array<int> dummy;
183  init(source, target, false, dummy, plist);
184 }
185 
186 template <class LocalOrdinal, class GlobalOrdinal, class Node>
188  Import(const Teuchos::RCP<const map_type>& source,
189  const Teuchos::RCP<const map_type>& target,
190  const Teuchos::RCP<Teuchos::FancyOStream>& out,
191  const Teuchos::RCP<Teuchos::ParameterList>& plist)
192  : base_type(source, target, out, plist, "Import") {
193  Teuchos::Array<int> dummy;
194  init(source, target, false, dummy, plist);
195 }
196 
197 template <class LocalOrdinal, class GlobalOrdinal, class Node>
199  Import(const Teuchos::RCP<const map_type>& source,
200  const Teuchos::RCP<const map_type>& target,
201  Teuchos::Array<int>& remotePIDs,
202  const Teuchos::RCP<Teuchos::ParameterList>& plist)
203  : base_type(source, target, Teuchos::null, plist, "Import") {
204  init(source, target, true, remotePIDs, plist);
205 }
206 
207 template <class LocalOrdinal, class GlobalOrdinal, class Node>
210  : base_type(rhs) {}
211 
212 template <class LocalOrdinal, class GlobalOrdinal, class Node>
215  : base_type(exporter, typename base_type::reverse_tag()) {}
216 
217 // cblcbl
218 // This is the "createExpert" version of the constructor to be used with pid/gid pairs obtained from
219 // reverse communication
220 template <class LocalOrdinal, class GlobalOrdinal, class Node>
222  Import(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node>>& source,
223  const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node>>& target,
224  const Teuchos::ArrayView<int>& userRemotePIDs,
225  const Teuchos::ArrayView<const LocalOrdinal>& userExportLIDs,
226  const Teuchos::ArrayView<const int>& userExportPIDs,
227  const Teuchos::RCP<Teuchos::ParameterList>& plist,
228  const Teuchos::RCP<Teuchos::FancyOStream>& out)
229  : base_type(source, target, out, plist, "Import") {
230  using std::endl;
231  using Teuchos::arcp;
232  using Teuchos::Array;
233  using Teuchos::ArrayRCP;
234  using Teuchos::ArrayView;
235  using Teuchos::as;
236  using Teuchos::null;
237  using Teuchos::rcp;
238  using ::Tpetra::Details::makeDualViewFromArrayView;
239  using LO = LocalOrdinal;
240  using GO = GlobalOrdinal;
241  using size_type = Teuchos::Array<int>::size_type;
242 
243  std::unique_ptr<std::string> prefix;
244  if (this->verbose()) {
245  auto comm = source.is_null() ? Teuchos::null : source->getComm();
246  const int myRank = comm.is_null() ? -1 : comm->getRank();
247  std::ostringstream os;
248  os << "Proc " << myRank << ": Tpetra::Import createExpert ctor: ";
249  prefix = std::unique_ptr<std::string>(new std::string(os.str()));
250  os << "Start" << endl;
251  this->verboseOutputStream() << os.str();
252  }
253 
254  ArrayView<const GO> sourceGIDs = source->getLocalElementList();
255  ArrayView<const GO> targetGIDs = target->getLocalElementList();
256 
257  Array<GO> tRemoteGIDs;
258  if (this->verbose()) {
259  std::ostringstream os;
260  os << *prefix << "Call setupSamePermuteRemote" << endl;
261  this->verboseOutputStream() << os.str();
262  }
263  setupSamePermuteRemote(tRemoteGIDs);
264 
265  if (this->verbose()) {
266  std::ostringstream os;
267  os << *prefix << "Sort & filter IDs" << endl;
268  this->verboseOutputStream() << os.str();
269  }
270 
271  auto tRemoteLIDs = this->TransferData_->remoteLIDs_.view_host();
272  this->TransferData_->remoteLIDs_.modify_host();
273  Teuchos::Array<int> tRemotePIDs(userRemotePIDs);
274 
275  if (this->verbose() && this->getNumRemoteIDs() > 0 && !source->isDistributed()) {
276  std::ostringstream os;
277  os << *prefix << "Target Map has remote LIDs but source Map is not "
278  "distributed. Importing to a submap of the target Map."
279  << endl;
280  this->verboseOutputStream() << os.str();
281  }
282  // FIXME (mfh 03 Feb 2019) I don't see this as "abuse"; it's
283  // perfectly valid Petra Object Model behavior.
284  TPETRA_ABUSE_WARNING(getNumRemoteIDs() > 0 && !source->isDistributed(),
285  std::runtime_error,
286  "::constructExpert: Target Map has remote LIDs but source Map "
287  "is not distributed. Importing to a submap of the target Map.");
288  TEUCHOS_TEST_FOR_EXCEPTION(tRemotePIDs.size() != tRemoteGIDs.size() ||
289  size_t(tRemoteGIDs.size()) != size_t(tRemoteLIDs.extent(0)),
290  std::runtime_error,
291  "Import::Import createExpert version: "
292  "Size mismatch on userRemotePIDs, remoteGIDs, and remoteLIDs "
293  "Array's to sort3.");
294 
295  sort3(tRemotePIDs.begin(),
296  tRemotePIDs.end(),
297  tRemoteGIDs.begin(),
298  tRemoteLIDs.data());
299 
300  // Get rid of IDs that don't exist in SourceMap
301  size_type cnt = 0;
302  size_type indexIntoRemotePIDs = tRemotePIDs.size();
303  for (size_type i = 0; i < indexIntoRemotePIDs; ++i) {
304  if (tRemotePIDs[i] == -1) {
305  ++cnt;
306  }
307  }
308 
309  if (cnt == 0) { // done modifying remoteLIDs_
310  this->TransferData_->remoteLIDs_.sync_device();
311  } else {
312  if (indexIntoRemotePIDs - cnt > 0) {
313  Array<GO> newRemoteGIDs(indexIntoRemotePIDs - cnt);
314  Array<LO> newRemoteLIDs(indexIntoRemotePIDs - cnt);
315  Array<int> newRemotePIDs(indexIntoRemotePIDs - cnt);
316  cnt = 0;
317 
318  for (size_type j = 0; j < indexIntoRemotePIDs; ++j)
319  if (tRemotePIDs[j] != -1) {
320  newRemoteGIDs[cnt] = tRemoteGIDs[j];
321  newRemotePIDs[cnt] = tRemotePIDs[j];
322  newRemoteLIDs[cnt] = target->getLocalElement(tRemoteGIDs[j]);
323  ++cnt;
324  }
325  indexIntoRemotePIDs = cnt;
326  tRemoteGIDs = newRemoteGIDs;
327  tRemotePIDs = newRemotePIDs;
328  makeDualViewFromArrayView(this->TransferData_->remoteLIDs_,
329  newRemoteLIDs().getConst(),
330  "remoteLIDs");
331  } else { // valid RemoteIDs empty
332  indexIntoRemotePIDs = 0;
333  tRemoteGIDs.clear();
334  tRemotePIDs.clear();
335  this->TransferData_->remoteLIDs_ = decltype(this->TransferData_->remoteLIDs_)();
336  }
337  }
338 
339  this->TransferData_->exportPIDs_ = Teuchos::Array<int>(userExportPIDs);
340  makeDualViewFromArrayView(this->TransferData_->exportLIDs_,
341  userExportLIDs, "exportLIDs");
342 
343  bool locallyComplete = true;
344  for (size_type i = 0; i < userExportPIDs.size() && locallyComplete; ++i) {
345  if (userExportPIDs[i] == -1) {
346  locallyComplete = false;
347  }
348  }
349  this->TransferData_->isLocallyComplete_ = locallyComplete;
350 
351  if (this->verbose()) {
352  std::ostringstream os;
353  os << *prefix << "locallyComplete: "
354  << (locallyComplete ? "true" : "false")
355  << "; call createFromSendsAndRecvs" << endl;
356  this->verboseOutputStream() << os.str();
357  }
358  {
359 #ifdef HAVE_TPETRA_MMM_TIMINGS
360  std::string mmm_prefix =
361  std::string("Tpetra ") + std::string(":iport_ctor:cFSAR ");
362 
363  auto MM3(*Teuchos::TimeMonitor::getNewTimer(mmm_prefix));
364 #endif
365  Distributor& distributor = this->TransferData_->distributor_;
366  distributor.createFromSendsAndRecvs(this->TransferData_->exportPIDs_, tRemotePIDs);
367  }
368 
369  this->detectRemoteExportLIDsContiguous();
370 
371  TEUCHOS_ASSERT(!this->TransferData_->permuteFromLIDs_.need_sync_device());
372  TEUCHOS_ASSERT(!this->TransferData_->permuteFromLIDs_.need_sync_host());
373  TEUCHOS_ASSERT(!this->TransferData_->permuteToLIDs_.need_sync_device());
374  TEUCHOS_ASSERT(!this->TransferData_->permuteToLIDs_.need_sync_host());
375  TEUCHOS_ASSERT(!this->TransferData_->remoteLIDs_.need_sync_device());
376  TEUCHOS_ASSERT(!this->TransferData_->remoteLIDs_.need_sync_host());
377  TEUCHOS_ASSERT(!this->TransferData_->exportLIDs_.need_sync_device());
378  TEUCHOS_ASSERT(!this->TransferData_->exportLIDs_.need_sync_host());
379 }
380 
381 template <class LocalOrdinal, class GlobalOrdinal, class Node>
383  Import(const Teuchos::RCP<const map_type>& source,
384  const Teuchos::RCP<const map_type>& target,
385  const size_t numSameIDs,
386  Teuchos::Array<LocalOrdinal>& permuteToLIDs,
387  Teuchos::Array<LocalOrdinal>& permuteFromLIDs,
388  Teuchos::Array<LocalOrdinal>& remoteLIDs,
389  Teuchos::Array<LocalOrdinal>& exportLIDs,
390  Teuchos::Array<int>& exportPIDs,
391  Distributor& distributor,
392  const Teuchos::RCP<Teuchos::FancyOStream>& out,
393  const Teuchos::RCP<Teuchos::ParameterList>& plist)
394  : base_type(source, target, out, plist, "Import") {
395  using std::endl;
396  using ::Tpetra::Details::makeDualViewFromArrayView;
397 
398  std::unique_ptr<std::string> prefix;
399  if (this->verbose()) {
400  auto comm = source.is_null() ? Teuchos::null : source->getComm();
401  const int myRank = comm.is_null() ? -1 : comm->getRank();
402  std::ostringstream os;
403  os << "Proc " << myRank << ": Tpetra::Import export ctor: ";
404  prefix = std::unique_ptr<std::string>(new std::string(os.str()));
405  os << "Start" << endl;
406  this->verboseOutputStream() << os.str();
407  }
408 
409  bool locallyComplete = true;
410  for (Teuchos::Array<int>::size_type i = 0; i < exportPIDs.size(); ++i) {
411  if (exportPIDs[i] == -1) {
412  locallyComplete = false;
413  }
414  }
415  if (this->verbose()) {
416  std::ostringstream os;
417  os << *prefix << "numSameIDs: " << numSameIDs << ", locallyComplete: "
418  << (locallyComplete ? "true" : "false") << endl;
419  this->verboseOutputStream() << os.str();
420  }
421 
422  this->TransferData_->isLocallyComplete_ = locallyComplete;
423  this->TransferData_->numSameIDs_ = numSameIDs;
424 
425  makeDualViewFromArrayView(this->TransferData_->permuteToLIDs_,
426  permuteToLIDs().getConst(),
427  "permuteToLIDs");
428  TEUCHOS_ASSERT(size_t(this->TransferData_->permuteToLIDs_.extent(0)) ==
429  size_t(permuteToLIDs.size()));
430  makeDualViewFromArrayView(this->TransferData_->permuteFromLIDs_,
431  permuteFromLIDs().getConst(),
432  "permuteFromLIDs");
433  TEUCHOS_ASSERT(size_t(this->TransferData_->permuteFromLIDs_.extent(0)) ==
434  size_t(permuteFromLIDs.size()));
435  makeDualViewFromArrayView(this->TransferData_->remoteLIDs_,
436  remoteLIDs().getConst(),
437  "remoteLIDs");
438  TEUCHOS_ASSERT(size_t(this->TransferData_->remoteLIDs_.extent(0)) ==
439  size_t(remoteLIDs.size()));
440  makeDualViewFromArrayView(this->TransferData_->exportLIDs_,
441  exportLIDs().getConst(),
442  "exportLIDs");
443  TEUCHOS_ASSERT(size_t(this->TransferData_->exportLIDs_.extent(0)) ==
444  size_t(exportLIDs.size()));
445  this->TransferData_->exportPIDs_.swap(exportPIDs);
446  this->TransferData_->distributor_.swap(distributor);
447 
448  this->detectRemoteExportLIDsContiguous();
449 
450  TEUCHOS_ASSERT(!this->TransferData_->permuteFromLIDs_.need_sync_device());
451  TEUCHOS_ASSERT(!this->TransferData_->permuteFromLIDs_.need_sync_host());
452  TEUCHOS_ASSERT(!this->TransferData_->permuteToLIDs_.need_sync_device());
453  TEUCHOS_ASSERT(!this->TransferData_->permuteToLIDs_.need_sync_host());
454  TEUCHOS_ASSERT(!this->TransferData_->remoteLIDs_.need_sync_device());
455  TEUCHOS_ASSERT(!this->TransferData_->remoteLIDs_.need_sync_host());
456  TEUCHOS_ASSERT(!this->TransferData_->exportLIDs_.need_sync_device());
457  TEUCHOS_ASSERT(!this->TransferData_->exportLIDs_.need_sync_host());
458 }
459 
460 namespace { // (anonymous)
461 
462 template <class LO, class GO, class NT>
463 struct ImportLocalSetupResult {
464  Teuchos::RCP<const ::Tpetra::Map<LO, GO, NT>> targetMap;
465  LO numSameIDs;
466  // std::vector<LO> permuteToLIDs; // users aren't supposed to have permutes
467  // std::vector<LO> permuteFromLIDs; // users aren't suppoosed to have permutes
468  std::vector<GO> remoteGIDs;
469  std::vector<LO> remoteLIDs;
470  std::vector<int> remotePIDs;
471  LO numPermutes; // users aren't supposed to have permutes
472 };
473 
474 template <class T>
475 void printArray(std::ostream& out, const T x[], const std::size_t N) {
476  out << "[";
477  for (std::size_t k = 0; k < N; ++k) {
478  out << x[k];
479  if (k + 1 < N) {
480  out << ", ";
481  }
482  }
483  out << "]";
484 }
485 
486 template <class LO, class GO, class NT>
487 ImportLocalSetupResult<LO, GO, NT>
488 setupSamePermuteRemoteFromUserGlobalIndexList(const ::Tpetra::Map<LO, GO, NT>& sourceMap,
489  const GO targetMapRemoteOrPermuteGlobalIndices[],
490  const int targetMapRemoteOrPermuteProcessRanks[],
491  const LO numTargetMapRemoteOrPermuteGlobalIndices,
492  const bool mayReorderTargetMapIndicesLocally,
493  Teuchos::FancyOStream* out, // only valid if verbose
494  const std::string* verboseHeader, // only valid if verbose
495  const bool verbose,
496  const bool debug) {
497  using std::endl;
498  const int myRank = sourceMap.getComm()->getRank();
499  ImportLocalSetupResult<LO, GO, NT> result;
500 
501  if (verbose) {
502  std::ostringstream os;
503  os << *verboseHeader << "- Import::setupSPR w/ remote GIDs & PIDs: " << endl
504  << *verboseHeader << " Input GIDs: ";
505  printArray(os, targetMapRemoteOrPermuteGlobalIndices, numTargetMapRemoteOrPermuteGlobalIndices);
506  os << endl
507  << " Input PIDs: ";
508  printArray(os, targetMapRemoteOrPermuteProcessRanks, numTargetMapRemoteOrPermuteGlobalIndices);
509  os << endl;
510  *out << os.str();
511  }
512 
513  // In debug mode, check whether any of the input GIDs are
514  // actually in the source Map on the calling process. That's an
515  // error, because it means duplicate GIDs on the calling
516  // process. Also check if any of the input PIDs are invalid.
517  if (debug) {
518  std::vector<GO> badGIDs;
519  std::vector<int> badPIDs;
520  const Teuchos::Comm<int>& comm = *(sourceMap.getComm());
521  const int numProcs = comm.getSize();
522 
523  for (LO k = 0; k < numTargetMapRemoteOrPermuteGlobalIndices; ++k) {
524  const GO tgtGID = targetMapRemoteOrPermuteGlobalIndices[k];
525  if (sourceMap.isNodeGlobalElement(tgtGID)) {
526  badGIDs.push_back(tgtGID);
527  }
528  const int tgtPID = targetMapRemoteOrPermuteProcessRanks[k];
529  if (tgtPID < 0 || tgtPID >= numProcs) {
530  badPIDs.push_back(tgtPID);
531  }
532  }
533 
534  std::array<int, 2> lclStatus{{badGIDs.size() == 0 ? 1 : 0,
535  badPIDs.size() == 0 ? 1 : 0}};
536  std::array<int, 2> gblStatus{{0, 0}}; // output argument
537  Teuchos::reduceAll<int, int>(comm, Teuchos::REDUCE_MIN, 2,
538  lclStatus.data(), gblStatus.data());
539  const bool good = gblStatus[0] == 1 && gblStatus[1] == 1;
540  // Don't actually print all the "bad" GIDs and/or PIDs unless
541  // in verbose mode, since there could be many of them.
542  if (verbose && gblStatus[0] != 1) {
543  std::ostringstream os;
544  os << *verboseHeader << "- Some input GIDs are already in the source Map: ";
545  printArray(os, badGIDs.data(), badGIDs.size());
546  os << endl;
547  ::Tpetra::Details::gathervPrint(*out, os.str(), comm);
548  }
549  if (verbose && gblStatus[0] != 1) {
550  std::ostringstream os;
551  os << *verboseHeader << "- Some input PIDs are invalid: ";
552  printArray(os, badPIDs.data(), badPIDs.size());
553  os << endl;
554  ::Tpetra::Details::gathervPrint(*out, os.str(), comm);
555  }
556 
557  if (!good) {
558  std::ostringstream os;
559  os << "Tpetra::Import constructor that takes remote GIDs and PIDs: ";
560  if (gblStatus[0] != 1) {
561  os << "Some input GIDs (global indices) are already in the source Map! ";
562  }
563  if (gblStatus[1] != 1) {
564  os << "Some input PIDs (process ranks) are invalid! ";
565  }
566  os << "Rerun with the environment variable TPETRA_VERBOSE=Tpetra::Import "
567  "to see what GIDs and/or PIDs are bad.";
568  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
569  }
570  }
571 
572  // Create list of GIDs to go into target Map. We need to copy
573  // the GIDs into this list anyway, so once we have them, we can
574  // sort the "remotes" in place.
575  const LO numLclSrcIDs = static_cast<LO>(sourceMap.getLocalNumElements());
576  const LO numLclTgtIDs = numLclSrcIDs + numTargetMapRemoteOrPermuteGlobalIndices;
577  if (verbose) {
578  std::ostringstream os;
579  os << *verboseHeader << "- Copy source Map GIDs into target Map GID list: "
580  "numLclSrcIDs="
581  << numLclSrcIDs
582  << ", numTargetMapRemoteOrPermuteGlobalIndices="
583  << numTargetMapRemoteOrPermuteGlobalIndices << endl;
584  *out << os.str();
585  }
586  std::vector<GO> tgtGIDs(numLclTgtIDs); // will go into target Map ctor
587  if (sourceMap.isContiguous()) {
588  GO curTgtGID = sourceMap.getMinGlobalIndex();
589  for (LO k = 0; k < numLclSrcIDs; ++k, ++curTgtGID) {
590  tgtGIDs[k] = curTgtGID;
591  }
592  } else { // avoid calling getLocalElementList on a contiguous Map
593  auto srcGIDs = sourceMap.getLocalElementList(); // Teuchos::ArrayView has a different
594  for (LO k = 0; k < numLclSrcIDs; ++k) { // iterator type, so can't std::copy
595  tgtGIDs[k] = srcGIDs[k];
596  }
597  }
598  std::copy(targetMapRemoteOrPermuteGlobalIndices,
599  targetMapRemoteOrPermuteGlobalIndices + numTargetMapRemoteOrPermuteGlobalIndices,
600  tgtGIDs.begin() + numLclSrcIDs);
601 
602  // Optionally, sort input by process rank, so that remotes
603  // coming from the same process are grouped together. Only sort
604  // remote GIDs. While doing so, detect permutes (input "remote"
605  // GIDs whose rank is the same as that of the calling process).
606  //
607  // Permutes are actually an error. We normally detect them in
608  // debug mode, but if we sort, we have a nearly free opportunity
609  // to do so. We may also safely ignore permutes as duplicates.
610  //
611  // NOTE: tgtPIDs only includes remotes, not source Map entries.
612  if (verbose) {
613  std::ostringstream os;
614  os << *verboseHeader << "- Sort by PID? "
615  << (mayReorderTargetMapIndicesLocally ? "true" : "false") << endl;
616  *out << os.str();
617  }
618  std::vector<int> tgtPIDs(targetMapRemoteOrPermuteProcessRanks,
619  targetMapRemoteOrPermuteProcessRanks + numTargetMapRemoteOrPermuteGlobalIndices);
620  result.numPermutes = 0;
621  if (mayReorderTargetMapIndicesLocally) {
622  Tpetra::sort2(tgtPIDs.begin(), tgtPIDs.end(), tgtGIDs.begin() + numLclSrcIDs);
623  auto range = std::equal_range(tgtPIDs.begin(), tgtPIDs.end(), myRank); // binary search
624  if (range.second > range.first) {
625  result.numPermutes = static_cast<LO>(range.second - range.first);
626  }
627  } else { // don't sort; linear search to count permutes
628  result.numPermutes = static_cast<LO>(std::count(tgtPIDs.begin(), tgtPIDs.end(), myRank));
629  }
630  // The _actual_ number of remotes.
631  const LO numRemotes = numTargetMapRemoteOrPermuteGlobalIndices - result.numPermutes;
632  result.numSameIDs = static_cast<LO>(sourceMap.getLocalNumElements());
633 
634  if (verbose) {
635  std::ostringstream os;
636  os << *verboseHeader << "- numSame=" << result.numSameIDs
637  << ", numPermutes=" << result.numPermutes
638  << ", numRemotes=" << numRemotes << endl;
639  *out << os.str();
640  }
641 
642  if (result.numPermutes == 0) {
643  if (verbose) {
644  std::ostringstream os;
645  os << *verboseHeader << "- No permutes" << endl;
646  *out << os.str();
647  }
648  result.remoteGIDs = std::vector<GO>(tgtGIDs.begin() + numLclSrcIDs, tgtGIDs.end());
649  result.remotePIDs.swap(tgtPIDs);
650  result.remoteLIDs.resize(numRemotes);
651  for (LO k = 0; k < numRemotes; ++k) {
652  const LO tgtLid = result.numSameIDs + k;
653  result.remoteLIDs[k] = tgtLid;
654  }
655  if (verbose) {
656  std::ostringstream os;
657  os << *verboseHeader << "- Remote GIDs: "
658  << Teuchos::toString(result.remoteGIDs) << endl;
659  os << *verboseHeader << "- Remote PIDs: "
660  << Teuchos::toString(result.remotePIDs) << endl;
661  os << *verboseHeader << "- Remote LIDs: "
662  << Teuchos::toString(result.remoteLIDs) << endl;
663  *out << os.str();
664  }
665  } else { // separate permutes from remotes
666  // This case doesn't need to be optimal; it just needs to be
667  // correct. Users really shouldn't give permutes to this
668  // Import constructor.
669  result.remoteGIDs.reserve(numRemotes);
670  result.remoteLIDs.reserve(numRemotes);
671  result.remotePIDs.reserve(numRemotes);
672  for (LO k = 0; k < numTargetMapRemoteOrPermuteGlobalIndices; ++k) {
673  const LO tgtLid = result.numSameIDs + k;
674  const GO tgtGid = tgtGIDs[numLclSrcIDs + k];
675  const int tgtPid = tgtPIDs[k];
676 
677  if (tgtPid != myRank) { // it's a remote
678  result.remoteGIDs.push_back(tgtGid);
679  result.remoteLIDs.push_back(tgtLid);
680  result.remotePIDs.push_back(tgtPid);
681  }
682  }
683  if (verbose) {
684  std::ostringstream os;
685  os << *verboseHeader << "- Some permutes" << endl;
686  *out << os.str();
687  }
688  }
689 
690  if (sourceMap.isDistributed()) {
691  if (verbose) {
692  std::ostringstream os;
693  os << *verboseHeader << "- Sort remotes by PID, as Import always does"
694  << endl
695  << *verboseHeader << "-- remotePIDs before: "
696  << Teuchos::toString(result.remotePIDs) << endl
697  << *verboseHeader << "-- remoteGIDs before: "
698  << Teuchos::toString(result.remoteGIDs) << endl
699  << *verboseHeader << "-- remoteLIDs before: "
700  << Teuchos::toString(result.remoteLIDs) << endl;
701  *out << os.str();
702  }
703  // Import always sorts these, regardless of what the user wanted.
704  sort3(result.remotePIDs.begin(),
705  result.remotePIDs.end(),
706  result.remoteGIDs.begin(),
707  result.remoteLIDs.begin());
708  if (verbose) {
709  std::ostringstream os;
710  os << *verboseHeader << "-- remotePIDs after: "
711  << Teuchos::toString(result.remotePIDs) << endl
712  << *verboseHeader << "-- remoteGIDs after: "
713  << Teuchos::toString(result.remoteGIDs) << endl
714  << *verboseHeader << "-- remoteLIDs after: "
715  << Teuchos::toString(result.remoteLIDs) << endl;
716  std::cerr << os.str();
717  }
718  }
719 
720  if (verbose) {
721  std::ostringstream os;
722  os << *verboseHeader << "- Make target Map" << endl;
723  *out << os.str();
724  }
725  using ::Teuchos::rcp;
726  typedef ::Tpetra::Map<LO, GO, NT> map_type;
727  typedef ::Tpetra::global_size_t GST;
728  const GST MAP_COMPUTES_GLOBAL_COUNT = ::Teuchos::OrdinalTraits<GST>::invalid();
729  result.targetMap = rcp(new map_type(MAP_COMPUTES_GLOBAL_COUNT,
730  tgtGIDs.data(),
731  numLclTgtIDs,
732  sourceMap.getIndexBase(),
733  sourceMap.getComm()));
734  if (verbose) {
735  std::ostringstream os;
736  os << *verboseHeader << "- Done with sameSPR..." << endl;
737  *out << os.str();
738  }
739  return result;
740 }
741 } // namespace
742 
743 template <class LocalOrdinal, class GlobalOrdinal, class Node>
745  Import(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node>>& sourceMap,
746  const GlobalOrdinal targetMapRemoteOrPermuteGlobalIndices[],
747  const int targetMapRemoteOrPermuteProcessRanks[],
748  const LocalOrdinal numTargetMapRemoteOrPermuteGlobalIndices,
749  const bool mayReorderTargetMapIndicesLocally,
750  const Teuchos::RCP<Teuchos::ParameterList>& plist,
751  const Teuchos::RCP<Teuchos::FancyOStream>& debugOutput)
752  : // Special case: target Map is null on base_type construction.
753  // It's worthwhile for invariants like out_ not being null.
754  // We'll set TransferData_ again below.
755  base_type(sourceMap, Teuchos::null, debugOutput, plist, "Import") {
756  using std::endl;
757  using Teuchos::ArrayView;
758  using Teuchos::getFancyOStream;
759  using Teuchos::RCP;
760  using Teuchos::rcp;
761  using Teuchos::rcpFromRef;
762  using ::Tpetra::Details::Behavior;
763  using ::Tpetra::Details::makeDualViewFromOwningHostView;
764  using ::Tpetra::Details::makeDualViewFromVector;
765  using ::Tpetra::Details::printDualView;
766  using ::Tpetra::Details::view_alloc_no_init;
767  typedef LocalOrdinal LO;
768  typedef GlobalOrdinal GO;
769  typedef Node NT;
770 
771  const bool debug = Behavior::debug("Import") ||
772  Behavior::debug("Tpetra::Import");
773 
774  std::unique_ptr<std::string> verbPfx;
775  if (this->verbose()) {
776  std::ostringstream os;
777  const int myRank = sourceMap->getComm()->getRank();
778  os << "Proc " << myRank << ": Tpetra::Import ctor from remotes: ";
779  verbPfx = std::unique_ptr<std::string>(new std::string(os.str()));
780  os << "mayReorder=" << (mayReorderTargetMapIndicesLocally ? "true" : "false")
781  << endl;
782  this->verboseOutputStream() << os.str();
783  }
784 
785  TEUCHOS_ASSERT(!this->TransferData_.is_null());
786  ImportLocalSetupResult<LO, GO, NT> localSetupResult =
787  setupSamePermuteRemoteFromUserGlobalIndexList<LO, GO, NT>(*sourceMap,
788  targetMapRemoteOrPermuteGlobalIndices,
789  targetMapRemoteOrPermuteProcessRanks,
790  numTargetMapRemoteOrPermuteGlobalIndices,
791  mayReorderTargetMapIndicesLocally,
792  this->TransferData_->out_.getRawPtr(),
793  verbPfx.get(),
794  this->verbose(),
795  debug);
796 
797  // Since we invoked the base_type constructor above, we know that
798  // out_ is nonnull, so we don't have to waste time creating it
799  // again.
800  using data_type = ImportExportData<LO, GO, NT>;
801  TEUCHOS_ASSERT(!this->TransferData_.is_null());
802  this->TransferData_ = rcp(new data_type(sourceMap,
803  localSetupResult.targetMap,
804  this->TransferData_->out_,
805  plist));
806  this->TransferData_->numSameIDs_ = localSetupResult.numSameIDs;
807  // Skip permutes; they are user error, because they duplicate
808  // non-remote indices.
809  makeDualViewFromVector(this->TransferData_->remoteLIDs_,
810  localSetupResult.remoteLIDs,
811  "remoteLIDs");
812  // "Is locally complete" for an Import means that all target Map
813  // indices on the calling process exist on at least one process
814  // (not necessarily this one) in the source Map. For this
815  // constructor, this is true if and only if all input target PIDs
816  // are valid PIDs in the communicator.
817  //
818  // FIXME (mfh 20 Feb 2018) For now, assume this is always true.
819  this->TransferData_->isLocallyComplete_ = true;
820 
821  Teuchos::Array<GO> exportGIDs;
822  if (sourceMap->isDistributed()) {
823  if (this->verbose()) {
824  std::ostringstream os;
825  os << *verbPfx << "Make Distributor (createFromRecvs)" << endl;
826  this->verboseOutputStream() << os.str();
827  }
828  ArrayView<const GO> remoteGIDs(localSetupResult.remoteGIDs.data(),
829  localSetupResult.remoteGIDs.size());
830  ArrayView<const int> remotePIDs(localSetupResult.remotePIDs.data(),
831  localSetupResult.remotePIDs.size());
832  // Call Distributor::createFromRecvs to turn the remote GIDs and
833  // their owning PIDs into a send-and-receive communication plan.
834  // remoteGIDs and remotePIDs are input; exportGIDs and
835  // exportPIDs are output arrays that createFromRecvs allocates.
836  Distributor& distributor = this->TransferData_->distributor_;
837  distributor.createFromRecvs(remoteGIDs,
838  remotePIDs,
839  exportGIDs,
840  this->TransferData_->exportPIDs_);
841  // Find the LIDs corresponding to the (outgoing) GIDs in
842  // exportGIDs. For sparse matrix-vector multiply, this tells
843  // the calling process how to index into the source vector to
844  // get the elements which it needs to send.
845  //
846  // NOTE (mfh 03 Mar 2014) This is now a candidate for a
847  // thread-parallel kernel, but only if using the new thread-safe
848  // Map implementation.
849  if (this->verbose()) {
850  std::ostringstream os;
851  os << *verbPfx << "Compute exportLIDs" << endl;
852  this->verboseOutputStream() << os.str();
853  }
854  using size_type = typename Teuchos::Array<GO>::size_type;
855  const size_type numExportIDs = exportGIDs.size();
856 
857  typename decltype(this->TransferData_->exportLIDs_)::t_host
858  exportLIDs(view_alloc_no_init("exportLIDs"), numExportIDs);
859 
860  for (size_type k = 0; k < numExportIDs; ++k) {
861  exportLIDs[k] = sourceMap->getLocalElement(exportGIDs[k]);
862  }
863  makeDualViewFromOwningHostView(this->TransferData_->exportLIDs_, exportLIDs);
864  }
865 
866  if (this->verbose()) {
867  std::ostringstream os;
868  os << *verbPfx;
869  printDualView(os, this->TransferData_->remoteLIDs_,
870  "ImportExportData::remoteLIDs_");
871  os << endl;
872  this->verboseOutputStream() << os.str();
873  }
874  if (this->verbose()) {
875  std::ostringstream os;
876  os << *verbPfx << "Done!" << endl;
877  this->verboseOutputStream() << os.str();
878  }
879 }
880 
881 template <class LocalOrdinal, class GlobalOrdinal, class Node>
883  describe(Teuchos::FancyOStream& out,
884  const Teuchos::EVerbosityLevel verbLevel) const {
885  // Call the base class' method. It does all the work.
886  this->describeImpl(out, "Tpetra::Import", verbLevel);
887 }
888 
889 template <class LocalOrdinal, class GlobalOrdinal, class Node>
891  print(std::ostream& os) const {
892  auto out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(os));
893  // "Print" traditionally meant "everything."
894  this->describe(*out, Teuchos::VERB_EXTREME);
895 }
896 
897 template <class LocalOrdinal, class GlobalOrdinal, class Node>
899  setupSamePermuteRemote(Teuchos::Array<GlobalOrdinal>& remoteGIDs) {
900  using Teuchos::arcp;
901  using Teuchos::Array;
902  using Teuchos::ArrayRCP;
903  using Teuchos::ArrayView;
904  using Teuchos::as;
905  using Teuchos::null;
906  using ::Tpetra::Details::makeDualViewFromOwningHostView;
907  using ::Tpetra::Details::ProfilingRegion;
908  using ::Tpetra::Details::view_alloc_no_init;
909  typedef LocalOrdinal LO;
910  typedef GlobalOrdinal GO;
911  typedef typename ArrayView<const GO>::size_type size_type;
912  ProfilingRegion regionExport("Tpetra::Import::setupSamePermuteRemote");
913 
914  const map_type& source = *(this->getSourceMap());
915  const map_type& target = *(this->getTargetMap());
916  ArrayView<const GO> sourceGIDs = source.getLocalElementList();
917  ArrayView<const GO> targetGIDs = target.getLocalElementList();
918 
919 #ifdef HAVE_TPETRA_DEBUG
920  ArrayView<const GO> rawSrcGids = sourceGIDs;
921  ArrayView<const GO> rawTgtGids = targetGIDs;
922 #else
923  const GO* const rawSrcGids = sourceGIDs.getRawPtr();
924  const GO* const rawTgtGids = targetGIDs.getRawPtr();
925 #endif // HAVE_TPETRA_DEBUG
926  const size_type numSrcGids = sourceGIDs.size();
927  const size_type numTgtGids = targetGIDs.size();
928  const size_type numGids = std::min(numSrcGids, numTgtGids);
929 
930  // Compute numSameIDs_: the number of initial GIDs that are the
931  // same (and occur in the same order) in both Maps. The point of
932  // numSameIDs_ is for the common case of an Import where all the
933  // overlapping GIDs are at the end of the target Map, but
934  // otherwise the source and target Maps are the same. This allows
935  // a fast contiguous copy for the initial "same IDs."
936  size_type numSameGids = 0;
937  for (; numSameGids < numGids && rawSrcGids[numSameGids] == rawTgtGids[numSameGids]; ++numSameGids) {
938  } // third clause of 'for' does everything
939  this->TransferData_->numSameIDs_ = numSameGids;
940 
941  // Compute permuteToLIDs_, permuteFromLIDs_, remoteGIDs, and
942  // remoteLIDs_. The first two arrays are IDs to be permuted, and
943  // the latter two arrays are IDs to be received ("imported"),
944  // called "remote" IDs.
945  //
946  // IDs to permute are in both the source and target Maps, which
947  // means we don't have to send or receive them, but we do have to
948  // rearrange (permute) them in general. IDs to receive are in the
949  // target Map, but not the source Map.
950 
951  // Iterate over the target Map's LIDs, since we only need to do
952  // GID -> LID lookups for the source Map.
953  const LO LINVALID = Teuchos::OrdinalTraits<LO>::invalid();
954  const LO numTgtLids = as<LO>(numTgtGids);
955  LO numPermutes = 0;
956 
957  for (LO tgtLid = numSameGids; tgtLid < numTgtLids; ++tgtLid) {
958  const GO curTargetGid = rawTgtGids[tgtLid];
959  // getLocalElement() returns LINVALID if the GID isn't in the
960  // source Map. This saves us a lookup (which
961  // isNodeGlobalElement() would do).
962  const LO srcLid = source.getLocalElement(curTargetGid);
963  if (srcLid != LINVALID) { // if source.isNodeGlobalElement (curTargetGid)
964  ++numPermutes;
965  }
966  }
967  const LO numRemotes = (numTgtLids - numSameGids) - numPermutes;
968 
969  using host_perm_type =
970  typename decltype(this->TransferData_->permuteToLIDs_)::t_host;
971  host_perm_type permuteToLIDs(view_alloc_no_init("permuteToLIDs"), numPermutes);
972  host_perm_type permuteFromLIDs(view_alloc_no_init("permuteFromLIDs"), numPermutes);
973  typename decltype(this->TransferData_->remoteLIDs_)::t_host remoteLIDs(view_alloc_no_init("permuteFromLIDs"), numRemotes);
974 
975  {
976  LO numPermutes2 = 0;
977  LO numRemotes2 = 0;
978  for (LO tgtLid = numSameGids; tgtLid < numTgtLids; ++tgtLid) {
979  const GO curTargetGid = rawTgtGids[tgtLid];
980  const LO srcLid = source.getLocalElement(curTargetGid);
981  if (srcLid != LINVALID) {
982  permuteToLIDs[numPermutes2] = tgtLid;
983  permuteFromLIDs[numPermutes2] = srcLid;
984  ++numPermutes2;
985  } else {
986  remoteGIDs.push_back(curTargetGid);
987  remoteLIDs[numRemotes2] = tgtLid;
988  ++numRemotes2;
989  }
990  }
991  TEUCHOS_ASSERT(numPermutes == numPermutes2);
992  TEUCHOS_ASSERT(numRemotes == numRemotes2);
993  TEUCHOS_ASSERT(size_t(numPermutes) + remoteGIDs.size() == size_t(numTgtLids - numSameGids));
994  }
995 
996  makeDualViewFromOwningHostView(this->TransferData_->permuteToLIDs_, permuteToLIDs);
997  makeDualViewFromOwningHostView(this->TransferData_->permuteFromLIDs_, permuteFromLIDs);
998  makeDualViewFromOwningHostView(this->TransferData_->remoteLIDs_, remoteLIDs);
999  if (remoteLIDs.extent(0) != 0 && !source.isDistributed()) {
1000  // This Import has remote LIDs, meaning that the target Map has
1001  // entries on this process that are not in the source Map on
1002  // this process. However, the source Map is not distributed
1003  // globally. This implies that this Import is not locally
1004  // complete on this process.
1005  this->TransferData_->isLocallyComplete_ = false;
1006  // mfh 12 Sep 2016: I disagree that this is "abuse"; it may be
1007  // correct behavior, depending on the circumstances.
1008  TPETRA_ABUSE_WARNING(true, std::runtime_error,
1009  "::setupSamePermuteRemote(): Target has "
1010  "remote LIDs but Source is not distributed globally. Importing to a "
1011  "submap of the target map.");
1012  }
1013 }
1014 
1015 template <class LocalOrdinal, class GlobalOrdinal, class Node>
1016 void Import<LocalOrdinal, GlobalOrdinal, Node>::
1017  setupExport(Teuchos::Array<GlobalOrdinal>& remoteGIDs,
1018  bool useRemotePIDs,
1019  Teuchos::Array<int>& userRemotePIDs,
1020  const Teuchos::RCP<Teuchos::ParameterList>& plist) {
1021  using std::endl;
1022  using Teuchos::Array;
1023  using Teuchos::ArrayView;
1024  using ::Tpetra::Details::makeDualViewFromOwningHostView;
1025  using ::Tpetra::Details::view_alloc_no_init;
1026  using GO = GlobalOrdinal;
1027  typedef typename Array<int>::difference_type size_type;
1028  const char tfecfFuncName[] = "setupExport: ";
1029  const char suffix[] = " Please report this bug to the Tpetra developers.";
1030 
1031  std::unique_ptr<std::string> prefix;
1032  if (this->verbose()) {
1033  auto srcMap = this->getSourceMap();
1034  auto comm = srcMap.is_null() ? Teuchos::null : srcMap->getComm();
1035  const int myRank = comm.is_null() ? -1 : comm->getRank();
1036  std::ostringstream os;
1037  os << "Proc " << myRank << ": Tpetra::Import::setupExport: ";
1038  prefix = std::unique_ptr<std::string>(new std::string(os.str()));
1039  os << "Start" << std::endl;
1040  this->verboseOutputStream() << os.str();
1041  }
1042 
1043  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(this->getSourceMap().is_null(), std::logic_error,
1044  "Source Map is null. " << suffix);
1045  const map_type& source = *(this->getSourceMap());
1046 
1047  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!useRemotePIDs && (userRemotePIDs.size() > 0), std::invalid_argument,
1048  "remotePIDs are non-empty but their use has not been requested.");
1049  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(userRemotePIDs.size() > 0 && remoteGIDs.size() != userRemotePIDs.size(),
1050  std::invalid_argument,
1051  "remotePIDs must either be of size zero or match "
1052  "the size of remoteGIDs.");
1053 
1054  // For each entry remoteGIDs[i], remoteProcIDs[i] will contain
1055  // the process ID of the process that owns that GID.
1056  ArrayView<GO> remoteGIDsView = remoteGIDs();
1057  ArrayView<int> remoteProcIDsView;
1058 
1059  // lookup == IDNotPresent means that the source Map wasn't able to
1060  // figure out to which processes one or more of the GIDs in the
1061  // given list of remote GIDs belong.
1062  //
1063  // The previous abuse warning said "The target Map has GIDs not
1064  // found in the source Map." This statement could be confusing,
1065  // because it doesn't refer to ownership by the current process,
1066  // but rather to ownership by _any_ process participating in the
1067  // Map. (It could not possibly refer to ownership by the current
1068  // process, since remoteGIDs is exactly the list of GIDs owned by
1069  // the target Map but not owned by the source Map. It was
1070  // constructed that way by setupSamePermuteRemote().)
1071  //
1072  // What this statement means is that the source and target Maps
1073  // don't contain the same set of GIDs globally (over all
1074  // processes). That is, there is at least one GID owned by some
1075  // process in the target Map, which is not owned by _any_ process
1076  // in the source Map.
1077  Array<int> newRemotePIDs;
1078  LookupStatus lookup = AllIDsPresent;
1079 
1080  if (!useRemotePIDs) {
1081  newRemotePIDs.resize(remoteGIDsView.size());
1082  if (this->verbose()) {
1083  std::ostringstream os;
1084  os << *prefix << "Call sourceMap.getRemoteIndexList" << endl;
1085  this->verboseOutputStream() << os.str();
1086  }
1087  lookup = source.getRemoteIndexList(remoteGIDsView, newRemotePIDs());
1088  }
1089  Array<int>& remoteProcIDs = useRemotePIDs ? userRemotePIDs : newRemotePIDs;
1090 
1091  if (lookup == IDNotPresent) {
1092  // There is at least one GID owned by the calling process in the
1093  // target Map, which is not owned by any process in the source
1094  // Map.
1095  this->TransferData_->isLocallyComplete_ = false;
1096 
1097  // mfh 12 Sep 2016: I disagree that this is "abuse"; it may be
1098  // correct behavior, depending on the circumstances.
1099  TPETRA_ABUSE_WARNING(true, std::runtime_error,
1100  "::setupExport(): the source Map wasn't "
1101  "able to figure out which process owns one or more of the GIDs in the "
1102  "list of remote GIDs. This probably means that there is at least one "
1103  "GID owned by some process in the target Map which is not owned by any"
1104  " process in the source Map. (That is, the source and target Maps do "
1105  "not contain the same set of GIDs globally.)");
1106 
1107  // Ignore remote GIDs that aren't owned by any process in the
1108  // source Map. getRemoteIndexList() gives each of these a
1109  // process ID of -1.
1110 
1111  const size_type numInvalidRemote =
1112  std::count_if(remoteProcIDs.begin(), remoteProcIDs.end(),
1113  std::bind(std::equal_to<int>(), -1, std::placeholders::_1));
1114  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numInvalidRemote == 0, std::logic_error,
1115  "Calling getRemoteIndexList "
1116  "on the source Map returned IDNotPresent, but none of the returned "
1117  "\"remote\" process ranks are -1. Please report this bug to the "
1118  "Tpetra developers.");
1119 
1120 #ifdef HAVE_TPETRA_MMM_TIMINGS
1121  using Teuchos::TimeMonitor;
1122  std::string label;
1123  if (!plist.is_null())
1124  label = plist->get("Timer Label", label);
1125  std::string prefix = std::string("Tpetra ") + label + std::string(":iport_ctor:setupExport:1 ");
1126  auto MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix)));
1127 #else
1128  (void)plist;
1129 #endif
1130 
1131  // If all of them are invalid, we can delete the whole array.
1132  const size_type totalNumRemote = this->getNumRemoteIDs();
1133  if (numInvalidRemote == totalNumRemote) {
1134  // all remotes are invalid; we have no remotes; we can delete the remotes
1135  remoteProcIDs.clear();
1136  remoteGIDs.clear(); // invalidates remoteGIDsView
1137  this->TransferData_->remoteLIDs_ =
1138  decltype(this->TransferData_->remoteLIDs_)();
1139  } else {
1140  // Some remotes are valid; we need to keep the valid ones.
1141  // Pack and resize remoteProcIDs, remoteGIDs, and remoteLIDs_.
1142  size_type numValidRemote = 0;
1143 #ifdef HAVE_TPETRA_DEBUG
1144  ArrayView<GO> remoteGIDsPtr = remoteGIDsView;
1145 #else
1146  GO* const remoteGIDsPtr = remoteGIDsView.getRawPtr();
1147 #endif // HAVE_TPETRA_DEBUG
1148 
1149  // Don't mark the DualView modified, since we'll reallocate it.
1150  auto remoteLIDs = this->TransferData_->remoteLIDs_.view_host();
1151 
1152  for (size_type r = 0; r < totalNumRemote; ++r) {
1153  // Pack in all the valid remote PIDs and GIDs.
1154  if (remoteProcIDs[r] != -1) {
1155  remoteProcIDs[numValidRemote] = remoteProcIDs[r];
1156  remoteGIDsPtr[numValidRemote] = remoteGIDsPtr[r];
1157  remoteLIDs[numValidRemote] = remoteLIDs[r];
1158  ++numValidRemote;
1159  }
1160  }
1161  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numValidRemote != totalNumRemote - numInvalidRemote,
1162  std::logic_error,
1163  "After removing invalid remote GIDs and packing "
1164  "the valid remote GIDs, numValidRemote = "
1165  << numValidRemote
1166  << " != totalNumRemote - numInvalidRemote = "
1167  << totalNumRemote - numInvalidRemote
1168  << ". Please report this bug to the Tpetra developers.");
1169 
1170  remoteProcIDs.resize(numValidRemote);
1171  remoteGIDs.resize(numValidRemote);
1172 
1173  Kokkos::resize(remoteLIDs, numValidRemote);
1174  this->TransferData_->remoteLIDs_ = decltype(this->TransferData_->remoteLIDs_)();
1175  makeDualViewFromOwningHostView(this->TransferData_->remoteLIDs_, remoteLIDs);
1176  }
1177  // Revalidate the view after clear or resize.
1178  remoteGIDsView = remoteGIDs();
1179  }
1180 
1181  // Sort remoteProcIDs in ascending order, and apply the resulting
1182  // permutation to remoteGIDs and remoteLIDs_. This ensures that
1183  // remoteProcIDs[i], remoteGIDs[i], and remoteLIDs_[i] all refer
1184  // to the same thing.
1185  {
1186  this->TransferData_->remoteLIDs_.modify_host();
1187  auto remoteLIDs = this->TransferData_->remoteLIDs_.view_host();
1188  sort3(remoteProcIDs.begin(),
1189  remoteProcIDs.end(),
1190  remoteGIDsView.getRawPtr(),
1191  remoteLIDs.data());
1192  this->TransferData_->remoteLIDs_.sync_device();
1193  }
1194 
1195  // Call the Distributor's createFromRecvs() method to turn the
1196  // remote GIDs and their owning processes into a send-and-receive
1197  // communication plan. remoteGIDs and remoteProcIDs_ are input;
1198  // exportGIDs and exportProcIDs_ are output arrays which are
1199  // allocated by createFromRecvs().
1200  Array<GO> exportGIDs;
1201 
1202 #ifdef HAVE_TPETRA_MMM_TIMINGS
1203  using Teuchos::TimeMonitor;
1204  std::string label;
1205  if (!plist.is_null())
1206  label = plist->get("Timer Label", label);
1207  std::string prefix2 = std::string("Tpetra ") + label + std::string(":iport_ctor:setupExport:3 ");
1208  auto MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix2)));
1209 #endif
1210 
1211  if (this->verbose()) {
1212  std::ostringstream os;
1213  os << *prefix << "Call createFromRecvs" << endl;
1214  this->verboseOutputStream() << endl;
1215  }
1216  this->TransferData_->distributor_.createFromRecvs(remoteGIDsView().getConst(),
1217  remoteProcIDs, exportGIDs,
1218  this->TransferData_->exportPIDs_);
1219 
1220  // Find the LIDs corresponding to the (outgoing) GIDs in
1221  // exportGIDs. For sparse matrix-vector multiply, this tells the
1222  // calling process how to index into the source vector to get the
1223  // elements which it needs to send.
1224 #ifdef HAVE_TPETRA_MMM_TIMINGS
1225  prefix2 = std::string("Tpetra ") + label + std::string(":iport_ctor:setupExport:4 ");
1226  MM = Teuchos::null;
1227  MM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(prefix2)));
1228 #endif
1229 
1230  // NOTE (mfh 03 Mar 2014) This is now a candidate for a
1231  // thread-parallel kernel, but only if using the new thread-safe
1232  // Map implementation.
1233  const size_type numExportIDs = exportGIDs.size();
1234  if (numExportIDs > 0) {
1235  typename decltype(this->TransferData_->exportLIDs_)::t_host
1236  exportLIDs(view_alloc_no_init("exportLIDs"), numExportIDs);
1237  ArrayView<const GO> expGIDs = exportGIDs();
1238 
1239  for (size_type k = 0; k < numExportIDs; ++k) {
1240  exportLIDs[k] = source.getLocalElement(expGIDs[k]);
1241  }
1242  makeDualViewFromOwningHostView(this->TransferData_->exportLIDs_, exportLIDs);
1243  }
1244 
1245  if (this->verbose()) {
1246  std::ostringstream os;
1247  os << *prefix << "Done!" << endl;
1248  this->verboseOutputStream() << os.str();
1249  }
1250 }
1251 
1252 template <class LocalOrdinal, class GlobalOrdinal, class Node>
1254  findUnionTargetGIDs(Teuchos::Array<GlobalOrdinal>& unionTgtGIDs,
1255  Teuchos::Array<std::pair<int, GlobalOrdinal>>& remotePGIDs,
1256  typename Teuchos::Array<GlobalOrdinal>::size_type& numSameGIDs,
1257  typename Teuchos::Array<GlobalOrdinal>::size_type& numPermuteGIDs,
1258  typename Teuchos::Array<GlobalOrdinal>::size_type& numRemoteGIDs,
1259  const Teuchos::ArrayView<const GlobalOrdinal>& sameGIDs1,
1260  const Teuchos::ArrayView<const GlobalOrdinal>& sameGIDs2,
1261  Teuchos::Array<GlobalOrdinal>& permuteGIDs1,
1262  Teuchos::Array<GlobalOrdinal>& permuteGIDs2,
1263  Teuchos::Array<GlobalOrdinal>& remoteGIDs1,
1264  Teuchos::Array<GlobalOrdinal>& remoteGIDs2,
1265  Teuchos::Array<int>& remotePIDs1,
1266  Teuchos::Array<int>& remotePIDs2) const {
1267  typedef GlobalOrdinal GO;
1268  typedef typename Teuchos::Array<GO>::size_type size_type;
1269 
1270  const size_type numSameGIDs1 = sameGIDs1.size();
1271  const size_type numSameGIDs2 = sameGIDs2.size();
1272 
1273  // Sort the permute GIDs
1274  std::sort(permuteGIDs1.begin(), permuteGIDs1.end());
1275  std::sort(permuteGIDs2.begin(), permuteGIDs2.end());
1276 
1277  // Get the union of the two target maps
1278  // Reserve the maximum possible size to guard against reallocations from
1279  // push_back operations.
1280  unionTgtGIDs.reserve(numSameGIDs1 + numSameGIDs2 +
1281  permuteGIDs1.size() + permuteGIDs2.size() +
1282  remoteGIDs1.size() + remoteGIDs2.size());
1283 
1284  // Copy the same GIDs to unionTgtGIDs. Cases for numSameGIDs1 !=
1285  // numSameGIDs2 must be treated separately.
1286  typename Teuchos::Array<GO>::iterator permuteGIDs1_end;
1287  typename Teuchos::Array<GO>::iterator permuteGIDs2_end;
1288  if (numSameGIDs2 > numSameGIDs1) {
1289  numSameGIDs = numSameGIDs2;
1290  permuteGIDs2_end = permuteGIDs2.end();
1291 
1292  // Copy the same GIDs from tgtGIDs to the union
1293  std::copy(sameGIDs2.begin(), sameGIDs2.end(), std::back_inserter(unionTgtGIDs));
1294 
1295  // Remove GIDs from permuteGIDs1 that have already been copied in to unionTgtGIDs
1296  // set_difference allows the last (output) argument to alias the first.
1297  permuteGIDs1_end = std::set_difference(permuteGIDs1.begin(), permuteGIDs1.end(),
1298  unionTgtGIDs.begin() + numSameGIDs1, unionTgtGIDs.end(),
1299  permuteGIDs1.begin());
1300 
1301  } else {
1302  numSameGIDs = numSameGIDs1;
1303  permuteGIDs1_end = permuteGIDs1.end();
1304 
1305  // Copy the same GIDs from tgtGIDs to the union
1306  std::copy(sameGIDs1.begin(), sameGIDs1.end(), std::back_inserter(unionTgtGIDs));
1307 
1308  // Remove GIDs from permuteGIDs2 that have already been copied in to unionTgtGIDs
1309  // set_difference allows the last (output) argument to alias the first.
1310  permuteGIDs2_end = std::set_difference(permuteGIDs2.begin(), permuteGIDs2.end(),
1311  unionTgtGIDs.begin() + numSameGIDs2, unionTgtGIDs.end(),
1312  permuteGIDs2.begin());
1313  }
1314 
1315  // Get the union of the permute GIDs and push it back on unionTgtGIDs
1316  std::set_union(permuteGIDs1.begin(), permuteGIDs1_end,
1317  permuteGIDs2.begin(), permuteGIDs2_end,
1318  std::back_inserter(unionTgtGIDs));
1319 
1320  // Sort the PID,GID pairs and find the unique set
1321  Teuchos::Array<std::pair<int, GO>> remotePGIDs1(remoteGIDs1.size());
1322  for (size_type k = 0; k < remoteGIDs1.size(); k++)
1323  remotePGIDs1[k] = std::make_pair(remotePIDs1[k], remoteGIDs1[k]);
1324  std::sort(remotePGIDs1.begin(), remotePGIDs1.end());
1325 
1326  Teuchos::Array<std::pair<int, GO>> remotePGIDs2(remoteGIDs2.size());
1327  for (size_type k = 0; k < remoteGIDs2.size(); k++)
1328  remotePGIDs2[k] = std::make_pair(remotePIDs2[k], remoteGIDs2[k]);
1329  std::sort(remotePGIDs2.begin(), remotePGIDs2.end());
1330 
1331  remotePGIDs.reserve(remotePGIDs1.size() + remotePGIDs2.size());
1332  std::merge(remotePGIDs1.begin(), remotePGIDs1.end(),
1333  remotePGIDs2.begin(), remotePGIDs2.end(),
1334  std::back_inserter(remotePGIDs));
1335  auto it = std::unique(remotePGIDs.begin(), remotePGIDs.end());
1336  remotePGIDs.resize(std::distance(remotePGIDs.begin(), it));
1337 
1338  // Finally, insert remote GIDs
1339  const size_type oldSize = unionTgtGIDs.size();
1340  unionTgtGIDs.resize(oldSize + remotePGIDs.size());
1341  for (size_type start = oldSize, k = 0; k < remotePGIDs.size(); k++)
1342  unionTgtGIDs[start + k] = remotePGIDs[k].second;
1343 
1344  // Compute output only quantities
1345  numRemoteGIDs = remotePGIDs.size();
1346  numPermuteGIDs = unionTgtGIDs.size() - numSameGIDs - numRemoteGIDs;
1347 
1348  return;
1349 }
1350 
1351 template <class LocalOrdinal, class GlobalOrdinal, class Node>
1352 Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node>>
1355  using Teuchos::Array;
1356  using Teuchos::ArrayView;
1357  using Teuchos::as;
1358  using Teuchos::Comm;
1359  using Teuchos::outArg;
1360  using Teuchos::RCP;
1361  using Teuchos::rcp;
1362  using Teuchos::REDUCE_MIN;
1363  using Teuchos::reduceAll;
1364  using ::Tpetra::Details::Behavior;
1365  using GST = Tpetra::global_size_t;
1366  using LO = LocalOrdinal;
1367  using GO = GlobalOrdinal;
1368  using import_type = Import<LO, GO, Node>;
1369  using size_type = typename Array<GO>::size_type;
1370 
1371 #ifdef HAVE_TPETRA_MMM_TIMINGS
1372  using Teuchos::TimeMonitor;
1373  std::string label = std::string("Tpetra::Import::setUnion");
1374  TimeMonitor MM(*TimeMonitor::getNewTimer(label));
1375 #endif
1376 
1377  RCP<const map_type> srcMap = this->getSourceMap();
1378  RCP<const map_type> tgtMap1 = this->getTargetMap();
1379  RCP<const map_type> tgtMap2 = rhs.getTargetMap();
1380  RCP<const Comm<int>> comm = srcMap->getComm();
1381 
1382  const bool debug = Behavior::debug("Import::setUnion") ||
1383  Behavior::debug("Tpetra::Import::setUnion");
1384 
1385  if (debug) {
1386  TEUCHOS_TEST_FOR_EXCEPTION(!srcMap->isSameAs(*(rhs.getSourceMap())), std::invalid_argument,
1387  "Tpetra::Import::setUnion: The source Map of the input Import must be the "
1388  "same as (in the sense of Map::isSameAs) the source Map of this Import.");
1389  const Comm<int>& comm1 = *(tgtMap1->getComm());
1390  const Comm<int>& comm2 = *(tgtMap2->getComm());
1391  TEUCHOS_TEST_FOR_EXCEPTION(!::Tpetra::Details::congruent(comm1, comm2),
1392  std::invalid_argument,
1393  "Tpetra::Import::setUnion: "
1394  "The target Maps must have congruent communicators.");
1395  }
1396 
1397  // It's probably worth the one all-reduce to check whether the two
1398  // Maps are the same. If so, we can just return a copy of *this.
1399  // isSameAs() bypasses the all-reduce if the pointers are equal.
1400  if (tgtMap1->isSameAs(*tgtMap2)) {
1401  return rcp(new import_type(*this));
1402  }
1403 
1404  // Alas, the two target Maps are not the same. That means we have
1405  // to compute their union, and the union Import object.
1406 
1407  // Get the same GIDs (same GIDs are a subview of the first numSame target
1408  // GIDs)
1409  const size_type numSameGIDs1 = this->getNumSameIDs();
1410  ArrayView<const GO> sameGIDs1 = (tgtMap1->getLocalElementList())(0, numSameGIDs1);
1411 
1412  const size_type numSameGIDs2 = rhs.getNumSameIDs();
1413  ArrayView<const GO> sameGIDs2 = (tgtMap2->getLocalElementList())(0, numSameGIDs2);
1414 
1415  // Get permute GIDs
1416  ArrayView<const LO> permuteToLIDs1 = this->getPermuteToLIDs();
1417  Array<GO> permuteGIDs1(permuteToLIDs1.size());
1418  for (size_type k = 0; k < permuteGIDs1.size(); k++)
1419  permuteGIDs1[k] = tgtMap1->getGlobalElement(permuteToLIDs1[k]);
1420 
1421  ArrayView<const LO> permuteToLIDs2 = rhs.getPermuteToLIDs();
1422  Array<GO> permuteGIDs2(permuteToLIDs2.size());
1423  for (size_type k = 0; k < permuteGIDs2.size(); k++)
1424  permuteGIDs2[k] = tgtMap2->getGlobalElement(permuteToLIDs2[k]);
1425 
1426  // Get remote GIDs
1427  ArrayView<const LO> remoteLIDs1 = this->getRemoteLIDs();
1428  Array<GO> remoteGIDs1(remoteLIDs1.size());
1429  for (size_type k = 0; k < remoteLIDs1.size(); k++)
1430  remoteGIDs1[k] = this->getTargetMap()->getGlobalElement(remoteLIDs1[k]);
1431 
1432  ArrayView<const LO> remoteLIDs2 = rhs.getRemoteLIDs();
1433  Array<GO> remoteGIDs2(remoteLIDs2.size());
1434  for (size_type k = 0; k < remoteLIDs2.size(); k++)
1435  remoteGIDs2[k] = rhs.getTargetMap()->getGlobalElement(remoteLIDs2[k]);
1436 
1437  // Get remote PIDs
1438  Array<int> remotePIDs1;
1439  Tpetra::Import_Util::getRemotePIDs(*this, remotePIDs1);
1440 
1441  Array<int> remotePIDs2;
1442  Tpetra::Import_Util::getRemotePIDs(rhs, remotePIDs2);
1443 
1444  // Get the union of the target GIDs
1445  Array<GO> unionTgtGIDs;
1446  Array<std::pair<int, GO>> remotePGIDs;
1447  size_type numSameIDsUnion, numPermuteIDsUnion, numRemoteIDsUnion;
1448 
1449  findUnionTargetGIDs(unionTgtGIDs, remotePGIDs,
1450  numSameIDsUnion, numPermuteIDsUnion, numRemoteIDsUnion,
1451  sameGIDs1, sameGIDs2, permuteGIDs1, permuteGIDs2,
1452  remoteGIDs1, remoteGIDs2, remotePIDs1, remotePIDs2);
1453 
1454  // Extract GIDs and compute LIDS, PIDs for the remotes in the union
1455  Array<LO> remoteLIDsUnion(numRemoteIDsUnion);
1456  Array<GO> remoteGIDsUnion(numRemoteIDsUnion);
1457  Array<int> remotePIDsUnion(numRemoteIDsUnion);
1458  const size_type unionRemoteIDsStart = numSameIDsUnion + numPermuteIDsUnion;
1459  for (size_type k = 0; k < numRemoteIDsUnion; ++k) {
1460  remoteLIDsUnion[k] = unionRemoteIDsStart + k;
1461  remotePIDsUnion[k] = remotePGIDs[k].first;
1462  remoteGIDsUnion[k] = remotePGIDs[k].second;
1463  }
1464 
1465  // Compute the permute-to LIDs (in the union target Map).
1466  // Convert the permute GIDs to permute-from LIDs in the source Map.
1467  Array<LO> permuteToLIDsUnion(numPermuteIDsUnion);
1468  Array<LO> permuteFromLIDsUnion(numPermuteIDsUnion);
1469 
1470  for (size_type k = 0; k < numPermuteIDsUnion; ++k) {
1471  size_type idx = numSameIDsUnion + k;
1472  permuteToLIDsUnion[k] = static_cast<LO>(idx);
1473  permuteFromLIDsUnion[k] = srcMap->getLocalElement(unionTgtGIDs[idx]);
1474  }
1475 
1476 #ifdef HAVE_TPETRA_MMM_TIMINGS
1477  MM.disableTimer(label);
1478  label = "Tpetra::Import::setUnion : Construct Target Map";
1479  TimeMonitor MM2(*TimeMonitor::getNewTimer(label));
1480 #endif
1481 
1482  // Create the union target Map.
1483  const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid();
1484  const GO indexBaseUnion = std::min(tgtMap1->getIndexBase(), tgtMap2->getIndexBase());
1485  RCP<const map_type> unionTgtMap =
1486  rcp(new map_type(INVALID, unionTgtGIDs(), indexBaseUnion, comm));
1487 
1488 #ifdef HAVE_TPETRA_MMM_TIMINGS
1489  MM2.disableTimer(label);
1490  label = "Tpetra::Import::setUnion : Export GIDs";
1491  TimeMonitor MM3(*TimeMonitor::getNewTimer(label));
1492 #endif
1493 
1494  // Thus far, we have computed the following in the union Import:
1495  // - numSameIDs
1496  // - numPermuteIDs and permuteFromLIDs
1497  // - numRemoteIDs, remoteGIDs, remoteLIDs, and remotePIDs
1498  //
1499  // Now it's time to compute the export IDs and initialize the
1500  // Distributor.
1501 
1502  Array<GO> exportGIDsUnion;
1503  Array<LO> exportLIDsUnion;
1504  Array<int> exportPIDsUnion;
1505  Distributor distributor(comm, this->TransferData_->out_);
1506 
1507 #ifdef TPETRA_IMPORT_SETUNION_USE_CREATE_FROM_SENDS
1508  // Compute the export IDs without communication, by merging the
1509  // lists of (export LID, export PID) pairs from the two input
1510  // Import objects. The export LIDs in both input Import objects
1511  // are LIDs in the source Map. Then, use the export PIDs to
1512  // initialize the Distributor via createFromSends.
1513 
1514  // const size_type numExportIDs1 = this->getNumExportIDs ();
1515  ArrayView<const LO> exportLIDs1 = this->getExportLIDs();
1516  ArrayView<const LO> exportPIDs1 = this->getExportPIDs();
1517 
1518  // const size_type numExportIDs2 = rhs.getNumExportIDs ();
1519  ArrayView<const LO> exportLIDs2 = rhs.getExportLIDs();
1520  ArrayView<const LO> exportPIDs2 = rhs.getExportPIDs();
1521 
1522  // We have to keep the export LIDs in PID-sorted order, then merge
1523  // them. So, first key-value merge (LID,PID) pairs, treating PIDs
1524  // as values, merging values by replacement. Then, sort the
1525  // (LID,PID) pairs again by PID.
1526 
1527  // Sort (LID,PID) pairs by LID for the later merge, and make
1528  // each sequence unique by LID.
1529  Array<LO> exportLIDs1Copy(exportLIDs1.begin(), exportLIDs1.end());
1530  Array<int> exportPIDs1Copy(exportLIDs1.begin(), exportLIDs1.end());
1531  sort2(exportLIDs1Copy.begin(), exportLIDs1Copy.end(),
1532  exportPIDs1Copy.begin());
1533  typename ArrayView<LO>::iterator exportLIDs1_end = exportLIDs1Copy.end();
1534  typename ArrayView<LO>::iterator exportPIDs1_end = exportPIDs1Copy.end();
1535  merge2(exportLIDs1_end, exportPIDs1_end,
1536  exportLIDs1Copy.begin(), exportLIDs1_end,
1537  exportPIDs1Copy.begin(), exportPIDs1_end,
1538  project1st<LO, LO>());
1539 
1540  Array<LO> exportLIDs2Copy(exportLIDs2.begin(), exportLIDs2.end());
1541  Array<int> exportPIDs2Copy(exportLIDs2.begin(), exportLIDs2.end());
1542  sort2(exportLIDs2Copy.begin(), exportLIDs2Copy.end(),
1543  exportPIDs2Copy.begin());
1544  typename ArrayView<LO>::iterator exportLIDs2_end = exportLIDs2Copy.end();
1545  typename ArrayView<LO>::iterator exportPIDs2_end = exportPIDs2Copy.end();
1546  merge2(exportLIDs2_end, exportPIDs2_end,
1547  exportLIDs2Copy.begin(), exportLIDs2_end,
1548  exportPIDs2Copy.begin(), exportPIDs2_end,
1549  project1st<LO, LO>());
1550 
1551  // Merge export (LID,PID) pairs. In this merge operation, the
1552  // LIDs are the "keys" and the PIDs their "values." We combine
1553  // the "values" (PIDs) in the pairs by replacement, rather than
1554  // by adding them together.
1555  keyValueMerge(exportLIDs1Copy.begin(), exportLIDs1Copy.end(),
1556  exportPIDs1Copy.begin(), exportPIDs1Copy.end(),
1557  exportLIDs2Copy.begin(), exportLIDs2Copy.end(),
1558  exportPIDs2Copy.begin(), exportPIDs2Copy.end(),
1559  std::back_inserter(exportLIDsUnion),
1560  std::back_inserter(exportPIDsUnion),
1562 
1563  // Resort the merged (LID,PID) pairs by PID.
1564  sort2(exportPIDsUnion.begin(), exportPIDsUnion.end(),
1565  exportLIDsUnion.begin());
1566 
1567  // Initialize the Distributor. Using createFromSends instead of
1568  // createFromRecvs avoids the initialization and use of a
1569  // temporary Distributor object.
1570  (void)distributor.createFromSends(exportPIDsUnion().getConst());
1571 #else // NOT TPETRA_IMPORT_SETUNION_USE_CREATE_FROM_SENDS
1572 
1573  // Call the Distributor's createFromRecvs() method to turn the
1574  // remote GIDs and their owning processes into a send-and-receive
1575  // communication plan. remoteGIDsUnion and remotePIDsUnion are
1576  // input; exportGIDsUnion and exportPIDsUnion are output arrays
1577  // which are allocated by createFromRecvs().
1578  distributor.createFromRecvs(remoteGIDsUnion().getConst(),
1579  remotePIDsUnion().getConst(),
1580  exportGIDsUnion, exportPIDsUnion);
1581 
1582  // Find the (source Map) LIDs corresponding to the export GIDs.
1583  const size_type numExportIDsUnion = exportGIDsUnion.size();
1584  exportLIDsUnion.resize(numExportIDsUnion);
1585  for (size_type k = 0; k < numExportIDsUnion; ++k) {
1586  exportLIDsUnion[k] = srcMap->getLocalElement(exportGIDsUnion[k]);
1587  }
1588 #endif // TPETRA_IMPORT_SETUNION_USE_CREATE_FROM_SENDS
1589 
1590  // Create and return the union Import. This uses the "expert" constructor
1591 #ifdef HAVE_TPETRA_MMM_TIMINGS
1592  MM3.disableTimer(label);
1593  label = "Tpetra::Import::setUnion : Construct Import";
1594  TimeMonitor MM4(*TimeMonitor::getNewTimer(label));
1595 #endif
1596  RCP<const import_type> unionImport =
1597  rcp(new import_type(srcMap, unionTgtMap,
1598  as<size_t>(numSameIDsUnion),
1599  permuteToLIDsUnion, permuteFromLIDsUnion,
1600  remoteLIDsUnion, exportLIDsUnion,
1601  exportPIDsUnion, distributor,
1602  this->TransferData_->out_));
1603  return unionImport;
1604 }
1605 
1606 template <class LocalOrdinal, class GlobalOrdinal, class Node>
1607 Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node>>
1609  setUnion() const {
1610  using Teuchos::Array;
1611  using Teuchos::ArrayView;
1612  using Teuchos::as;
1613  using Teuchos::Comm;
1614  using Teuchos::outArg;
1615  using Teuchos::RCP;
1616  using Teuchos::rcp;
1617  using Teuchos::REDUCE_MIN;
1618  using Teuchos::reduceAll;
1619  typedef LocalOrdinal LO;
1620  typedef GlobalOrdinal GO;
1621  Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node>> unionImport;
1622  RCP<const map_type> srcMap = this->getSourceMap();
1623  RCP<const map_type> tgtMap = this->getTargetMap();
1624  RCP<const Comm<int>> comm = srcMap->getComm();
1625 
1626  ArrayView<const GO> srcGIDs = srcMap->getLocalElementList();
1627  ArrayView<const GO> tgtGIDs = tgtMap->getLocalElementList();
1628 
1629  // All elements in srcMap will be in the "new" target map, so...
1630  size_t numSameIDsNew = srcMap->getLocalNumElements();
1631  size_t numRemoteIDsNew = this->getNumRemoteIDs();
1632  Array<LO> permuteToLIDsNew, permuteFromLIDsNew; // empty on purpose
1633 
1634  // Grab some old data
1635  ArrayView<const LO> remoteLIDsOld = this->getRemoteLIDs();
1636  ArrayView<const LO> exportLIDsOld = this->getExportLIDs();
1637 
1638  // Build up the new map (same part)
1639  Array<GO> GIDs(numSameIDsNew + numRemoteIDsNew);
1640  for (size_t i = 0; i < numSameIDsNew; i++)
1641  GIDs[i] = srcGIDs[i];
1642 
1643  // Build up the new map (remote part) and remotes list
1644  Array<LO> remoteLIDsNew(numRemoteIDsNew);
1645  for (size_t i = 0; i < numRemoteIDsNew; i++) {
1646  GIDs[numSameIDsNew + i] = tgtGIDs[remoteLIDsOld[i]];
1647  remoteLIDsNew[i] = numSameIDsNew + i;
1648  }
1649 
1650  // Build the new target map
1651  GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
1652  RCP<const map_type> targetMapNew =
1653  rcp(new map_type(GO_INVALID, GIDs, tgtMap->getIndexBase(),
1654  tgtMap->getComm()));
1655 
1656  // Exports are trivial (since the sourcemap doesn't change)
1657  Array<int> exportPIDsnew(this->getExportPIDs());
1658  Array<LO> exportLIDsnew(this->getExportLIDs());
1659 
1660  // Copy the Distributor (due to how the Import constructor works)
1661  Distributor D(this->getDistributor());
1662 
1663  // Build the importer using the "expert" constructor
1664  unionImport = rcp(new Import<LocalOrdinal, GlobalOrdinal, Node>(srcMap,
1665  targetMapNew,
1666  numSameIDsNew,
1667  permuteToLIDsNew,
1668  permuteFromLIDsNew,
1669  remoteLIDsNew,
1670  exportLIDsnew,
1671  exportPIDsnew, D));
1672 
1673  return unionImport;
1674 }
1675 
1676 template <class LocalOrdinal, class GlobalOrdinal, class Node>
1677 Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node>>
1679  createRemoteOnlyImport(const Teuchos::RCP<const map_type>& remoteTarget) const {
1680  using std::endl;
1681  using Teuchos::outArg;
1682  using Teuchos::rcp;
1683  using Teuchos::REDUCE_MIN;
1684  using Teuchos::reduceAll;
1685  using ::Tpetra::Details::Behavior;
1686  using ::Tpetra::Details::gathervPrint;
1687  using LO = LocalOrdinal;
1688  using GO = GlobalOrdinal;
1689  using import_type = Import<LocalOrdinal, GlobalOrdinal, Node>;
1690 
1691  const char funcPrefix[] = "Tpetra::createRemoteOnlyImport: ";
1692  int lclSuccess = 1;
1693  int gblSuccess = 1;
1694  const bool debug = Behavior::debug();
1695 
1696  const size_t NumRemotes = this->getNumRemoteIDs();
1697  std::unique_ptr<std::string> procPrefix;
1698  Teuchos::RCP<const Teuchos::Comm<int>> comm;
1699  if (debug) {
1700  comm = remoteTarget.is_null() ? Teuchos::null : remoteTarget->getComm();
1701  std::ostringstream os;
1702  os << "Proc ";
1703  if (comm.is_null()) {
1704  os << "?";
1705  } else {
1706  os << comm->getRank();
1707  }
1708  os << ": ";
1709  procPrefix = std::unique_ptr<std::string>(new std::string(os.str()));
1710  }
1711 
1712  if (debug) {
1713  std::ostringstream lclErr;
1714  if (remoteTarget.is_null()) {
1715  lclSuccess = -1;
1716  } else if (NumRemotes != remoteTarget->getLocalNumElements()) {
1717  lclSuccess = 0;
1718  lclErr << *procPrefix << "getNumRemoteIDs() = " << NumRemotes
1719  << " != remoteTarget->getLocalNumElements() = "
1720  << remoteTarget->getLocalNumElements() << "." << endl;
1721  }
1722 
1723  if (comm.is_null()) {
1724  lclSuccess = gblSuccess;
1725  } else {
1726  reduceAll(*comm, REDUCE_MIN, lclSuccess, outArg(gblSuccess));
1727  }
1728  TEUCHOS_TEST_FOR_EXCEPTION(gblSuccess == -1, std::invalid_argument, funcPrefix << "Input target Map is null on at least one process.");
1729 
1730  if (gblSuccess != 1) {
1731  if (comm.is_null()) {
1732  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, lclErr.str());
1733  } else {
1734  std::ostringstream gblErr;
1735  gblErr << funcPrefix << endl;
1736  gathervPrint(gblErr, lclErr.str(), *comm);
1737  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, gblErr.str());
1738  }
1739  }
1740  }
1741 
1742  // Compute the new Remote LIDs
1743  Teuchos::ArrayView<const LO> oldRemoteLIDs = this->getRemoteLIDs();
1744  Teuchos::Array<LO> newRemoteLIDs(NumRemotes);
1745  const map_type& tgtMap = *(this->getTargetMap());
1746  size_t badCount = 0;
1747 
1748  std::unique_ptr<std::vector<size_t>> badIndices;
1749  if (debug) {
1750  badIndices = std::unique_ptr<std::vector<size_t>>(new std::vector<size_t>);
1751  }
1752 
1753  for (size_t i = 0; i < NumRemotes; ++i) {
1754  const LO oldLclInd = oldRemoteLIDs[i];
1755  if (oldLclInd == Teuchos::OrdinalTraits<LO>::invalid()) {
1756  ++badCount;
1757  if (debug) {
1758  badIndices->push_back(i);
1759  }
1760  continue;
1761  }
1762  const GO gblInd = tgtMap.getGlobalElement(oldLclInd);
1763  if (gblInd == Teuchos::OrdinalTraits<GO>::invalid()) {
1764  ++badCount;
1765  if (debug) {
1766  badIndices->push_back(i);
1767  }
1768  continue;
1769  }
1770  const LO newLclInd = remoteTarget->getLocalElement(gblInd);
1771  if (newLclInd == Teuchos::OrdinalTraits<LO>::invalid()) {
1772  ++badCount;
1773  if (debug) {
1774  badIndices->push_back(i);
1775  }
1776  continue;
1777  }
1778  newRemoteLIDs[i] = newLclInd;
1779  // Now we make sure these guys are in sorted order (AztecOO-ML ordering)
1780  if (i > 0 && newRemoteLIDs[i] < newRemoteLIDs[i - 1]) {
1781  ++badCount;
1782  if (debug) {
1783  badIndices->push_back(i);
1784  }
1785  }
1786  }
1787 
1788  if (badCount != 0) {
1789  lclSuccess = 0;
1790  }
1791 
1792  if (debug) {
1793  if (comm.is_null()) {
1794  lclSuccess = gblSuccess;
1795  } else {
1796  reduceAll(*comm, REDUCE_MIN, lclSuccess, outArg(gblSuccess));
1797  }
1798  std::ostringstream lclErr;
1799  if (lclSuccess != 1) {
1800  lclErr << *procPrefix << "Count of bad indices: " << badCount
1801  << ", bad indices: [";
1802  // TODO (mfh 04 Sep 2019) Limit the maximum number of bad
1803  // indices to print.
1804  for (size_t k = 0; k < badCount; ++k) {
1805  const size_t badIndex = (*badIndices)[k];
1806  lclErr << "(" << badIndex << ","
1807  << oldRemoteLIDs[badIndex] << ")";
1808  if (k + size_t(1) < badCount) {
1809  lclErr << ", ";
1810  }
1811  }
1812  lclErr << "]" << endl;
1813  }
1814 
1815  if (gblSuccess != 1) {
1816  std::ostringstream gblErr;
1817  gblErr << funcPrefix << "this->getRemoteLIDs() has \"bad\" "
1818  "indices on one or more processes. \"Bad\" means that the "
1819  "indices are invalid, they don't exist in the target Map, "
1820  "they don't exist in remoteTarget, or they are not in "
1821  "sorted order. In what follows, I will show the \"bad\" "
1822  "indices as (k, LID) pairs, where k is the zero-based "
1823  "index of the LID in this->getRemoteLIDs()."
1824  << endl;
1825  if (comm.is_null()) {
1826  gblErr << lclErr.str();
1827  } else {
1828  gathervPrint(gblErr, lclErr.str(), *comm);
1829  }
1830  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, gblErr.str());
1831  }
1832  } else { // not debug
1833  TEUCHOS_TEST_FOR_EXCEPTION(lclSuccess == 0, std::runtime_error, funcPrefix << "this->getRemoteLIDs() has " << badCount << "ind" << (badCount == 1 ? "ex" : "ices") << " \"bad\" indices on this process." << endl);
1834  }
1835 
1836  // Copy ExportPIDs and such
1837  // NOTE: Be careful: The "Expert" Import constructor we use does a "swap"
1838  // for most of the LID/PID lists and the Distributor, meaning it
1839  // ruins the existing object if we pass things in directly. Hence
1840  // we copy them first.
1841  Teuchos::Array<int> newExportPIDs(this->getExportPIDs());
1842  Teuchos::Array<LO> newExportLIDs(this->getExportLIDs());
1843  Teuchos::Array<LO> dummy;
1844  Distributor newDistor(this->getDistributor());
1845 
1846  return rcp(new import_type(this->getSourceMap(), remoteTarget,
1847  static_cast<size_t>(0), dummy, dummy,
1848  newRemoteLIDs, newExportLIDs,
1849  newExportPIDs, newDistor));
1850 }
1851 
1852 } // namespace Tpetra
1853 
1854 #define TPETRA_IMPORT_CLASS_INSTANT(LO, GO, NODE) \
1855  template class Import<LO, GO, NODE>;
1856 
1857 // Explicit instantiation macro.
1858 // Only invoke this when in the Tpetra namespace.
1859 // Most users do not need to use this.
1860 //
1861 // LO: The local ordinal type.
1862 // GO: The global ordinal type.
1863 // NODE: The Kokkos Node type.
1864 #define TPETRA_IMPORT_INSTANT(LO, GO, NODE) \
1865  TPETRA_IMPORT_CLASS_INSTANT(LO, GO, NODE)
1866 
1867 #endif // TPETRA_IMPORT_DEF_HPP
bool congruent(const Teuchos::Comm< int > &comm1, const Teuchos::Comm< int > &comm2)
Whether the two communicators are congruent.
Definition: Tpetra_Util.cpp:33
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
void merge2(IT1 &indResultOut, IT2 &valResultOut, IT1 indBeg, IT1 indEnd, IT2 valBeg, IT2)
Merge values in place, additively, with the same index.
Teuchos::FancyOStream & verboseOutputStream() const
Valid (nonnull) output stream for verbose output.
Teuchos::RCP< const Import< LocalOrdinal, GlobalOrdinal, Node > > createRemoteOnlyImport(const Teuchos::RCP< const map_type > &remoteTarget) const
Returns an importer that contains only the remote entries of this.
Declaration of Tpetra::Details::Profiling, a scope guard for Kokkos Profiling.
LookupStatus
Return status of Map remote index lookup (getRemoteIndexList()).
Declaration of a function that prints strings from each process.
void sort3(const IT1 &first1, const IT1 &last1, const IT2 &first2, const IT3 &first3, const bool stableSort=false)
Sort the first array, and apply the same permutation to the second and third arrays.
Teuchos::ArrayView< const int > getExportPIDs() const
List of processes to which entries will be sent.
Teuchos::ArrayView< const LocalOrdinal > getPermuteToLIDs() const
List of local IDs in the target Map that are permuted.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Describe this object in a human-readable way to the given output stream.
Teuchos::RCP< const map_type > getTargetMap() const
The target Map used to construct this Export or Import.
size_t getNumSameIDs() const
Number of initial identical IDs.
void sort(View &view, const size_t &size)
Convenience wrapper for std::sort for host-accessible views.
size_t getNumRemoteIDs() const
Number of entries not on the calling process.
size_t createFromSends(const Teuchos::ArrayView< const int > &exportProcIDs)
Set up Distributor using list of process ranks to which this process will send.
Teuchos::RCP< const Import< LocalOrdinal, GlobalOrdinal, Node > > setUnion() const
Return the union of this Import this-&gt;getSourceMap()
size_t global_size_t
Global size_t object.
Import(const Teuchos::RCP< const map_type > &source, const Teuchos::RCP< const map_type > &target)
Construct an Import from the source and target Maps.
Teuchos::ArrayView< const LocalOrdinal > getRemoteLIDs() const
List of entries in the target Map to receive from other processes.
global_ordinal_type getGlobalElement(local_ordinal_type localIndex) const
The global index corresponding to the given local index.
void keyValueMerge(KeyInputIterType keyBeg1, KeyInputIterType keyEnd1, ValueInputIterType valBeg1, ValueInputIterType valEnd1, KeyInputIterType keyBeg2, KeyInputIterType keyEnd2, ValueInputIterType valBeg2, ValueInputIterType valEnd2, KeyOutputIterType keyOut, ValueOutputIterType valOut, BinaryFunction f)
Merge two sorted (by keys) sequences of unique (key,value) pairs by combining pairs with equal keys...
void createFromSendsAndRecvs(const Teuchos::ArrayView< const int > &exportProcIDs, const Teuchos::ArrayView< const int > &remoteProcIDs)
Set up Distributor using list of process ranks to which to send, and list of process ranks from which...
auto view_alloc_no_init(const std::string &label) -> decltype(Kokkos::view_alloc(label, Kokkos::WithoutInitializing))
Use in place of the string label as the first argument of Kokkos::View&#39;s constructor, in case you want to allocate without initializing.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
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.
#define TPETRA_ABUSE_WARNING(throw_exception_test, Exception, msg)
Handle an abuse warning, according to HAVE_TPETRA_THROW_ABUSE_WARNINGS and HAVE_TPETRA_PRINT_ABUSE_WA...
Implementation detail of Import and Export.
Sets up and executes a communication plan for a Tpetra DistObject.
Teuchos::RCP< ImportExportData< LocalOrdinal, GlobalOrdinal, Node > > TransferData_
All the data needed for executing the Export communication plan.
Teuchos::RCP< const map_type > getSourceMap() const
The source Map used to construct this Export or Import.
void findUnionTargetGIDs(Teuchos::Array< GlobalOrdinal > &unionTgtGIDs, Teuchos::Array< std::pair< int, GlobalOrdinal >> &remotePGIDs, typename Teuchos::Array< GlobalOrdinal >::size_type &numSameGIDs, typename Teuchos::Array< GlobalOrdinal >::size_type &numPermuteGIDs, typename Teuchos::Array< GlobalOrdinal >::size_type &numRemoteGIDs, const Teuchos::ArrayView< const GlobalOrdinal > &sameGIDs1, const Teuchos::ArrayView< const GlobalOrdinal > &sameGIDs2, Teuchos::Array< GlobalOrdinal > &permuteGIDs1, Teuchos::Array< GlobalOrdinal > &permuteGIDs2, Teuchos::Array< GlobalOrdinal > &remoteGIDs1, Teuchos::Array< GlobalOrdinal > &remoteGIDs2, Teuchos::Array< int > &remotePIDs1, Teuchos::Array< int > &remotePIDs2) const
Find the union of the target IDs from two Import objects.
void start()
Start the deep_copy counter.
Teuchos::ArrayView< const LocalOrdinal > getExportLIDs() const
List of entries in the source Map that will be sent to other processes.
A parallel distribution of indices over processes.
Internal functions and macros designed for use with Tpetra::Import and Tpetra::Export objects...
Stand-alone utility functions and macros.
void makeDualViewFromOwningHostView(Kokkos::DualView< ElementType *, DeviceType > &dv, const typename Kokkos::DualView< ElementType *, DeviceType >::t_host &hostView)
Initialize dv such that its host View is hostView.
bool verbose() const
Whether to print verbose debugging output.
virtual void print(std::ostream &os) const
Print the Import&#39;s data to the given output stream.
void createFromRecvs(const Teuchos::ArrayView< const Ordinal > &remoteIDs, const Teuchos::ArrayView< const int > &remoteProcIDs, Teuchos::Array< Ordinal > &exportIDs, Teuchos::Array< int > &exportProcIDs)
Set up Distributor using list of process ranks from which to receive.
Binary function that returns its first argument.
Declaration of Tpetra::Details::Behavior, a class that describes Tpetra&#39;s behavior.