Tpetra parallel linear algebra  Version of the Day
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>
65 
66 namespace Tpetra {
67 namespace Details {
68 
75  template<class MapType>
76  class OptColMap {
77  public:
78  typedef MapType map_type;
79  typedef typename MapType::local_ordinal_type local_ordinal_type;
80  typedef typename MapType::global_ordinal_type global_ordinal_type;
81  typedef typename MapType::node_type node_type;
82  typedef Import<local_ordinal_type,
83  global_ordinal_type,
84  node_type> import_type;
85 
120  static std::pair<map_type, Teuchos::RCP<import_type> >
121  make (std::ostream& errStream,
122  bool& lclErr,
123  const map_type& domMap,
124  const map_type& colMap,
125  const import_type* oldImport,
126  const bool makeImport)
127  {
128  using Teuchos::Array;
129  using Teuchos::ArrayView;
130  using Teuchos::RCP;
131  using Teuchos::rcp;
132  using std::endl;
133  typedef local_ordinal_type LO;
134  typedef global_ordinal_type GO;
135  const char prefix[] = "Tpetra::makeOptimizedColMapAndImport: ";
136  std::ostream& err = errStream;
137 
138  (void) oldImport; // We don't currently use this argument.
139 
140  RCP<const Teuchos::Comm<int> > comm = colMap.getComm ();
141  const LO colMapMinLid = colMap.getMinLocalIndex ();
142  const LO colMapMaxLid = colMap.getMaxLocalIndex ();
143 
144  // Count the numbers of GIDs in colMap that are in and not in
145  // domMap on the calling process. Check for zero indices on the
146  // calling process first, because if it's true, then we shouldn't
147  // trust [getMinLocalIndex(), getMaxLocalIndex()] to return a
148  // correct range.
149  LO numOwnedGids = 0;
150  LO numRemoteGids = 0;
151  if (colMap.getNodeNumElements () != 0) {
152  for (LO colMapLid = colMapMinLid; colMapLid <= colMapMaxLid; ++colMapLid) {
153  const GO colMapGid = colMap.getGlobalElement (colMapLid);
154  if (domMap.isNodeLocalElement (colMapGid)) {
155  ++numOwnedGids;
156  } else {
157  ++numRemoteGids;
158  }
159  }
160  }
161 
162  // Put all colMap GIDs on the calling process in a single array.
163  // Owned GIDs go in front, and remote GIDs at the end.
164  Array<GO> allGids (numOwnedGids + numRemoteGids);
165  ArrayView<GO> ownedGids = allGids.view (0, numOwnedGids);
166  ArrayView<GO> remoteGids = allGids.view (numOwnedGids, numRemoteGids);
167 
168  // Fill ownedGids and remoteGids (and therefore allGids). We use
169  // two loops, one to count (above) and one to fill (here), in
170  // order to avoid dynamic memory allocation during the loop (in
171  // this case, lots of calls to push_back()). That will simplify
172  // use of Kokkos to parallelize these loops later.
173  LO ownedPos = 0;
174  LO remotePos = 0;
175  if (colMap.getNodeNumElements () != 0) {
176  for (LO colMapLid = colMapMinLid; colMapLid <= colMapMaxLid; ++colMapLid) {
177  const GO colMapGid = colMap.getGlobalElement (colMapLid);
178  if (domMap.isNodeLocalElement (colMapGid)) {
179  ownedGids[ownedPos++] = colMapGid;
180  } else {
181  remoteGids[remotePos++] = colMapGid;
182  }
183  }
184  }
185 
186  // If, for some reason, the running count doesn't match the
187  // orignal count, fill in any remaining GID spots with an
188  // obviously invalid value. We don't want to stop yet, because
189  // other processes might not have noticed this error; Map
190  // construction is a collective, so we can't stop now.
191  if (ownedPos != numOwnedGids) {
192  lclErr = true;
193  err << prefix << "On Process " << comm->getRank () << ", ownedPos = "
194  << ownedPos << " != numOwnedGids = " << numOwnedGids << endl;
195  for (LO colMapLid = ownedPos; colMapLid < numOwnedGids; ++colMapLid) {
196  ownedGids[colMapLid] = Teuchos::OrdinalTraits<GO>::invalid ();
197  }
198  }
199  if (remotePos != numRemoteGids) {
200  lclErr = true;
201  err << prefix << "On Process " << comm->getRank () << ", remotePos = "
202  << remotePos << " != numRemoteGids = " << numRemoteGids << endl;
203  for (LO colMapLid = remotePos; colMapLid < numRemoteGids; ++colMapLid) {
204  remoteGids[colMapLid] = Teuchos::OrdinalTraits<GO>::invalid ();
205  }
206  }
207 
208  // Figure out what processes own what GIDs in the domain Map.
209  // Initialize the output array of remote PIDs with the "invalid
210  // process rank" -1, to help us test whether getRemoteIndexList
211  // did its job.
212  Array<int> remotePids (numRemoteGids, -1);
213  Array<LO> remoteLids;
214  if (makeImport) {
215  remoteLids.resize (numRemoteGids);
216  std::fill (remoteLids.begin (), remoteLids.end (),
217  Teuchos::OrdinalTraits<LO>::invalid ());
218  }
219  LookupStatus lookupStatus;
220  if (makeImport) {
221  lookupStatus = domMap.getRemoteIndexList (remoteGids, remotePids (),
222  remoteLids ());
223  } else {
224  lookupStatus = domMap.getRemoteIndexList (remoteGids, remotePids ());
225  }
226 
227  // If any process returns IDNotPresent, then at least one of the
228  // remote indices was not present in the domain Map. This means
229  // that the Import object cannot be constructed, because of
230  // incongruity between the column Map and domain Map. This means
231  // that either the column Map or domain Map, or both, is
232  // incorrect.
233  const bool getRemoteIndexListFailed = (lookupStatus == IDNotPresent);
234  if (getRemoteIndexListFailed) {
235  lclErr = true;
236  err << prefix << "On Process " << comm->getRank () << ", some indices "
237  "in the input colMap (the original column Map) are not in domMap (the "
238  "domain Map). Either these indices or the domain Map is invalid. "
239  "Likely cause: For a nonsquare matrix, you must give the domain and "
240  "range Maps as input to fillComplete." << endl;
241  }
242 
243  // Check that getRemoteIndexList actually worked, by making sure
244  // that none of the remote PIDs are -1.
245  for (LO k = 0; k < numRemoteGids; ++k) {
246  bool foundInvalidPid = false;
247  if (remotePids[k] == -1) {
248  foundInvalidPid = true;
249  break;
250  }
251  if (foundInvalidPid) {
252  lclErr = true;
253  err << prefix << "On Process " << comm->getRank () << ", "
254  "getRemoteIndexList returned -1 for the process ranks of "
255  "one or more GIDs on this process." << endl;
256  }
257  }
258 
259  // Sort incoming remote column Map indices so that all columns
260  // coming from a given remote process are contiguous. This means
261  // the Import's Distributor doesn't need to reorder data.
262  if (makeImport) {
263  sort2 (remotePids.begin (), remotePids.end (), remoteGids.begin ());
264  }
265  else {
266  sort3 (remotePids.begin (), remotePids.end (),
267  remoteGids.begin (),
268  remoteLids.begin ());
269  }
270  // Make the new column Map.
271  MapType newColMap (colMap.getGlobalNumElements (), allGids (),
272  colMap.getIndexBase (), comm, colMap.getNode ());
273  // Optionally, make the new Import object.
274  RCP<import_type> imp;
275  if (makeImport) {
276  imp = rcp (new import_type (rcp (new map_type (domMap)),
277  rcp (new map_type (newColMap))));
278  // FIXME (mfh 06 Jul 2014) This constructor throws a runtime
279  // error, so I'm not using it for now.
280  //
281  // imp = rcp (new import_type (domMap, newColMap, remoteGids,
282  // remotePids (), remoteLids (),
283  // Teuchos::null, Teuchos::null));
284  }
285  return std::make_pair (newColMap, imp);
286  }
287  };
288 
308  template<class MapType>
309  MapType
310  makeOptimizedColMap (std::ostream& errStream,
311  bool& lclErr,
312  const MapType& domMap,
313  const MapType& colMap)
314  {
315  typedef typename MapType::local_ordinal_type LO;
316  typedef typename MapType::global_ordinal_type GO;
317  typedef typename MapType::node_type NT;
318  typedef ::Tpetra::Import<LO, GO, NT> import_type;
319 
320  const bool makeImport = false;
321  std::pair<MapType, Teuchos::RCP<import_type> > ret =
322  OptColMap<MapType>::make (errStream, lclErr, domMap, colMap,
323  NULL, makeImport);
324  return ret.first;
325  }
326 
381  template<class MapType>
382  std::pair<MapType, Teuchos::RCP<typename OptColMap<MapType>::import_type> >
383  makeOptimizedColMapAndImport (std::ostream& errStream,
384  bool& lclErr,
385  const MapType& domMap,
386  const MapType& colMap,
387  const typename OptColMap<MapType>::import_type* oldImport,
388  const bool makeImport)
389  {
390  return OptColMap<MapType>::make (errStream, lclErr, domMap, colMap,
391  oldImport, makeImport);
392  }
393 
394 } // namespace Details
395 } // namespace Tpetra
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
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.
LookupStatus
Return status of Map remote index lookup (getRemoteIndexList()).
MapType makeOptimizedColMap(std::ostream &errStream, bool &lclErr, const MapType &domMap, const MapType &colMap)
Return an optimized reordering of the given column Map.
static std::pair< map_type, Teuchos::RCP< import_type > > make(std::ostream &errStream, bool &lclErr, const map_type &domMap, const map_type &colMap, const import_type *oldImport, const bool makeImport)
Return an optimized reordering of the given column Map. Optionally, recompute an Import from the inpu...
Implementation details of Tpetra.
Implementation detail of makeOptimizedColMap, and makeOptimizedColMapAndImport.
void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2)
Sort the first array, and apply the resulting permutation to the second array.
Stand-alone utility functions and macros.
std::pair< 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, const bool makeImport)
Return an optimized reordering of the given column Map. Optionally, recompute an Import from the inpu...