Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Tpetra_Details_makeOptimizedColMap.hpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Tpetra: Templated Linear Algebra Services Package
6// Copyright (2008) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42*/
43
61
62#include "Tpetra_Map.hpp"
63#include "Tpetra_Import.hpp"
64#include "Tpetra_Util.hpp"
66#include "Teuchos_FancyOStream.hpp"
67#include <memory>
68
69namespace Tpetra {
70namespace Details {
71
78 template<class MapType>
79 class OptColMap {
80 public:
81 using local_ordinal_type = typename MapType::local_ordinal_type;
82 using global_ordinal_type = typename MapType::global_ordinal_type;
83 using node_type = typename MapType::node_type;
86
87 static Teuchos::RCP<const map_type>
88 makeOptColMap (std::ostream& errStream,
89 bool& lclErr,
90 const map_type& domMap,
91 const map_type& colMap,
92 const import_type* /* oldImport */)
93 {
94 using ::Tpetra::Details::Behavior;
95 using Teuchos::Array;
96 using Teuchos::ArrayView;
97 using Teuchos::FancyOStream;
98 using Teuchos::getFancyOStream;
99 using Teuchos::RCP;
100 using Teuchos::rcp;
101 using Teuchos::rcpFromRef;
102 using std::endl;
103 using LO = local_ordinal_type;
104 using GO = global_ordinal_type;
105 const char prefix[] = "Tpetra::Details::makeOptimizedColMap: ";
106
107 RCP<const Teuchos::Comm<int> > comm = colMap.getComm ();
108 std::ostream& err = errStream;
109
110 const bool verbose = Behavior::verbose ("Tpetra::Details::makeOptimizedColMap");
111
112 RCP<FancyOStream> outPtr = getFancyOStream (rcpFromRef (std::cerr));
113 TEUCHOS_TEST_FOR_EXCEPTION
114 (outPtr.is_null (), std::logic_error,
115 "outPtr is null; this should never happen!");
116 FancyOStream& out = *outPtr;
117 Teuchos::OSTab tab1 (out);
118
119 std::unique_ptr<std::string> verboseHeader;
120 if (verbose) {
121 std::ostringstream os;
122 const int myRank = comm->getRank ();
123 os << "Proc " << myRank << ": ";
124 verboseHeader = std::unique_ptr<std::string> (new std::string (os.str ()));
125 }
126 if (verbose) {
127 std::ostringstream os;
128 os << *verboseHeader << "Tpetra::Details::makeOptimizedColMap" << endl;
129 out << os.str ();
130 }
131
132 if (verbose) {
133 std::ostringstream os;
134 os << *verboseHeader << "Domain Map GIDs: [";
135 const LO domMapLclNumInds = static_cast<LO> (domMap.getLocalNumElements ());
136 for (LO lid = 0; lid < domMapLclNumInds; ++lid) {
137 const GO gid = domMap.getGlobalElement (lid);
138 os << gid;
139 if (lid + LO (1) < domMapLclNumInds) {
140 os << ", ";
141 }
142 }
143 os << "]" << endl;
144 out << os.str ();
145 }
146
147 const LO colMapLclNumInds = static_cast<LO> (colMap.getLocalNumElements ());
148
149 if (verbose) {
150 std::ostringstream os;
151 os << *verboseHeader << "Column Map GIDs: [";
152 for (LO lid = 0; lid < colMapLclNumInds; ++lid) {
153 const GO gid = colMap.getGlobalElement (lid);
154 os << gid;
155 if (lid + LO (1) < colMapLclNumInds) {
156 os << ", ";
157 }
158 }
159 os << "]" << endl;
160 out << os.str ();
161 }
162
163 // Count remote GIDs.
164 LO numOwnedGids = 0;
165 LO numRemoteGids = 0;
166 for (LO colMapLid = 0; colMapLid < colMapLclNumInds; ++colMapLid) {
167 const GO colMapGid = colMap.getGlobalElement (colMapLid);
168 if (domMap.isNodeGlobalElement (colMapGid)) {
169 ++numOwnedGids;
170 }
171 else {
172 ++numRemoteGids;
173 }
174 }
175
176 if (verbose) {
177 std::ostringstream os;
178 os << *verboseHeader << "- numOwnedGids: " << numOwnedGids << endl
179 << *verboseHeader << "- numRemoteGids: " << numRemoteGids << endl;
180 out << os.str ();
181 }
182
183 // Put all colMap GIDs on the calling process in a single array.
184 // Owned GIDs go in front, and remote GIDs at the end.
185 Array<GO> allGids (numOwnedGids + numRemoteGids);
186 ArrayView<GO> ownedGids = allGids.view (0, numOwnedGids);
187 ArrayView<GO> remoteGids = allGids.view (numOwnedGids, numRemoteGids);
188
189 // Fill ownedGids and remoteGids (and therefore allGids). We use
190 // two loops, one to count (above) and one to fill (here), in
191 // order to avoid dynamic memory allocation during the loop (in
192 // this case, lots of calls to push_back()). That will simplify
193 // use of Kokkos to parallelize these loops later.
194 LO ownedPos = 0;
195 LO remotePos = 0;
196 for (LO colMapLid = 0; colMapLid < colMapLclNumInds; ++colMapLid) {
197 const GO colMapGid = colMap.getGlobalElement (colMapLid);
198 if (domMap.isNodeGlobalElement (colMapGid)) {
199 ownedGids[ownedPos++] = colMapGid;
200 }
201 else {
202 remoteGids[remotePos++] = colMapGid;
203 }
204 }
205
206 // If, for some reason, the running count doesn't match the
207 // orignal count, fill in any remaining GID spots with an
208 // obviously invalid value. We don't want to stop yet, because
209 // other processes might not have noticed this error; Map
210 // construction is a collective, so we can't stop now.
211 if (ownedPos != numOwnedGids) {
212 lclErr = true;
213 err << prefix << "On Process " << comm->getRank () << ", ownedPos = "
214 << ownedPos << " != numOwnedGids = " << numOwnedGids << endl;
215 for (LO colMapLid = ownedPos; colMapLid < numOwnedGids; ++colMapLid) {
216 ownedGids[colMapLid] = Teuchos::OrdinalTraits<GO>::invalid ();
217 }
218 }
219 if (remotePos != numRemoteGids) {
220 lclErr = true;
221 err << prefix << "On Process " << comm->getRank () << ", remotePos = "
222 << remotePos << " != numRemoteGids = " << numRemoteGids << endl;
223 for (LO colMapLid = remotePos; colMapLid < numRemoteGids; ++colMapLid) {
224 remoteGids[colMapLid] = Teuchos::OrdinalTraits<GO>::invalid ();
225 }
226 }
227
228 // Figure out what processes own what GIDs in the domain Map.
229 // Initialize the output array of remote PIDs with the "invalid
230 // process rank" -1, to help us test whether getRemoteIndexList
231 // did its job.
232 Array<int> remotePids (numRemoteGids, -1);
233 const LookupStatus lookupStatus =
234 domMap.getRemoteIndexList (remoteGids, remotePids ());
235
236 // If any process returns IDNotPresent, then at least one of the
237 // remote indices was not present in the domain Map. This means
238 // that the Import object cannot be constructed, because of
239 // incongruity between the column Map and domain Map. This
240 // means that either the column Map or domain Map, or both, is
241 // incorrect.
242 const bool getRemoteIndexListFailed = (lookupStatus == IDNotPresent);
243 if (getRemoteIndexListFailed) {
244 lclErr = true;
245 err << prefix << "On Process " << comm->getRank () << ", some indices "
246 "in the input colMap (the original column Map) are not in domMap (the "
247 "domain Map). Either these indices or the domain Map is invalid. "
248 "Likely cause: For a nonsquare matrix, you must give the domain and "
249 "range Maps as input to fillComplete." << endl;
250 }
251
252 // Check that getRemoteIndexList actually worked, by making sure
253 // that none of the remote PIDs are -1.
254 for (LO k = 0; k < numRemoteGids; ++k) {
255 bool foundInvalidPid = false;
256 if (remotePids[k] == -1) {
257 foundInvalidPid = true;
258 break;
259 }
260 if (foundInvalidPid) {
261 lclErr = true;
262 err << prefix << "On Process " << comm->getRank () << ", "
263 "getRemoteIndexList returned -1 for the process ranks of "
264 "one or more GIDs on this process." << endl;
265 }
266 }
267
268 if (verbose) {
269 std::ostringstream os;
270 os << *verboseHeader << "- Before sort2:" << endl
271 << *verboseHeader << "-- ownedGids: " << Teuchos::toString (ownedGids) << endl
272 << *verboseHeader << "-- remoteGids: " << Teuchos::toString (remoteGids) << endl
273 << *verboseHeader << "-- allGids: " << Teuchos::toString (allGids ()) << endl;
274 out << os.str ();
275 }
276 using Tpetra::sort2;
277 sort2 (remotePids.begin (), remotePids.end (), remoteGids.begin ());
278 if (verbose) {
279 std::ostringstream os;
280 os << *verboseHeader << "- After sort2:" << endl
281 << *verboseHeader << "-- ownedGids: " << Teuchos::toString (ownedGids) << endl
282 << *verboseHeader << "-- remoteGids: " << Teuchos::toString (remoteGids) << endl
283 << *verboseHeader << "-- allGids: " << Teuchos::toString (allGids ()) << endl;
284 out << os.str ();
285 }
286
287 auto optColMap = rcp (new map_type (colMap.getGlobalNumElements (),
288 allGids (),
289 colMap.getIndexBase (),
290 comm));
291 if (verbose) {
292 std::ostringstream os;
293 os << *verboseHeader << "Tpetra::Details::makeOptimizedColMap: Done" << endl;
294 out << os.str ();
295 }
296 return optColMap;
297 }
298
329 static std::pair<Teuchos::RCP<const map_type>,
330 Teuchos::RCP<import_type> >
331 makeOptColMapAndImport (std::ostream& errStream,
332 bool& lclErr,
333 const map_type& domMap,
334 const map_type& colMap,
335 const import_type* oldImport)
336 {
337 using Teuchos::RCP;
338 using Teuchos::rcp;
339
340 // mfh 15 May 2018: For now, just call makeOptColMap, and use
341 // the conventional two-Map (source and target) Import
342 // constructor.
343 RCP<const map_type> newColMap =
344 makeOptColMap (errStream, lclErr, domMap, colMap, oldImport);
345 RCP<import_type> imp (new import_type (rcp (new map_type (domMap)), newColMap));
346
347 // FIXME (mfh 06 Jul 2014) This constructor throws a runtime
348 // error, so I'm not using it for now.
349 //
350 // imp = rcp (new import_type (domMap, newColMap, remoteGids,
351 // remotePids (), remoteLids (),
352 // Teuchos::null, Teuchos::null));
353 return std::make_pair (newColMap, imp);
354 }
355 };
356
381 template<class MapType>
382 Teuchos::RCP<const MapType>
383 makeOptimizedColMap (std::ostream& errStream,
384 bool& lclErr,
385 const MapType& domMap,
386 const MapType& colMap,
387 const Tpetra::Import<
388 typename MapType::local_ordinal_type,
389 typename MapType::global_ordinal_type,
390 typename MapType::node_type>* oldImport = nullptr)
391 {
392 using map_type = ::Tpetra::Map<
393 typename MapType::local_ordinal_type,
394 typename MapType::global_ordinal_type,
395 typename MapType::node_type>;
396 using impl_type = OptColMap<map_type>;
397 auto mapPtr = impl_type::makeOptColMap (errStream, lclErr,
398 domMap, colMap, oldImport);
399 return mapPtr;
400 }
401
458 template<class MapType>
459 std::pair<Teuchos::RCP<const MapType>,
460 Teuchos::RCP<typename OptColMap<MapType>::import_type> >
461 makeOptimizedColMapAndImport (std::ostream& errStream,
462 bool& lclErr,
463 const MapType& domMap,
464 const MapType& colMap,
465 const typename OptColMap<MapType>::import_type* oldImport = nullptr)
466 {
467 using local_ordinal_type = typename MapType::local_ordinal_type;
468 using global_ordinal_type = typename MapType::global_ordinal_type;
469 using node_type = typename MapType::node_type;
471 using impl_type = OptColMap<map_type>;
472
473 auto mapAndImp = impl_type::makeOptColMapAndImport (errStream, lclErr, domMap, colMap, oldImport);
474 return std::make_pair (mapAndImp.first, mapAndImp.second);
475 }
476
477} // namespace Details
478} // namespace Tpetra
Declaration of Tpetra::Details::Behavior, a class that describes Tpetra's behavior.
Stand-alone utility functions and macros.
static bool verbose()
Whether Tpetra is in verbose mode.
Implementation detail of makeOptimizedColMap, and makeOptimizedColMapAndImport.
static std::pair< Teuchos::RCP< const map_type >, Teuchos::RCP< import_type > > makeOptColMapAndImport(std::ostream &errStream, bool &lclErr, const map_type &domMap, const map_type &colMap, const import_type *oldImport)
Return an optimized reordering of the given column Map. Optionally, recompute an Import from the inpu...
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
A parallel distribution of indices over processes.
global_ordinal_type getGlobalElement(local_ordinal_type localIndex) const
The global index corresponding to the given local index.
LookupStatus getRemoteIndexList(const Teuchos::ArrayView< const global_ordinal_type > &GIDList, const Teuchos::ArrayView< int > &nodeIDList, const Teuchos::ArrayView< local_ordinal_type > &LIDList) const
Return the process ranks and corresponding local indices for the given global indices.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
global_ordinal_type getIndexBase() const
The index base for this Map.
global_size_t getGlobalNumElements() const
The number of elements in this Map.
bool isNodeGlobalElement(global_ordinal_type globalIndex) const
Whether the given global index is owned by this Map on the calling process.
size_t getLocalNumElements() const
The number of elements belonging to the calling process.
Implementation details of Tpetra.
std::pair< Teuchos::RCP< const MapType >, Teuchos::RCP< typename OptColMap< MapType >::import_type > > makeOptimizedColMapAndImport(std::ostream &errStream, bool &lclErr, const MapType &domMap, const MapType &colMap, const typename OptColMap< MapType >::import_type *oldImport=nullptr)
Return an optimized reordering of the given column Map. Optionally, recompute an Import from the inpu...
Teuchos::RCP< const MapType > makeOptimizedColMap(std::ostream &errStream, bool &lclErr, const MapType &domMap, const MapType &colMap, const Tpetra::Import< typename MapType::local_ordinal_type, typename MapType::global_ordinal_type, typename MapType::node_type > *oldImport=nullptr)
Return an optimized reordering of the given column Map.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
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()).