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