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 <Teuchos_as.hpp>
53 
54 
55 namespace {
56  // Default value of Import's "Debug" parameter.
57  const bool tpetraImportDebugDefault = false;
58 } // namespace (anonymous)
59 
60 namespace Tpetra {
61 
62  template <class LocalOrdinal, class GlobalOrdinal, class Node>
63  void
64  Import<LocalOrdinal,GlobalOrdinal,Node>::
65  setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
66  {
67  bool debug = tpetraImportDebugDefault;
68  if (! plist.is_null ()) {
69  try {
70  debug = plist->get<bool> ("Debug");
71  } catch (Teuchos::Exceptions::InvalidParameter&) {}
72  }
73  debug_ = debug;
74  ImportData_->distributor_.setParameterList (plist);
75  }
76 
77  template <class LocalOrdinal, class GlobalOrdinal, class Node>
78  void
80  init (const Teuchos::RCP<const map_type>& source,
81  const Teuchos::RCP<const map_type>& target,
82  bool useRemotePIDs,
83  Teuchos::Array<int> & remotePIDs,
84  const Teuchos::RCP<Teuchos::ParameterList>& plist)
85  {
86  using Teuchos::null;
87  using Teuchos::Ptr;
88  using Teuchos::rcp;
89  using std::endl;
91 
92  // Read "Debug" parameter from the input ParameterList.
93  bool debug = tpetraImportDebugDefault;
94  if (! plist.is_null ()) {
95  try {
96  debug = plist->get<bool> ("Debug");
97  } catch (Teuchos::Exceptions::InvalidParameter&) {}
98  }
99  debug_ = debug;
100 
101  if (! out_.is_null ()) {
102  out_->pushTab ();
103  }
104  if (debug_) {
105  std::ostringstream os;
106  const int myRank = source->getComm ()->getRank ();
107  os << myRank << ": Import ctor" << endl;
108  *out_ << os.str ();
109  }
110  ImportData_ = rcp (new data_type (source, target, out_, plist));
111 
112  Array<GlobalOrdinal> remoteGIDs;
113  setupSamePermuteRemote (remoteGIDs);
114  if (debug_) {
115  std::ostringstream os;
116  const int myRank = source->getComm ()->getRank ();
117  os << myRank << ": Import ctor: "
118  << "setupSamePermuteRemote done" << endl;
119  *out_ << os.str ();
120  }
121  if (source->isDistributed ()) {
122  setupExport (remoteGIDs,useRemotePIDs,remotePIDs);
123  }
124  if (debug_) {
125  std::ostringstream os;
126  const int myRank = source->getComm ()->getRank ();
127  os << myRank << ": Import ctor: done" << endl;
128  *out_ << os.str ();
129  }
130  if (! out_.is_null ()) {
131  out_->popTab ();
132  }
133  }
134 
135  template <class LocalOrdinal, class GlobalOrdinal, class Node>
137  Import (const RCP<const map_type >& source,
138  const RCP<const map_type >& target) :
139  out_ (Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr))),
140  debug_ (tpetraImportDebugDefault)
141  {
142  Teuchos::Array<int> dummy;
143  init (source, target, false, dummy, Teuchos::null);
144  }
145 
146  template <class LocalOrdinal, class GlobalOrdinal, class Node>
148  Import (const RCP<const map_type >& source,
149  const RCP<const map_type >& target,
150  const RCP<Teuchos::FancyOStream>& out) :
151  out_ (out),
152  debug_ (tpetraImportDebugDefault)
153  {
154  Teuchos::Array<int> dummy;
155  init (source, target, false, dummy, Teuchos::null);
156  }
157 
158  template <class LocalOrdinal, class GlobalOrdinal, class Node>
160  Import (const Teuchos::RCP<const map_type >& source,
161  const Teuchos::RCP<const map_type >& target,
162  const Teuchos::RCP<Teuchos::ParameterList>& plist) :
163  out_ (Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr))),
164  debug_ (tpetraImportDebugDefault)
165  {
166  Teuchos::Array<int> dummy;
167  init (source, target, false, dummy, plist);
168  }
169 
170  template <class LocalOrdinal, class GlobalOrdinal, class Node>
172  Import (const Teuchos::RCP<const map_type >& source,
173  const Teuchos::RCP<const map_type >& target,
174  const RCP<Teuchos::FancyOStream>& out,
175  const Teuchos::RCP<Teuchos::ParameterList>& plist) :
176  out_ (out),
177  debug_ (tpetraImportDebugDefault)
178  {
179  Teuchos::Array<int> dummy;
180  init (source, target, false, dummy, plist);
181  }
182 
183  template <class LocalOrdinal, class GlobalOrdinal, class Node>
185  Import (const Teuchos::RCP<const map_type >& source,
186  const Teuchos::RCP<const map_type >& target,
187  Teuchos::Array<int> & remotePIDs) :
188  debug_ (tpetraImportDebugDefault)
189  {
190  init (source, target, true, remotePIDs, Teuchos::null);
191  }
192 
193 
194  template <class LocalOrdinal, class GlobalOrdinal, class Node>
197  : ImportData_ (rhs.ImportData_)
198  , out_ (rhs.out_)
199  , debug_ (rhs.debug_)
200  {}
201 
202  template <class LocalOrdinal, class GlobalOrdinal, class Node>
205  : out_ (exporter.out_)
206  , debug_ (exporter.debug_)
207  {
208  if (! exporter.ExportData_.is_null ()) {
209  ImportData_ = exporter.ExportData_->reverseClone ();
210  }
211  }
212 
213  template <class LocalOrdinal, class GlobalOrdinal, class Node>
215  Import (const Teuchos::RCP<const map_type>& source,
216  const Teuchos::RCP<const map_type>& target,
217  const size_t numSameIDs,
218  Teuchos::Array<LocalOrdinal>& permuteToLIDs,
219  Teuchos::Array<LocalOrdinal>& permuteFromLIDs,
220  Teuchos::Array<LocalOrdinal>& remoteLIDs,
221  Teuchos::Array<LocalOrdinal>& exportLIDs,
222  Teuchos::Array<int>& exportPIDs,
223  Distributor& distributor,
224  const Teuchos::RCP<Teuchos::FancyOStream>& out,
225  const Teuchos::RCP<Teuchos::ParameterList>& plist) :
226  out_ (out.is_null () ? Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr)) : out),
227  debug_ (tpetraImportDebugDefault)
228  {
229  using Teuchos::null;
230  using Teuchos::Ptr;
231  using Teuchos::rcp;
232  using std::cerr;
233  using std::endl;
235 
236  // Read "Debug" parameter from the input ParameterList.
237  bool debug = tpetraImportDebugDefault;
238  if (! plist.is_null ()) {
239  try {
240  debug = plist->get<bool> ("Debug");
241  } catch (Teuchos::Exceptions::InvalidParameter&) {}
242  }
243  debug_ = debug;
244 
245  if (! out_.is_null ()) {
246  out_->pushTab ();
247  }
248  if (debug_ && ! out_.is_null ()) {
249  std::ostringstream os;
250  const int myRank = source->getComm ()->getRank ();
251  os << myRank << ": Import expert ctor" << endl;
252  *out_ << os.str ();
253  }
254  ImportData_ = rcp (new data_type (source, target, out_, plist));
255 
256  ImportData_->numSameIDs_ = numSameIDs;
257  ImportData_->permuteToLIDs_.swap (permuteToLIDs);
258  ImportData_->permuteFromLIDs_.swap (permuteFromLIDs);
259  ImportData_->remoteLIDs_.swap (remoteLIDs);
260  ImportData_->distributor_.swap (distributor);
261  ImportData_->exportLIDs_.swap (exportLIDs);
262  ImportData_->exportPIDs_.swap (exportPIDs);
263  }
264 
265  template <class LocalOrdinal, class GlobalOrdinal, class Node>
267  {}
268 
269  template <class LocalOrdinal, class GlobalOrdinal, class Node>
271  return ImportData_->numSameIDs_;
272  }
273 
274  template <class LocalOrdinal, class GlobalOrdinal, class Node>
276  return ImportData_->permuteFromLIDs_.size();
277  }
278 
279  template <class LocalOrdinal, class GlobalOrdinal, class Node>
280  ArrayView<const LocalOrdinal>
282  return ImportData_->permuteFromLIDs_();
283  }
284 
285  template <class LocalOrdinal, class GlobalOrdinal, class Node>
286  ArrayView<const LocalOrdinal>
288  return ImportData_->permuteToLIDs_();
289  }
290 
291  template <class LocalOrdinal, class GlobalOrdinal, class Node>
293  return ImportData_->remoteLIDs_.size();
294  }
295 
296  template <class LocalOrdinal, class GlobalOrdinal, class Node>
297  ArrayView<const LocalOrdinal>
299  return ImportData_->remoteLIDs_();
300  }
301 
302  template <class LocalOrdinal, class GlobalOrdinal, class Node>
304  return ImportData_->exportLIDs_.size();
305  }
306 
307  template <class LocalOrdinal, class GlobalOrdinal, class Node>
308  ArrayView<const LocalOrdinal>
310  return ImportData_->exportLIDs_();
311  }
312 
313  template <class LocalOrdinal, class GlobalOrdinal, class Node>
314  ArrayView<const int>
316  return ImportData_->exportPIDs_();
317  }
318 
319  template <class LocalOrdinal, class GlobalOrdinal, class Node>
320  Teuchos::RCP<const typename Import<LocalOrdinal,GlobalOrdinal,Node>::map_type>
322  return ImportData_->source_;
323  }
324 
325  template <class LocalOrdinal, class GlobalOrdinal, class Node>
326  Teuchos::RCP<const typename Import<LocalOrdinal,GlobalOrdinal,Node>::map_type>
328  return ImportData_->target_;
329  }
330 
331  template <class LocalOrdinal, class GlobalOrdinal, class Node>
332  Distributor &
334  return ImportData_->distributor_;
335  }
336 
337  template <class LocalOrdinal, class GlobalOrdinal, class Node>
341  if (&rhs != this) {
342  ImportData_ = rhs.ImportData_;
343  }
344  return *this;
345  }
346 
347  template <class LocalOrdinal, class GlobalOrdinal, class Node>
349  print (std::ostream& os) const
350  {
351  using Teuchos::Comm;
352  using Teuchos::getFancyOStream;
353  using Teuchos::RCP;
354  using Teuchos::rcpFromRef;
355  using Teuchos::toString;
356  using std::endl;
357 
358  RCP<const Comm<int> > comm = getSourceMap ()->getComm ();
359  const int myImageID = comm->getRank ();
360  const int numImages = comm->getSize ();
361  for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
362  if (myImageID == imageCtr) {
363  os << endl;
364  if (myImageID == 0) { // I'm the root node (only output this info once)
365  os << "Import Data Members:" << endl;
366  }
367  os << "Image ID : " << myImageID << endl;
368 
369  os << "permuteFromLIDs: " << toString (getPermuteFromLIDs ()) << endl;
370  os << "permuteToLIDs : " << toString (getPermuteToLIDs ()) << endl;
371  os << "remoteLIDs : " << toString (getRemoteLIDs ()) << endl;
372  os << "exportLIDs : " << toString (getExportLIDs ()) << endl;
373  os << "exportPIDs : " << toString (getExportPIDs ()) << endl;
374 
375  os << "numSameIDs : " << getNumSameIDs () << endl;
376  os << "numPermuteIDs : " << getNumPermuteIDs () << endl;
377  os << "numRemoteIDs : " << getNumRemoteIDs () << endl;
378  os << "numExportIDs : " << getNumExportIDs () << endl;
379  }
380  // A few global barriers give output a chance to complete.
381  comm->barrier();
382  comm->barrier();
383  comm->barrier();
384  }
385 
386  const bool printMaps = false;
387  if (printMaps) {
388  if (myImageID == 0) {
389  os << endl << endl << "Source Map:" << endl << std::flush;
390  }
391  comm->barrier();
392  os << *getSourceMap();
393  comm->barrier();
394 
395  if (myImageID == 0) {
396  os << endl << endl << "Target Map:" << endl << std::flush;
397  }
398  comm->barrier();
399  os << *getTargetMap();
400  comm->barrier();
401  }
402 
403  // It's also helpful for debugging to print the Distributor
404  // object. Epetra_Import::Print() does this, so we can do a
405  // side-by-side comparison.
406  if (myImageID == 0) {
407  os << endl << endl << "Distributor:" << endl << std::flush;
408  }
409  comm->barrier();
410  getDistributor().describe (*(getFancyOStream (rcpFromRef (os))),
411  Teuchos::VERB_EXTREME);
412  }
413 
414 
415  template <class LocalOrdinal, class GlobalOrdinal, class Node>
416  void
418  setupSamePermuteRemote (Teuchos::Array<GlobalOrdinal>& remoteGIDs)
419  {
420  using Teuchos::arcp;
421  using Teuchos::Array;
422  using Teuchos::ArrayRCP;
423  using Teuchos::ArrayView;
424  using Teuchos::as;
425  using Teuchos::null;
426  typedef LocalOrdinal LO;
427  typedef GlobalOrdinal GO;
428  typedef typename ArrayView<const GO>::size_type size_type;
429  const map_type& source = * (getSourceMap ());
430  const map_type& target = * (getTargetMap ());
431  ArrayView<const GO> sourceGIDs = source.getNodeElementList ();
432  ArrayView<const GO> targetGIDs = target.getNodeElementList ();
433 
434 #ifdef HAVE_TPETRA_DEBUG
435  ArrayView<const GO> rawSrcGids = sourceGIDs;
436  ArrayView<const GO> rawTgtGids = targetGIDs;
437 #else
438  const GO* const rawSrcGids = sourceGIDs.getRawPtr ();
439  const GO* const rawTgtGids = targetGIDs.getRawPtr ();
440 #endif // HAVE_TPETRA_DEBUG
441  const size_type numSrcGids = sourceGIDs.size ();
442  const size_type numTgtGids = targetGIDs.size ();
443  const size_type numGids = std::min (numSrcGids, numTgtGids);
444 
445  // Compute numSameIDs_: the number of initial GIDs that are the
446  // same (and occur in the same order) in both Maps. The point of
447  // numSameIDs_ is for the common case of an Import where all the
448  // overlapping GIDs are at the end of the target Map, but
449  // otherwise the source and target Maps are the same. This allows
450  // a fast contiguous copy for the initial "same IDs."
451  size_type numSameGids = 0;
452  for ( ; numSameGids < numGids && rawSrcGids[numSameGids] == rawTgtGids[numSameGids]; ++numSameGids)
453  {} // third clause of 'for' does everything
454  ImportData_->numSameIDs_ = numSameGids;
455 
456  // Compute permuteToLIDs_, permuteFromLIDs_, remoteGIDs, and
457  // remoteLIDs_. The first two arrays are IDs to be permuted, and
458  // the latter two arrays are IDs to be received ("imported"),
459  // called "remote" IDs.
460  //
461  // IDs to permute are in both the source and target Maps, which
462  // means we don't have to send or receive them, but we do have to
463  // rearrange (permute) them in general. IDs to receive are in the
464  // target Map, but not the source Map.
465 
466  Array<LO>& permuteToLIDs = ImportData_->permuteToLIDs_;
467  Array<LO>& permuteFromLIDs = ImportData_->permuteFromLIDs_;
468  Array<LO>& remoteLIDs = ImportData_->remoteLIDs_;
469  const LO LINVALID = Teuchos::OrdinalTraits<LO>::invalid ();
470  const LO numTgtLids = as<LO> (numTgtGids);
471  // Iterate over the target Map's LIDs, since we only need to do
472  // GID -> LID lookups for the source Map.
473  for (LO tgtLid = numSameGids; tgtLid < numTgtLids; ++tgtLid) {
474  const GO curTargetGid = rawTgtGids[tgtLid];
475  // getLocalElement() returns LINVALID if the GID isn't in the source Map.
476  // This saves us a lookup (which isNodeGlobalElement() would do).
477  const LO srcLid = source.getLocalElement (curTargetGid);
478  if (srcLid != LINVALID) { // if source.isNodeGlobalElement (curTargetGid)
479  permuteToLIDs.push_back (tgtLid);
480  permuteFromLIDs.push_back (srcLid);
481  } else {
482  remoteGIDs.push_back (curTargetGid);
483  remoteLIDs.push_back (tgtLid);
484  }
485  }
486 
488  getNumRemoteIDs() > 0 && ! source.isDistributed(),
489  std::runtime_error,
490  "::setupSamePermuteRemote(): Target has remote LIDs but Source is not "
491  "distributed globally." << std::endl
492  << "Importing to a submap of the target map.");
493  }
494 
495 
496  template <class LocalOrdinal, class GlobalOrdinal, class Node>
498  setupExport (Teuchos::Array<GlobalOrdinal>& remoteGIDs,
499  bool useRemotePIDs,
500  Teuchos::Array<int>& userRemotePIDs)
501  {
502  using Teuchos::arcp;
503  using Teuchos::Array;
504  using Teuchos::ArrayRCP;
505  using Teuchos::ArrayView;
506  using Teuchos::null;
507  using std::endl;
508  typedef LocalOrdinal LO;
509  typedef GlobalOrdinal GO;
510  typedef typename Array<int>::difference_type size_type;
511 
512  TEUCHOS_TEST_FOR_EXCEPTION(
513  getSourceMap ().is_null (), std::logic_error, "Tpetra::Import::"
514  "setupExport: Source Map is null. Please report this bug to the Tpetra "
515  "developers.");
516  const map_type& source = * (getSourceMap ());
517 
518  Teuchos::OSTab tab (out_);
519 
520  // if (debug_ && ! out_.is_null ()) {
521  // std::ostringstream os;
522  // const int myRank = source.getComm ()->getRank ();
523  // os << myRank << ": Import::setupExport:" << endl;
524  // }
525 
526  // Sanity checks
527  TEUCHOS_TEST_FOR_EXCEPTION(
528  ! useRemotePIDs && (userRemotePIDs.size() > 0), std::invalid_argument,
529  "Tpetra::Import::setupExport: remotePIDs are non-empty but their use has "
530  "not been requested.");
531  TEUCHOS_TEST_FOR_EXCEPTION(
532  userRemotePIDs.size () > 0 && remoteGIDs.size () != userRemotePIDs.size (),
533  std::invalid_argument, "Tpetra::Import::setupExport: remotePIDs must "
534  "either be of size zero or match the size of remoteGIDs.");
535 
536  // For each entry remoteGIDs[i], remoteProcIDs[i] will contain
537  // the process ID of the process that owns that GID.
538  ArrayView<GO> remoteGIDsView = remoteGIDs ();
539  ArrayView<int> remoteProcIDsView;
540 
541  // lookup == IDNotPresent means that the source Map wasn't able to
542  // figure out to which processes one or more of the GIDs in the
543  // given list of remote GIDs belong.
544  //
545  // The previous abuse warning said "The target Map has GIDs not
546  // found in the source Map." This statement could be confusing,
547  // because it doesn't refer to ownership by the current process,
548  // but rather to ownership by _any_ process participating in the
549  // Map. (It could not possibly refer to ownership by the current
550  // process, since remoteGIDs is exactly the list of GIDs owned by
551  // the target Map but not owned by the source Map. It was
552  // constructed that way by setupSamePermuteRemote().)
553  //
554  // What this statement means is that the source and target Maps
555  // don't contain the same set of GIDs globally (over all
556  // processes). That is, there is at least one GID owned by some
557  // process in the target Map, which is not owned by _any_ process
558  // in the source Map.
559  Array<int> newRemotePIDs;
560  LookupStatus lookup = AllIDsPresent;
561 
562  if (! useRemotePIDs) {
563  newRemotePIDs.resize (remoteGIDsView.size ());
564  if (debug_ && ! out_.is_null ()) {
565  std::ostringstream os;
566  const int myRank = source.getComm ()->getRank ();
567  os << myRank << ": Import::setupExport: about to call "
568  "getRemoteIndexList on source Map" << endl;
569  *out_ << os.str ();
570  }
571  lookup = source.getRemoteIndexList (remoteGIDsView, newRemotePIDs ());
572  }
573  Array<int>& remoteProcIDs = useRemotePIDs ? userRemotePIDs : newRemotePIDs;
574 
575  TPETRA_ABUSE_WARNING( lookup == IDNotPresent, std::runtime_error,
576  "::setupExport(): the source Map wasn't able to figure out which process "
577  "owns one or more of the GIDs in the list of remote GIDs. This probably "
578  "means that there is at least one GID owned by some process in the target"
579  " Map which is not owned by any process in the source Map. (That is, the"
580  " source and target Maps do not contain the same set of GIDs globally.)");
581 
582  // Ignore remote GIDs that aren't owned by any process in the
583  // source Map. getRemoteIndexList() gives each of these a process
584  // ID of -1.
585  if (lookup == IDNotPresent) {
586  const size_type numInvalidRemote =
587  std::count_if (remoteProcIDs.begin (), remoteProcIDs.end (),
588  std::bind1st (std::equal_to<int> (), -1));
589  // If all of them are invalid, we can delete the whole array.
590  const size_type totalNumRemote = getNumRemoteIDs ();
591  if (numInvalidRemote == totalNumRemote) {
592  // all remotes are invalid; we have no remotes; we can delete the remotes
593  remoteProcIDs.clear ();
594  remoteGIDs.clear (); // This invalidates the view remoteGIDsView
595  ImportData_->remoteLIDs_.clear();
596  }
597  else {
598  // Some remotes are valid; we need to keep the valid ones.
599  // Pack and resize remoteProcIDs, remoteGIDs, and remoteLIDs_.
600  size_type numValidRemote = 0;
601 #ifdef HAVE_TPETRA_DEBUG
602  ArrayView<GlobalOrdinal> remoteGIDsPtr = remoteGIDsView;
603 #else
604  GlobalOrdinal* const remoteGIDsPtr = remoteGIDsView.getRawPtr ();
605 #endif // HAVE_TPETRA_DEBUG
606  for (size_type r = 0; r < totalNumRemote; ++r) {
607  // Pack in all the valid remote PIDs and GIDs.
608  if (remoteProcIDs[r] != -1) {
609  remoteProcIDs[numValidRemote] = remoteProcIDs[r];
610  remoteGIDsPtr[numValidRemote] = remoteGIDsPtr[r];
611  ImportData_->remoteLIDs_[numValidRemote] = ImportData_->remoteLIDs_[r];
612  ++numValidRemote;
613  }
614  }
615  TEUCHOS_TEST_FOR_EXCEPTION(
616  numValidRemote != totalNumRemote - numInvalidRemote, std::logic_error,
617  "Tpetra::Import::setupExport(): After removing invalid remote GIDs and"
618  " packing the valid remote GIDs, numValidRemote = " << numValidRemote
619  << " != totalNumRemote - numInvalidRemote = "
620  << totalNumRemote - numInvalidRemote
621  << ". Please report this bug to the Tpetra developers.");
622 
623  remoteProcIDs.resize (numValidRemote);
624  remoteGIDs.resize (numValidRemote);
625  ImportData_->remoteLIDs_.resize (numValidRemote);
626  }
627  // Revalidate the view after clear or resize.
628  remoteGIDsView = remoteGIDs ();
629  }
630 
631  // Sort remoteProcIDs in ascending order, and apply the resulting
632  // permutation to remoteGIDs and remoteLIDs_. This ensures that
633  // remoteProcIDs[i], remoteGIDs[i], and remoteLIDs_[i] all refer
634  // to the same thing.
635  sort3 (remoteProcIDs.begin (),
636  remoteProcIDs.end (),
637  remoteGIDsView.begin (),
638  ImportData_->remoteLIDs_.begin ());
639 
640  // Call the Distributor's createFromRecvs() method to turn the
641  // remote GIDs and their owning processes into a send-and-receive
642  // communication plan. remoteGIDs and remoteProcIDs_ are input;
643  // exportGIDs and exportProcIDs_ are output arrays which are
644  // allocated by createFromRecvs().
645  Array<GO> exportGIDs;
646  ImportData_->distributor_.createFromRecvs (remoteGIDsView ().getConst (),
647  remoteProcIDs, exportGIDs,
648  ImportData_->exportPIDs_);
649  // if (debug_ && ! out_.is_null ()) {
650  // std::ostringstream os;
651  // const int myRank = source.getComm ()->getRank ();
652  // os << myRank << ": Import::setupExport: Getting LIDs" << endl;
653  // *out_ << os.str ();
654  // }
655 
656  // Find the LIDs corresponding to the (outgoing) GIDs in
657  // exportGIDs. For sparse matrix-vector multiply, this tells the
658  // calling process how to index into the source vector to get the
659  // elements which it needs to send.
660  //
661  // NOTE (mfh 03 Mar 2014) This is now a candidate for a
662  // thread-parallel kernel, but only if using the new thread-safe
663  // Map implementation.
664  const size_type numExportIDs = exportGIDs.size ();
665  if (numExportIDs > 0) {
666  ImportData_->exportLIDs_.resize (numExportIDs);
667  ArrayView<const GO> expGIDs = exportGIDs ();
668  ArrayView<LO> expLIDs = ImportData_->exportLIDs_ ();
669  for (size_type k = 0; k < numExportIDs; ++k) {
670  expLIDs[k] = source.getLocalElement (expGIDs[k]);
671  }
672  }
673 
674  if (debug_ && ! out_.is_null ()) {
675  std::ostringstream os;
676  const int myRank = source.getComm ()->getRank ();
677  os << myRank << ": Import::setupExport: done" << endl;
678  *out_ << os.str ();
679  }
680  }
681 
682 
683  template <class LocalOrdinal, class GlobalOrdinal, class Node>
684  Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> >
687  {
688  typedef Tpetra::global_size_t GST;
689  using Teuchos::Array;
690  using Teuchos::ArrayView;
691  using Teuchos::as;
692  using Teuchos::Comm;
693  using Teuchos::RCP;
694  using Teuchos::rcp;
695  using Teuchos::outArg;
696  using Teuchos::REDUCE_MIN;
697  using Teuchos::reduceAll;
698  typedef LocalOrdinal LO;
699  typedef GlobalOrdinal GO;
700  typedef Import<LO, GO, Node> import_type;
701  typedef typename Array<GO>::size_type size_type;
702 
703 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
704  using Teuchos::toString;
705  using std::cerr;
706  using std::endl;
707 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
708 
709  RCP<const map_type> srcMap = this->getSourceMap ();
710  RCP<const map_type> tgtMap1 = this->getTargetMap ();
711  RCP<const map_type> tgtMap2 = rhs.getTargetMap ();
712  RCP<const Comm<int> > comm = srcMap->getComm ();
713 
714 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
715  const int myRank = comm->getRank ();
716 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
717 
718 #ifdef HAVE_TPETRA_DEBUG
719  TEUCHOS_TEST_FOR_EXCEPTION(
720  ! srcMap->isSameAs (* (rhs.getSourceMap ())), std::invalid_argument,
721  "Tpetra::Import::setUnion: The source Map of the input Import must be the "
722  "same as (in the sense of Map::isSameAs) the source Map of this Import.");
723  TEUCHOS_TEST_FOR_EXCEPTION(
724  ! Details::congruent (* (tgtMap1->getComm ()), * (tgtMap2->getComm ())),
725  std::invalid_argument, "Tpetra::Import::setUnion: "
726  "The target Maps must have congruent communicators.");
727 #endif // HAVE_TPETRA_DEBUG
728 
729  // It's probably worth the one all-reduce to check whether the two
730  // Maps are the same. If so, we can just return a copy of *this.
731  // isSameAs() bypasses the all-reduce if the pointers are equal.
732  if (tgtMap1->isSameAs (*tgtMap2)) {
733  return rcp (new import_type (*this));
734  }
735  // Alas, the two target Maps are not the same. That means we have
736  // to compute their union, and the union Import object.
737 
738  ArrayView<const GO> srcGIDs = srcMap->getNodeElementList ();
739  ArrayView<const GO> tgtGIDs1 = tgtMap1->getNodeElementList ();
740  ArrayView<const GO> tgtGIDs2 = tgtMap2->getNodeElementList ();
741 
742 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
743  comm->barrier ();
744  {
745  std::ostringstream os;
746  os << myRank << ": srcGIDs: " << toString (srcGIDs) << endl;
747  os << myRank << ": tgtGIDs1: " << toString (tgtGIDs1) << endl;
748  os << myRank << ": tgtGIDs2: " << toString (tgtGIDs2) << endl;
749  cerr << os.str ();
750  }
751  comm->barrier ();
752 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
753 
754 
755  // Fill this as we go with the union target Map's GIDs, in the
756  // desired order. We'll need them for the Map constructor.
757  Array<GO> unionTgtGIDs;
758  // Upper bound on the number of union target Map GIDs. This
759  // happens to be strict, but doesn't have to be. Setting some
760  // reasonable upper bound avoids reallocation in loops that do
761  // push_back operations.
762  unionTgtGIDs.reserve (tgtGIDs1.size () + tgtGIDs2.size ());
763 
764 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
765  if (myRank == 0) {
766  cerr << endl;
767  }
768  comm->barrier ();
769  comm->barrier ();
770  comm->barrier ();
771  cerr << myRank << ": Computing \"same\" GIDs" << endl;
772 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
773 
774  // Compute the initial sequence of "same" GIDs in the union
775  // import. The number of "same" GIDs in the union is the maximum
776  // of the lengths of this in the two inputs.
777 
778  const size_type numSameGIDs1 = this->getNumSameIDs ();
779  const size_type numSameGIDs2 = rhs.getNumSameIDs ();
780  ArrayView<const GO> sameGIDs1 = tgtGIDs1 (0, numSameGIDs1);
781  ArrayView<const GO> sameGIDs2 = tgtGIDs2 (0, numSameGIDs2);
782 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
783  {
784  std::ostringstream os;
785  os << myRank << ": same IDs for target Map 1: " << toString (sameGIDs1) << endl;
786  os << myRank << ": same IDs for target Map 2: " << toString (sameGIDs2) << endl;
787  cerr << os.str ();
788  }
789 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
790  // For the input target Map with fewer "same" GIDs, that Map's
791  // permute IDs could include some of the other input target Map's
792  // "same" GIDs. We have to make sure not to include them twice.
793  // To do so, keep a view of them for now, and remove them (via set
794  // intersection) from the permute ID list. Keeping track of which
795  // GID set had the max number of "same" IDs avoids unnecessary set
796  // intersection operations later.
797  ArrayView<const GO> doubleCountedSameGIDs;
798  size_type numSameIDsUnion;
799  bool tgtMap1HadMaxSameGIDs;
800  if (numSameGIDs1 >= numSameGIDs2) {
801  tgtMap1HadMaxSameGIDs = true;
802  numSameIDsUnion = numSameGIDs1;
803  std::copy (sameGIDs1.begin (), sameGIDs1.end (), std::back_inserter (unionTgtGIDs));
804  // There could be GIDs in target Map 2 that are not included in
805  // the "same" IDs, but are included in Import 2's permute IDs.
806  // Keep track of them so we don't double-count them when
807  // building the list of permute IDs.
808  doubleCountedSameGIDs = tgtGIDs1 (numSameGIDs2, numSameGIDs1 - numSameGIDs2);
809  } else {
810  tgtMap1HadMaxSameGIDs = false;
811  numSameIDsUnion = numSameGIDs2;
812  std::copy (sameGIDs2.begin (), sameGIDs2.end (), std::back_inserter (unionTgtGIDs));
813  // There could be GIDs in target Map 1 that are not included in
814  // the "same" IDs, but are included in Import 1's permute IDs.
815  // Keep track of them so we don't double-count them when
816  // building the list of permute IDs.
817  doubleCountedSameGIDs = tgtGIDs2 (numSameGIDs1, numSameGIDs2 - numSameGIDs1);
818  }
819 
820 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
821  {
822  std::ostringstream os;
823  os << myRank << ": union Map's same GIDs: " << toString (unionTgtGIDs ()) << endl;
824  os << myRank << ": doubleCountedSameGIDs: " << toString (doubleCountedSameGIDs) << endl;
825  cerr << os.str ();
826  }
827  if (myRank == 0) {
828  cerr << endl;
829  }
830  comm->barrier ();
831  comm->barrier ();
832  comm->barrier ();
833  cerr << myRank << ": Computing permute IDs" << endl;
834 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
835 
836  // Each input Import knows its permute-from LIDs (in the source
837  // Map) and permute-to LIDs (in its target Map). We will have to
838  // reassign LIDs in the union target Map, but we can use these
839  // permute-to LIDs to construct the union Import's permute-to IDs.
840  Array<LO> permuteFromLIDsUnion;
841  Array<LO> permuteToLIDsUnion;
842  LO curTgtLid = as<LO> (numSameIDsUnion);
843  {
844  // Permute-to LIDs in the two input target Maps.
845  ArrayView<const LO> permuteToLIDs1 = this->getPermuteToLIDs ();
846  ArrayView<const LO> permuteToLIDs2 = rhs.getPermuteToLIDs ();
847  const size_type numPermuteIDs1 = this->getNumPermuteIDs ();
848  const size_type numPermuteIDs2 = rhs.getNumPermuteIDs ();
849 
850 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
851  cerr << myRank << ": Converting permute-to LIDs to GIDs" << endl;
852 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
853 
854  // Convert the permute-to LID lists to GIDs, so that we can
855  // later reassign LIDs for the (output) union target Map.
856  Array<GO> permuteGIDs1 (numPermuteIDs1);
857  for (size_type k = 0; k < numPermuteIDs1; ++k) {
858  permuteGIDs1[k] = tgtMap1->getGlobalElement (permuteToLIDs1[k]);
859  }
860  Array<GO> permuteGIDs2 (numPermuteIDs2);
861  for (size_type k = 0; k < numPermuteIDs2; ++k) {
862  permuteGIDs2[k] = tgtMap2->getGlobalElement (permuteToLIDs2[k]);
863  }
864 
865 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
866  {
867  std::ostringstream os;
868  os << myRank << ": permuteGIDs1: " << toString (permuteGIDs1) << endl;
869  os << myRank << ": permuteGIDs2: " << toString (permuteGIDs2) << endl;
870  cerr << os.str ();
871  }
872  cerr << myRank << ": Sorting and merging permute GID lists" << endl;
873 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
874 
875  // Sort the two permute GID lists, remove the GIDs that we don't
876  // want to double-count, and merge the result into the global
877  // list of GIDs in the target union Map.
878  std::sort (permuteGIDs1.begin (), permuteGIDs1.end ());
879  std::sort (permuteGIDs2.begin (), permuteGIDs2.end ());
880 
881  typename Array<GO>::iterator permuteGIDs1_beg = permuteGIDs1.begin ();
882  typename Array<GO>::iterator permuteGIDs1_end = permuteGIDs1.end ();
883  typename Array<GO>::iterator permuteGIDs2_beg = permuteGIDs2.begin ();
884  typename Array<GO>::iterator permuteGIDs2_end = permuteGIDs2.end ();
885  if (tgtMap1HadMaxSameGIDs) {
886  // This operation allows the last (output) argument to alias the first.
887  permuteGIDs2_end =
888  std::set_difference(permuteGIDs2_beg,
889  permuteGIDs2_end,
890  doubleCountedSameGIDs.begin (),
891  doubleCountedSameGIDs.end (),
892  permuteGIDs2_beg);
893 
894 
895  } else {
896  // This operation allows the last (output) argument to alias the first.
897  permuteGIDs1_end =
898  std::set_difference(permuteGIDs1_beg,
899  permuteGIDs1_end,
900  doubleCountedSameGIDs.begin (),
901  doubleCountedSameGIDs.end (),
902  permuteGIDs1_beg);
903 
904  }
905  std::set_union (permuteGIDs1_beg, permuteGIDs1_end,
906  permuteGIDs2_beg, permuteGIDs2_end,
907  std::back_inserter (unionTgtGIDs));
908 
909 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
910  {
911  std::ostringstream os;
912  if(tgtMap1HadMaxSameGIDs) os << myRank << ": tgtMap1HadMaxSameGIDs == true"<<endl;
913  else os << myRank << ": tgtMap1HadMaxSameGIDs == false"<<endl;
914 
915  os << myRank << ": reduced permuteGIDs1: {";
916  for(typename Array<GO>::iterator k = permuteGIDs1_beg; k != permuteGIDs1_end; k++)
917  os<<*k<<", ";
918  os<<"}"<<endl;
919  os << myRank << ": reduced permuteGIDs2: {";
920  for(typename Array<GO>::iterator k = permuteGIDs2_beg; k != permuteGIDs2_end; k++)
921  os<<*k<<", ";
922  os<<"}"<<endl;
923  cerr << os.str ();
924  }
925 #endif
926  const size_type numPermuteIDsUnion =
927  unionTgtGIDs.size () - numSameIDsUnion;
928  ArrayView<const GO> permuteGIDsUnion =
929  unionTgtGIDs (numSameIDsUnion, numPermuteIDsUnion).getConst ();
930 
931 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
932  {
933  std::ostringstream os;
934  os << myRank << ": permuteGIDsUnion: " << toString (permuteGIDsUnion) << endl;
935  cerr << os.str ();
936  }
937 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
938 
939 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
940  comm->barrier ();
941  cerr << myRank << ": Computing permute-to LIDs" << endl;
942 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
943 
944  // Compute the permute-to LIDs (in the union target Map).
945  permuteToLIDsUnion.resize (numPermuteIDsUnion);
946  for (size_type k = 0; k < numPermuteIDsUnion; ++k) {
947  permuteToLIDsUnion[k] = curTgtLid++;
948  }
949 
950 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
951  {
952  std::ostringstream os;
953  os << myRank << ": permuteToLIDsUnion: " << toString (permuteToLIDsUnion) << endl;
954  cerr << os.str ();
955  }
956 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
957 
958 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
959  comm->barrier ();
960  cerr << myRank << ": Computing permute-from LIDs" << endl;
961 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
962 
963  // Convert the permute GIDs to permute-from LIDs in the source Map.
964  permuteFromLIDsUnion.resize (numPermuteIDsUnion);
965  for (size_type k = 0; k < numPermuteIDsUnion; ++k) {
966  permuteFromLIDsUnion[k] = srcMap->getLocalElement (permuteGIDsUnion[k]);
967  }
968 
969 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
970  {
971  std::ostringstream os;
972  os << myRank << ": permuteFromLIDsUnion: " << toString (permuteFromLIDsUnion) << endl;
973  cerr << os.str ();
974  }
975 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
976  }// end permutes
977 
978 
979 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
980  {
981  std::ostringstream os;
982  os << myRank << ": unionTgtGIDs after permutes: "
983  << toString (unionTgtGIDs ()) << endl;
984  cerr << os.str ();
985  }
986 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
987 
988  // Thus far, we have computed the following in the union Import:
989  // - getNumSameIDs()
990  // - getNumPermuteIDs()
991  // - getPermuteFromLIDs ()
992  // - getPermuteToLIDs ()
993  //
994  // Now it's time to compute the remote IDs. By definition, none
995  // of these IDs are in the source Map (on the calling process), so
996  // they can't possibly overlap with any of the "same" or permute
997  // IDs in either target Map.
998  //
999  // After the first numSameIDsUnion IDs, we get to control the
1000  // order of GIDs in the union target Map. We'll put the permute
1001  // IDs first (which we already did above) and the remote IDs last
1002  // (which we are about to do). We'll sort the remote IDs by
1003  // process rank, so that Distributor doesn't have to pack buffers.
1004  // (That way, doPosts() will always take the "fast path" on all
1005  // processes.)
1006 
1007 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1008  if (myRank == 0) {
1009  cerr << endl;
1010  }
1011  comm->barrier ();
1012  comm->barrier ();
1013  comm->barrier ();
1014  cerr << myRank << ": Computing remote IDs" << endl;
1015 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1016 
1017  Array<GO> remoteGIDsUnion;
1018  Array<int> remotePIDsUnion;
1019  Array<LO> remoteLIDsUnion;
1020  size_type numRemoteIDsUnion = 0;
1021  {
1022  // Distributor::createFromRecvs takes remote IDs and PIDs as
1023  // input, and computes exportIDs and exportPIDs. The easiest
1024  // way to get the remote PIDs is to imitate setupExport by using
1025  // getRemoteIndexList(). We could try to get them out of the
1026  // Distributor via getImagesFrom(), but Distributor reorders
1027  // them in some not entirely transparent way.
1028 
1029  // Grab the remoteLIDs
1030  ArrayView<const LO> remoteLIDs1 = this->getRemoteLIDs();
1031  ArrayView<const LO> remoteLIDs2 = rhs.getRemoteLIDs();
1032 
1033  // Grab the remotePIDs
1034  Array<int> remotePIDs1, remotePIDs2;
1035  Tpetra::Import_Util::getRemotePIDs(*this,remotePIDs1);
1036  Tpetra::Import_Util::getRemotePIDs(rhs,remotePIDs2);
1037 
1038  // Put the (PID,GID) into std:pairs to make for easier sorting
1039  Array<std::pair<int,GO> > remotePGIDs1, remotePGIDs2,remotePGUnion;
1040  remotePGIDs1.resize(remotePIDs1.size());
1041  remotePGIDs2.resize(remotePIDs2.size());
1042 
1043  for(size_type k=0; k < remotePIDs1.size(); k++)
1044  remotePGIDs1[k] = std::pair<int,GO>(remotePIDs1[k],this->getTargetMap()->getGlobalElement(remoteLIDs1[k]));
1045 
1046  for(size_type k=0; k < remotePIDs2.size(); k++)
1047  remotePGIDs2[k] = std::pair<int,GO>(remotePIDs2[k],rhs.getTargetMap()->getGlobalElement(remoteLIDs2[k]));
1048 
1049 
1050  // Sort and merge the (PID,GID) pairs (with the LIDs along for the ride at least for the sort)
1051  std::sort(remotePGIDs1.begin(), remotePGIDs1.end());
1052  std::sort(remotePGIDs2.begin(), remotePGIDs2.end());
1053  std::merge(remotePGIDs1.begin(), remotePGIDs1.end(),
1054  remotePGIDs2.begin(), remotePGIDs2.end(),
1055  std::back_inserter(remotePGUnion));
1056  typename Array<std::pair<int,GO> >::iterator it = std::unique(remotePGUnion.begin(),remotePGUnion.end());
1057  remotePGUnion.resize(std::distance(remotePGUnion.begin(),it));
1058 
1059  // Assign the remote LIDs in order; copy out
1060  numRemoteIDsUnion = remotePGUnion.size();
1061  remoteLIDsUnion.resize(numRemoteIDsUnion);
1062  remotePIDsUnion.resize(numRemoteIDsUnion);
1063  remoteGIDsUnion.resize(numRemoteIDsUnion);
1064 
1065  for (size_type k = 0; k < numRemoteIDsUnion; ++k) {
1066  remoteLIDsUnion[k] = curTgtLid++;
1067  remotePIDsUnion[k] = remotePGUnion[k].first;
1068  remoteGIDsUnion[k] = remotePGUnion[k].second;
1069  }
1070 
1071  // Update the unionTgtGIDs
1072  const size_type oldSize = unionTgtGIDs.size();
1073  unionTgtGIDs.resize(oldSize + numRemoteIDsUnion);
1074  for(size_type k=0; k<numRemoteIDsUnion; k++)
1075  unionTgtGIDs[oldSize+k] = remoteGIDsUnion[k];
1076 
1077 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1078  {
1079  // For debugging purposes only
1080  Array<GO> remoteGIDs1(remotePIDs1.size());
1081  Array<GO> remoteGIDs2(remotePIDs2.size());
1082  for(size_type k=0; k < remotePIDs1.size(); k++)
1083  remoteGIDs1[k] = this->getTargetMap()->getGlobalElement(remoteLIDs1[k]);
1084  for(size_type k=0; k < remotePIDs2.size(); k++)
1085  remoteGIDs2[k] = rhs.getTargetMap()->getGlobalElement(remoteLIDs2[k]);
1086 
1087  std::ostringstream os;
1088  os << myRank << ": remoteGIDs1 : " << toString (remoteGIDs1 ()) << endl;
1089  os << myRank << ": remotePIDs1 : " << toString (remotePIDs1 ()) << endl;
1090  os << myRank << ": remoteGIDs2 : " << toString (remoteGIDs2 ()) << endl;
1091  os << myRank << ": remotePIDs2 : " << toString (remotePIDs2 ()) << endl;
1092  cerr << os.str ();
1093  }
1094 #endif
1095  }//end remotes
1096 
1097 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1098  {
1099  std::ostringstream os;
1100  os << myRank << ": remoteGIDsUnion sorted: " << toString (remoteGIDsUnion ()) << endl;
1101  os << myRank << ": remotePIDsUnion sorted: " << toString (remotePIDsUnion ()) << endl;
1102  os << myRank << ": remoteLIDsUnion sorted: " << toString (remoteLIDsUnion ()) << endl;
1103  cerr << os.str ();
1104  }
1105 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1106 
1107 
1108 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1109  {
1110  std::ostringstream os;
1111  os << myRank << ": unionTgtGIDs after remotes: "
1112  << toString (unionTgtGIDs ()) << endl;
1113  cerr << os.str ();
1114  }
1115 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1116 
1117  // Find the union target Map's index base, which must also be the
1118  // union target Map's global min GID. Thus, by definition, it
1119  // must be the minimum of the two input target Maps' index bases.
1120  // We already know these, so we don't have to do another
1121  // all-reduce to find it.
1122  const GO indexBaseUnion =
1123  std::min (tgtMap1->getIndexBase (), tgtMap2->getIndexBase ());
1124 
1125 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1126  cerr << myRank << "Creating union target Map" << endl;
1127 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1128 
1129  // Create the union target Map.
1130  //
1131  // mfh 01 May 2013, 28 Feb 2014: It might be handy to have a Map
1132  // constructor that takes the global min and max GIDs; that would
1133  // obviate the need for Map to compute them. On the other hand,
1134  // for signed GlobalOrdinal, this version of Map's constructor
1135  // already computes as few all-reduces as possible (not including
1136  // optimizations that might be possible if one were to fold in
1137  // Directory construction). The constructor must do two
1138  // all-reduces:
1139  //
1140  // 1. Get the global number of GIDs (since the first argument is
1141  // INVALID, we're asking Map to compute this for us)
1142  //
1143  // 2. Figure out three things: global min and max GID, and
1144  // whether the Map is distributed or locally replicated.
1145  //
1146  // #2 above happens in one all-reduce (of three integers).
1147  // #Figuring out whether the Map is distributed or locally
1148  // #replicated requires knowing the global number of GIDs.
1149  const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();
1150  RCP<const map_type> unionTgtMap =
1151  rcp (new map_type (INVALID, unionTgtGIDs (), indexBaseUnion,
1152  comm, srcMap->getNode ()));
1153 
1154 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1155  if (myRank == 0) {
1156  cerr << endl;
1157  }
1158  comm->barrier ();
1159  comm->barrier ();
1160  comm->barrier ();
1161  cerr << myRank << ": Computing export IDs and Distributor" << endl;
1162 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1163 
1164  // Thus far, we have computed the following in the union Import:
1165  // - numSameIDs
1166  // - numPermuteIDs and permuteFromLIDs
1167  // - numRemoteIDs, remoteGIDs, remoteLIDs, and remotePIDs
1168  //
1169  // Now it's time to compute the export IDs and initialize the
1170  // Distributor.
1171 
1172  Array<GO> exportGIDsUnion;
1173  Array<LO> exportLIDsUnion;
1174  Array<int> exportPIDsUnion;
1175  Distributor distributor (comm, this->out_);
1176 
1177 #ifdef TPETRA_IMPORT_SETUNION_USE_CREATE_FROM_SENDS
1178  // Compute the export IDs without communication, by merging the
1179  // lists of (export LID, export PID) pairs from the two input
1180  // Import objects. The export LIDs in both input Import objects
1181  // are LIDs in the source Map. Then, use the export PIDs to
1182  // initialize the Distributor via createFromSends.
1183 
1184  // const size_type numExportIDs1 = this->getNumExportIDs ();
1185  ArrayView<const LO> exportLIDs1 = this->getExportLIDs ();
1186  ArrayView<const LO> exportPIDs1 = this->getExportPIDs ();
1187 
1188  // const size_type numExportIDs2 = rhs.getNumExportIDs ();
1189  ArrayView<const LO> exportLIDs2 = rhs.getExportLIDs ();
1190  ArrayView<const LO> exportPIDs2 = rhs.getExportPIDs ();
1191 
1192  // We have to keep the export LIDs in PID-sorted order, then merge
1193  // them. So, first key-value merge (LID,PID) pairs, treating PIDs
1194  // as values, merging values by replacement. Then, sort the
1195  // (LID,PID) pairs again by PID.
1196 
1197  // Sort (LID,PID) pairs by LID for the later merge, and make
1198  // each sequence unique by LID.
1199  Array<LO> exportLIDs1Copy (exportLIDs1.begin (), exportLIDs1.end ());
1200  Array<int> exportPIDs1Copy (exportLIDs1.begin (), exportLIDs1.end ());
1201  sort2 (exportLIDs1Copy.begin (), exportLIDs1Copy.end (),
1202  exportPIDs1Copy.begin ());
1203  typename ArrayView<LO>::iterator exportLIDs1_end = exportLIDs1Copy.end ();
1204  typename ArrayView<LO>::iterator exportPIDs1_end = exportPIDs1Copy.end ();
1205  merge2 (exportLIDs1_end, exportPIDs1_end,
1206  exportLIDs1Copy.begin (), exportLIDs1_end,
1207  exportPIDs1Copy.begin (), exportPIDs1_end,
1208  project1st<LO, LO> ());
1209 
1210  Array<LO> exportLIDs2Copy (exportLIDs2.begin (), exportLIDs2.end ());
1211  Array<int> exportPIDs2Copy (exportLIDs2.begin (), exportLIDs2.end ());
1212  sort2 (exportLIDs2Copy.begin (), exportLIDs2Copy.end (),
1213  exportPIDs2Copy.begin ());
1214  typename ArrayView<LO>::iterator exportLIDs2_end = exportLIDs2Copy.end ();
1215  typename ArrayView<LO>::iterator exportPIDs2_end = exportPIDs2Copy.end ();
1216  merge2 (exportLIDs2_end, exportPIDs2_end,
1217  exportLIDs2Copy.begin (), exportLIDs2_end,
1218  exportPIDs2Copy.begin (), exportPIDs2_end,
1219  project1st<LO, LO> ());
1220 
1221  // Merge export (LID,PID) pairs. In this merge operation, the
1222  // LIDs are the "keys" and the PIDs their "values." We combine
1223  // the "values" (PIDs) in the pairs by replacement, rather than
1224  // by adding them together.
1225  keyValueMerge (exportLIDs1Copy.begin (), exportLIDs1Copy.end (),
1226  exportPIDs1Copy.begin (), exportPIDs1Copy.end (),
1227  exportLIDs2Copy.begin (), exportLIDs2Copy.end (),
1228  exportPIDs2Copy.begin (), exportPIDs2Copy.end (),
1229  std::back_inserter (exportLIDsUnion),
1230  std::back_inserter (exportPIDsUnion),
1232 
1233  // Resort the merged (LID,PID) pairs by PID.
1234  sort2 (exportPIDsUnion.begin (), exportPIDsUnion.end (),
1235  exportLIDsUnion.begin ());
1236 
1237  // Initialize the Distributor. Using createFromSends instead of
1238  // createFromRecvs avoids the initialization and use of a
1239  // temporary Distributor object.
1240  (void) distributor.createFromSends (exportPIDsUnion ().getConst ());
1241 #else // NOT TPETRA_IMPORT_SETUNION_USE_CREATE_FROM_SENDS
1242 
1243  // Call the Distributor's createFromRecvs() method to turn the
1244  // remote GIDs and their owning processes into a send-and-receive
1245  // communication plan. remoteGIDsUnion and remotePIDsUnion are
1246  // input; exportGIDsUnion and exportPIDsUnion are output arrays
1247  // which are allocated by createFromRecvs().
1248  distributor.createFromRecvs (remoteGIDsUnion().getConst (),
1249  remotePIDsUnion ().getConst (),
1250  exportGIDsUnion, exportPIDsUnion);
1251 
1252  // Find the (source Map) LIDs corresponding to the export GIDs.
1253  const size_type numExportIDsUnion = exportGIDsUnion.size ();
1254  exportLIDsUnion.resize (numExportIDsUnion);
1255  for (size_type k = 0; k < numExportIDsUnion; ++k) {
1256  exportLIDsUnion[k] = srcMap->getLocalElement (exportGIDsUnion[k]);
1257  }
1258 #endif // TPETRA_IMPORT_SETUNION_USE_CREATE_FROM_SENDS
1259 
1260 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1261  {
1262  std::ostringstream os;
1263  os << myRank << ": exportGIDsUnion: " << toString (exportGIDsUnion ()) << endl;
1264  os << myRank << ": exportPIDsUnion: " << toString (exportPIDsUnion ()) << endl;
1265  os << myRank << ": exportLIDsUnion: " << toString (exportLIDsUnion ()) << endl;
1266  cerr << os.str ();
1267  }
1268  comm->barrier ();
1269  cerr << "Creating union Import" << endl;
1270 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1271 
1272  // Create and return the union Import.
1273  RCP<const import_type> unionImport =
1274  rcp (new import_type (srcMap, unionTgtMap,
1275  as<size_t> (numSameIDsUnion),
1276  permuteToLIDsUnion, permuteFromLIDsUnion,
1277  remoteLIDsUnion, exportLIDsUnion,
1278  exportPIDsUnion, distributor, this->out_));
1279 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1280  comm->barrier ();
1281  cerr << "Created union Import; done!" << endl;
1282 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1283 
1284  return unionImport;
1285  }
1286 
1287 
1288 
1289  template <class LocalOrdinal, class GlobalOrdinal, class Node>
1290  Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> >
1292  setUnion () const
1293  {
1294  using Teuchos::Array;
1295  using Teuchos::ArrayView;
1296  using Teuchos::as;
1297  using Teuchos::Comm;
1298  using Teuchos::RCP;
1299  using Teuchos::rcp;
1300  using Teuchos::outArg;
1301  using Teuchos::REDUCE_MIN;
1302  using Teuchos::reduceAll;
1303  typedef LocalOrdinal LO;
1304  typedef GlobalOrdinal GO;
1305  Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > unionImport;
1306  RCP<const map_type> srcMap = this->getSourceMap ();
1307  RCP<const map_type> tgtMap = this->getTargetMap ();
1308  RCP<const Comm<int> > comm = srcMap->getComm ();
1309 
1310 #ifdef HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1311  const int myRank = comm->getRank ();
1312 #endif // HAVE_TPETRA_IMPORT_SETUNION_EXTRA_DEBUG_OUTPUT
1313 
1314  ArrayView<const GO> srcGIDs = srcMap->getNodeElementList ();
1315  ArrayView<const GO> tgtGIDs = tgtMap->getNodeElementList ();
1316 
1317  // All elements in srcMap will be in the "new" target map, so...
1318  size_t numSameIDsNew = srcMap->getNodeNumElements();
1319  size_t numRemoteIDsNew = getNumRemoteIDs();
1320  Array<LO> permuteToLIDsNew, permuteFromLIDsNew; // empty on purpose
1321 
1322  // Grab some old data
1323  ArrayView<const LO> remoteLIDsOld = getRemoteLIDs();
1324  ArrayView<const LO> exportLIDsOld = getExportLIDs();
1325 
1326  // Build up the new map (same part)
1327  Array<GO> GIDs(numSameIDsNew + numRemoteIDsNew);
1328  for(size_t i=0; i<numSameIDsNew; i++)
1329  GIDs[i] = srcGIDs[i];
1330 
1331  // Build up the new map (remote part) and remotes list
1332  Array<LO> remoteLIDsNew(numRemoteIDsNew);
1333  for(size_t i=0; i<numRemoteIDsNew; i++) {
1334  GIDs[numSameIDsNew + i] = tgtGIDs[remoteLIDsOld[i]];
1335  remoteLIDsNew[i] = numSameIDsNew+i;
1336  }
1337 
1338  // Build the new target map
1339  GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
1340  RCP<const map_type> targetMapNew = rcp(new map_type(GO_INVALID,GIDs,tgtMap->getIndexBase(),tgtMap->getComm(),tgtMap->getNode()));
1341 
1342  // Exports are trivial (since the sourcemap doesn't change)
1343  Array<int> exportPIDsnew(getExportPIDs());
1344  Array<LO> exportLIDsnew(getExportLIDs());
1345 
1346  // Copy the Distributor (due to how the Import constructor works)
1348 
1349  // Build the importer
1350  unionImport = rcp(new Import<LocalOrdinal, GlobalOrdinal, Node>(srcMap,targetMapNew,numSameIDsNew,permuteToLIDsNew,permuteFromLIDsNew,
1351  remoteLIDsNew,exportLIDsnew,exportPIDsnew,D));
1352 
1353  return unionImport;
1354  }
1355 
1356 
1357 
1358 
1359  template <class LocalOrdinal, class GlobalOrdinal, class Node>
1360  Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> >
1362  createRemoteOnlyImport(const Teuchos::RCP<const map_type>& remoteTarget) const
1363  {
1364  using Teuchos::rcp;
1365  typedef Import<LocalOrdinal,GlobalOrdinal,Node> import_type;
1366 
1367  const size_t NumRemotes = getNumRemoteIDs ();
1368  TEUCHOS_TEST_FOR_EXCEPTION(
1369  NumRemotes != remoteTarget->getNodeNumElements (),
1370  std::runtime_error, "Tpetra::createRemoteOnlyImport: "
1371  "remoteTarget map ID count doesn't match.");
1372 
1373  // Compute the new Remote LIDs
1374  Teuchos::ArrayView<const LocalOrdinal> oldRemoteLIDs = getRemoteLIDs ();
1375  Teuchos::Array<LocalOrdinal> newRemoteLIDs (NumRemotes);
1376  for (size_t i = 0; i < NumRemotes; ++i) {
1377  newRemoteLIDs[i] = remoteTarget->getLocalElement (getTargetMap ()->getGlobalElement (oldRemoteLIDs[i]));
1378  // Now we make sure these guys are in sorted order (AztecOO-ML ordering)
1379  TEUCHOS_TEST_FOR_EXCEPTION(
1380  i > 0 && newRemoteLIDs[i] < newRemoteLIDs[i-1],
1381  std::runtime_error, "Tpetra::createRemoteOnlyImport: "
1382  "this and remoteTarget order don't match.");
1383  }
1384 
1385  // Copy ExportPIDs and such
1386  // NOTE: Be careful: The Import constructor we use does a "swap"
1387  // for most of the LID/PID lists and the Distributor, meaning it
1388  // ruins the existing object if we pass things in directly. Hence
1389  // we copy them first.
1390  Teuchos::Array<int> newExportPIDs (getExportPIDs ());
1391  Teuchos::Array<LocalOrdinal> newExportLIDs (getExportLIDs ());
1392  Teuchos::Array<LocalOrdinal> dummy;
1393  Distributor newDistor (getDistributor ());
1394 
1395  return rcp (new import_type (getSourceMap (), remoteTarget,
1396  static_cast<size_t> (0), dummy, dummy,
1397  newRemoteLIDs, newExportLIDs,
1398  newExportPIDs, newDistor));
1399  }
1400 
1401 } // namespace Tpetra
1402 
1403 #define TPETRA_IMPORT_CLASS_INSTANT(LO, GO, NODE) \
1404  \
1405  template class Import< LO , GO , NODE >;
1406 
1407 // Explicit instantiation macro.
1408 // Only invoke this when in the Tpetra namespace.
1409 // Most users do not need to use this.
1410 //
1411 // LO: The local ordinal type.
1412 // GO: The global ordinal type.
1413 // NODE: The Kokkos Node type.
1414 #define TPETRA_IMPORT_INSTANT(LO, GO, NODE) \
1415  TPETRA_IMPORT_CLASS_INSTANT(LO, GO, NODE)
1416 
1417 #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:65
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
::Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > map_type
The specialization of Map used by this class.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
ArrayView< const int > getExportPIDs() const
List of processes to which entries will be sent.
ArrayView< const LocalOrdinal > getRemoteLIDs() const
List of entries in the target Map to receive from other processes.
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.
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.
LookupStatus
Return status of Map remote index lookup (getRemoteIndexList()).
size_t getNumExportIDs() const
Number of entries that must be sent by the calling process to other processes.
bool isDistributed() const
Whether this Map is globally distributed or locally replicated.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
virtual ~Import()
Destructor.
Teuchos::RCP< const Import< LocalOrdinal, GlobalOrdinal, Node > > setUnion() const
Return the union of this Import this->getSourceMap()
size_t global_size_t
Global size_t object.
void merge2(IT1 &indResultOut, IT2 &valResultOut, IT1 indBeg, IT1 indEnd, IT2 valBeg, IT2 valEnd)
Merge values in place, additively, with the same index.
Import(const Teuchos::RCP< const map_type > &source, const Teuchos::RCP< const map_type > &target)
Construct an Import from the source and target Maps.
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...
ArrayView< const LocalOrdinal > getPermuteFromLIDs() const
List of local IDs in the source Map that are permuted.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
#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...
Import< LocalOrdinal, GlobalOrdinal, Node > & operator=(const Import< LocalOrdinal, GlobalOrdinal, Node > &Source)
Assignment operator.
Implementation detail of Import and Export.
Sets up and executes a communication plan for a Tpetra DistObject.
size_t createFromSends(const ArrayView< const int > &exportNodeIDs)
Set up Distributor using list of process ranks to which this process will send.
LookupStatus getRemoteIndexList(const Teuchos::ArrayView< const GlobalOrdinal > &GIDList, const Teuchos::ArrayView< int > &nodeIDList, const Teuchos::ArrayView< LocalOrdinal > &LIDList) const
Return the process ranks and corresponding local indices for the given global indices.
ArrayView< const LocalOrdinal > getPermuteToLIDs() const
List of local IDs in the target Map that are permuted.
ArrayView< const LocalOrdinal > getExportLIDs() const
List of entries in the source Map that will be sent to other processes.
LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const
The local index corresponding to the given global index.
Teuchos::ArrayView< const GlobalOrdinal > getNodeElementList() const
Return a view of the global indices owned by this process.
void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2)
Sort the first array, and apply the resulting permutation to the second array.
Describes a parallel distribution of objects over processes.
Internal functions and macros designed for use with Tpetra::Import and Tpetra::Export objects...
Stand-alone utility functions and macros.
size_t getNumPermuteIDs() const
Number of IDs to permute but not to communicate.
size_t getNumSameIDs() const
Number of initial identical IDs.
Teuchos::RCP< const map_type > getSourceMap() const
The Source Map used to construct this Import object.
Distributor & getDistributor() const
The Distributor that this Import object uses to move data.
virtual void print(std::ostream &os) const
Print the Import&#39;s data to the given output stream.
size_t getNumRemoteIDs() const
Number of entries not on the calling process.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream.
void createFromRecvs(const ArrayView< const Ordinal > &remoteIDs, const ArrayView< const int > &remoteNodeIDs, Array< Ordinal > &exportIDs, Array< int > &exportNodeIDs)
Set up Distributor using list of process ranks from which to receive.
Teuchos::RCP< const map_type > getTargetMap() const
The Target Map used to construct this Import object.
Binary function that returns its first argument.