Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
69 namespace Tpetra {
70 namespace 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.getNodeNumElements ());
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.getNodeNumElements ());
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
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
LookupStatus
Return status of Map remote index lookup (getRemoteIndexList()).
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
Implementation detail of makeOptimizedColMap, and makeOptimizedColMapAndImport.
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...
global_ordinal_type getIndexBase() const
The index base for this Map.
global_ordinal_type getGlobalElement(local_ordinal_type localIndex) const
The global index corresponding to the given local index.
bool isNodeGlobalElement(global_ordinal_type globalIndex) const
Whether the given global index is owned by this Map on the calling process.
static bool verbose()
Whether Tpetra is in verbose mode.
void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2)
Sort the first array, and apply the resulting permutation to the second array.
size_t getNodeNumElements() const
The number of elements belonging to the calling process.
A parallel distribution of indices over processes.
Stand-alone utility functions and macros.
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.
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...
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.
global_size_t getGlobalNumElements() const
The number of elements in this Map.
Declaration of Tpetra::Details::Behavior, a class that describes Tpetra&#39;s behavior.