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