Xpetra_MapUtils.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
6 // Copyright 2012 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
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 #ifndef PACKAGES_XPETRA_SUP_MAP_UTILS_HPP_
48 #define PACKAGES_XPETRA_SUP_MAP_UTILS_HPP_
49 
50 #include "Xpetra_ConfigDefs.hpp"
51 
52 #include "Xpetra_Map.hpp"
53 
54 namespace Xpetra {
55 
63 template <class LocalOrdinal,
64  class GlobalOrdinal,
65  class Node>
66 class MapUtils {
67 #undef XPETRA_MAPUTILS_SHORT
69 
70 public:
71 
86  static Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > concatenateMaps(const std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > & subMaps) {
87 
88  // merge submaps to global map
89  std::vector<GlobalOrdinal> gids;
90  for(size_t tt = 0; tt<subMaps.size(); ++tt) {
91  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > subMap = subMaps[tt];
92  for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(subMap->getNodeNumElements()); ++l) {
93  GlobalOrdinal gid = subMap->getGlobalElement(l);
94  gids.push_back(gid);
95  }
96  }
97 
98  const GlobalOrdinal INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
99  //std::sort(gids.begin(), gids.end());
100  //gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
101  Teuchos::ArrayView<GlobalOrdinal> gidsView(&gids[0], gids.size());
102  Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullMap = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(subMaps[0]->lib(), INVALID, gidsView, subMaps[0]->getIndexBase(), subMaps[0]->getComm());
103  return fullMap;
104  }
105 
121  static Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > shrinkMapGIDs(
124  TEUCHOS_TEST_FOR_EXCEPTION(nonOvlInput.getNodeNumElements() > input.getNodeNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::shrinkMapGIDs: the non-overlapping map must not have more local ids than the overlapping map.")
125  TEUCHOS_TEST_FOR_EXCEPTION(nonOvlInput.getMaxAllGlobalIndex() != input.getMaxAllGlobalIndex(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::shrinkMapGIDs: the maximum GIDs of the overlapping and non-overlapping maps must be the same.")
126 
127  RCP< const Teuchos::Comm<int> > comm = input.getComm();
128 
129  // we expect input to be the potentially overlapping map associated with nonOvlInput as the non-overlapping
130  // map with the same GIDs over all processors (e.g. column map and domain map). We use the nonOvlInput map
131  // to determine which GIDs are owned by which processor.
132 
133  // calculate offset for new global Ids
134  std::vector<int> myGIDs(comm->getSize(),0);
135  std::vector<int> numGIDs(comm->getSize(),0);
136  myGIDs[comm->getRank()] = nonOvlInput.getNodeNumElements();
137  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myGIDs[0],&numGIDs[0]);
138  size_t gidOffset = 0;
139  for(int p = 0; p < comm->getRank(); p++) gidOffset += numGIDs[p];
140 
141  // we use nonOvlInput to assign the globally unique shrinked GIDs and communicate them to input.
142  std::map<const GlobalOrdinal, GlobalOrdinal> origGID2newGID;
143  for(size_t i = 0; i < nonOvlInput.getNodeNumElements(); i++) {
144  origGID2newGID[nonOvlInput.getGlobalElement(i)] = Teuchos::as<GlobalOrdinal>(i) + Teuchos::as<GlobalOrdinal>(gidOffset);
145  }
146  // build an overlapping version of mySpecialMap
147  Teuchos::Array<GlobalOrdinal> ovlUnknownStatusGids;
148  Teuchos::Array<GlobalOrdinal> ovlFoundStatusGids;
149  // loop over global column map of A and find all GIDs where it is not sure, whether they are special or not
150  for(size_t i = 0; i<input.getNodeNumElements(); i++) {
151  GlobalOrdinal gcid = input.getGlobalElement(i);
152  if( nonOvlInput.isNodeGlobalElement(gcid) == false) {
153  ovlUnknownStatusGids.push_back(gcid);
154  }
155  }
156 
157  // Communicate the number of DOFs on each processor
158  std::vector<int> myUnknownDofGIDs(comm->getSize(),0);
159  std::vector<int> numUnknownDofGIDs(comm->getSize(),0);
160  myUnknownDofGIDs[comm->getRank()] = ovlUnknownStatusGids.size();
161  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myUnknownDofGIDs[0],&numUnknownDofGIDs[0]);
162 
163  // create array containing all DOF GIDs
164  size_t cntUnknownDofGIDs = 0;
165  for(int p = 0; p < comm->getSize(); p++) cntUnknownDofGIDs += numUnknownDofGIDs[p];
166  std::vector<GlobalOrdinal> lUnknownDofGIDs(cntUnknownDofGIDs,0); // local version to be filled
167  std::vector<GlobalOrdinal> gUnknownDofGIDs(cntUnknownDofGIDs,0); // global version after communication
168  // calculate the offset and fill chunk of memory with local data on each processor
169  size_t cntUnknownOffset = 0;
170  for(int p = 0; p < comm->getRank(); p++) cntUnknownOffset += numUnknownDofGIDs[p];
171  for(size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
172  lUnknownDofGIDs[k+cntUnknownOffset] = ovlUnknownStatusGids[k];
173  }
174  if(cntUnknownDofGIDs > 0) // only perform communication if there are unknown DOF GIDs
175  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lUnknownDofGIDs[0],&gUnknownDofGIDs[0]);
176  std::vector<GlobalOrdinal> lTranslatedDofGIDs(cntUnknownDofGIDs,0); // local version to be filled
177  std::vector<GlobalOrdinal> gTranslatedDofGIDs(cntUnknownDofGIDs,0); // global version after communication
178  // loop through all GIDs with unknown status
179  for(size_t k=0; k < gUnknownDofGIDs.size(); k++) {
180  GlobalOrdinal curgid = gUnknownDofGIDs[k];
181  if(nonOvlInput.isNodeGlobalElement(curgid)) {
182  lTranslatedDofGIDs[k] = origGID2newGID[curgid]; // curgid is in special map (on this processor)
183  }
184  }
185  if(cntUnknownDofGIDs > 0) // only perform communication if there are unknown DOF GIDs
186  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lTranslatedDofGIDs[0],&gTranslatedDofGIDs[0]);
187 
188  for(size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
189  origGID2newGID[ovlUnknownStatusGids[k]] = gTranslatedDofGIDs[k+cntUnknownOffset];
190  }
191  Teuchos::Array<GlobalOrdinal> ovlDomainMapArray;
192  for(size_t i = 0; i<input.getNodeNumElements(); i++) {
193  GlobalOrdinal gcid = input.getGlobalElement(i);
194  ovlDomainMapArray.push_back(origGID2newGID[gcid]);
195  }
196  RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > ovlDomainMap =
198  (nonOvlInput.lib(),Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),ovlDomainMapArray(),0,comm);
199  return ovlDomainMap;
200  }
201 
218  static Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > transformThyra2XpetraGIDs(
221  const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& nonOvlReferenceInput) {
222  //TEUCHOS_TEST_FOR_EXCEPTION(nonOvlInput.getNodeNumElements() > input.getNodeNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::transformThyra2XpetraGIDs: the non-overlapping map must not have more local ids than the overlapping map.");
223  TEUCHOS_TEST_FOR_EXCEPTION(nonOvlInput.getNodeNumElements() != nonOvlReferenceInput.getNodeNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::transformThyra2XpetraGIDs: the number of local Xpetra reference GIDs and local Thyra GIDs of the non-overlapping maps must be the same!");
224  //TEUCHOS_TEST_FOR_EXCEPTION(nonOvlInput.getMaxAllGlobalIndex() != input.getMaxAllGlobalIndex(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::transformThyra2XpetraGIDs: the maximum GIDs of the overlapping and non-overlapping maps must be the same. nonOvlInput.getMaxAllGlobalIndex() = " << nonOvlInput.getMaxAllGlobalIndex() << " ovlInput.getMaxAllGlobalIndex() = " << input.getMaxAllGlobalIndex());
225 
226  RCP< const Teuchos::Comm<int> > comm = input.getComm();
227 
228  // fill translation map as far as possible
229  std::map<const GlobalOrdinal, GlobalOrdinal> thyra2xpetraGID;
230  for(size_t i = 0; i < nonOvlInput.getNodeNumElements(); i++) {
231  thyra2xpetraGID[nonOvlInput.getGlobalElement(i)] =
232  nonOvlReferenceInput.getGlobalElement(i);
233  }
234 
235  // find all GIDs of the overlapping Thyra map which are not owned by this proc
236  Teuchos::Array<GlobalOrdinal> ovlUnknownStatusGids;
237  // loop over global column map of A and find all GIDs where it is not sure, whether they are special or not
238  for(size_t i = 0; i<input.getNodeNumElements(); i++) {
239  GlobalOrdinal gcid = input.getGlobalElement(i);
240  if( nonOvlInput.isNodeGlobalElement(gcid) == false) {
241  ovlUnknownStatusGids.push_back(gcid);
242  }
243  }
244 
245  // Communicate the number of DOFs on each processor
246  std::vector<int> myUnknownDofGIDs(comm->getSize(),0);
247  std::vector<int> numUnknownDofGIDs(comm->getSize(),0);
248  myUnknownDofGIDs[comm->getRank()] = ovlUnknownStatusGids.size();
249  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myUnknownDofGIDs[0],&numUnknownDofGIDs[0]);
250 
251  // create array containing all DOF GIDs
252  size_t cntUnknownDofGIDs = 0;
253  for(int p = 0; p < comm->getSize(); p++) cntUnknownDofGIDs += numUnknownDofGIDs[p];
254  std::vector<GlobalOrdinal> lUnknownDofGIDs(cntUnknownDofGIDs,0); // local version to be filled
255  std::vector<GlobalOrdinal> gUnknownDofGIDs(cntUnknownDofGIDs,0); // global version after communication
256  // calculate the offset and fill chunk of memory with local data on each processor
257  size_t cntUnknownOffset = 0;
258  for(int p = 0; p < comm->getRank(); p++) cntUnknownOffset += numUnknownDofGIDs[p];
259  for(size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
260  lUnknownDofGIDs[k+cntUnknownOffset] = ovlUnknownStatusGids[k];
261  }
262  if(cntUnknownDofGIDs > 0) // only perform communication if there are unknown DOF GIDs
263  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lUnknownDofGIDs[0],&gUnknownDofGIDs[0]);
264  std::vector<GlobalOrdinal> lTranslatedDofGIDs(cntUnknownDofGIDs,0); // local version to be filled
265  std::vector<GlobalOrdinal> gTranslatedDofGIDs(cntUnknownDofGIDs,0); // global version after communication
266  // loop through all GIDs with unknown status
267  for(size_t k=0; k < gUnknownDofGIDs.size(); k++) {
268  GlobalOrdinal curgid = gUnknownDofGIDs[k];
269  if(nonOvlInput.isNodeGlobalElement(curgid)) {
270  lTranslatedDofGIDs[k] = thyra2xpetraGID[curgid];
271  }
272  }
273  if(cntUnknownDofGIDs > 0) // only perform communication if there are unknown DOF GIDs
274  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lTranslatedDofGIDs[0],&gTranslatedDofGIDs[0]);
275 
276  for(size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
277  thyra2xpetraGID[ovlUnknownStatusGids[k]] = gTranslatedDofGIDs[k+cntUnknownOffset];
278  }
279  Teuchos::Array<GlobalOrdinal> ovlDomainMapArray;
280  for(size_t i = 0; i<input.getNodeNumElements(); i++) {
281  GlobalOrdinal gcid = input.getGlobalElement(i);
282  ovlDomainMapArray.push_back(thyra2xpetraGID[gcid]);
283  }
284  RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > ovlDomainMap =
286  (nonOvlInput.lib(),Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),ovlDomainMapArray(),0,comm);
287 
288  TEUCHOS_TEST_FOR_EXCEPTION(input.getNodeNumElements() != ovlDomainMap->getNodeNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::transformThyra2XpetraGIDs: the number of local Thyra reference GIDs (overlapping) and local Xpetra GIDs (overlapping) must be the same!");
289  //TEUCHOS_TEST_FOR_EXCEPTION(nonOvlReferenceInput.getMaxAllGlobalIndex() != ovlDomainMap->getMaxAllGlobalIndex(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::transformThyra2XpetraGIDs: the maximum GIDs of the overlapping and non-overlapping Xpetra maps must be the same.");
290 
291  return ovlDomainMap;
292  }
293 
294 
295 };
296 
297 } // end namespace Xpetra
298 
299 #define XPETRA_MAPUTILS_SHORT
300 
301 #endif // PACKAGES_XPETRA_SUP_MAP_UTILS_HPP_
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const =0
Get this Map&#39;s Comm object.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &node=defaultArgNode())
Map constructor with Xpetra-defined contiguous uniform distribution.
Xpetra namespace
virtual size_t getNodeNumElements() const =0
The number of elements belonging to the calling process.
static Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > concatenateMaps(const std::vector< Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > > &subMaps)
Helper function to concatenate several maps.
virtual bool isNodeGlobalElement(GlobalOrdinal globalIndex) const =0
Whether the given global index is valid for this Map on this process.
static Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > transformThyra2XpetraGIDs(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &input, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlInput, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlReferenceInput)
replace set of global ids by new global ids
virtual GlobalOrdinal getMaxAllGlobalIndex() const =0
The maximum global index over all processes in the communicator.
virtual GlobalOrdinal getGlobalElement(LocalOrdinal localIndex) const =0
The global index corresponding to the given local index.
static Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > shrinkMapGIDs(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &input, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlInput)
Helper function to shrink the GIDs and generate a standard map whith GIDs starting at 0...
Exception throws to report incompatible objects (like maps).
virtual UnderlyingLib lib() const =0
Get the library used by this object (Tpetra or Epetra?)