Xpetra Version of the Day
Loading...
Searching...
No Matches
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_Exceptions.hpp"
53#include "Xpetra_Map.hpp"
54#include "Xpetra_MapFactory.hpp"
55
56namespace Xpetra {
57
58#ifndef DOXYGEN_SHOULD_SKIP_THIS
59 // forward declaration of BlockedMap, needed to prevent circular inclusions
60 template<class LO, class GO, class N> class BlockedMap;
61#endif
62
70template <class LocalOrdinal,
71 class GlobalOrdinal,
73class MapUtils {
74#undef XPETRA_MAPUTILS_SHORT
76
77public:
78
94
95 // ToDo Resolve header issues to allow for using this routing in Xpetra::BlockedMap.
96
97 // merge submaps to global map
98 std::vector<GlobalOrdinal> gids;
99 for(size_t tt = 0; tt<subMaps.size(); ++tt) {
101 Teuchos::ArrayView< const GlobalOrdinal > subMapGids = subMap->getLocalElementList();
102 gids.insert(gids.end(), subMapGids.begin(), subMapGids.end());
103 }
104
105 const GlobalOrdinal INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
106 //std::sort(gids.begin(), gids.end());
107 //gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
108 Teuchos::ArrayView<GlobalOrdinal> gidsView(&gids[0], gids.size());
109 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());
110 return fullMap;
111 }
112
131 {
134 "Xpetra::MatrixUtils::shrinkMapGIDs: the non-overlapping map must not have more local ids than the overlapping map.")
135
138 "Xpetra::MatrixUtils::shrinkMapGIDs: the maximum GIDs of the overlapping and non-overlapping maps must be the same.")
139
140 RCP< const Teuchos::Comm<int> > comm = input.getComm();
141
142 // we expect input to be the potentially overlapping map associated with nonOvlInput as the non-overlapping
143 // map with the same GIDs over all processors (e.g. column map and domain map). We use the nonOvlInput map
144 // to determine which GIDs are owned by which processor.
145
146 // calculate offset for new global Ids
147 std::vector<int> myGIDs(comm->getSize(),0);
148 std::vector<int> numGIDs(comm->getSize(),0);
149 myGIDs[comm->getRank()] = nonOvlInput.getLocalNumElements();
150 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myGIDs[0],&numGIDs[0]);
151 size_t gidOffset = 0;
152 for(int p = 0; p < comm->getRank(); p++) gidOffset += numGIDs[p];
153
154 // we use nonOvlInput to assign the globally unique shrinked GIDs and communicate them to input.
155 std::map<const GlobalOrdinal, GlobalOrdinal> origGID2newGID;
156 for(size_t i = 0; i < nonOvlInput.getLocalNumElements(); i++)
157 {
158 origGID2newGID[nonOvlInput.getGlobalElement(i)] = Teuchos::as<GlobalOrdinal>(i) + Teuchos::as<GlobalOrdinal>(gidOffset);
159 }
160 // build an overlapping version of mySpecialMap
161 Teuchos::Array<GlobalOrdinal> ovlUnknownStatusGids;
162 Teuchos::Array<GlobalOrdinal> ovlFoundStatusGids;
163 // loop over global column map of A and find all GIDs where it is not sure, whether they are special or not
164 for(size_t i = 0; i<input.getLocalNumElements(); i++) {
165 GlobalOrdinal gcid = input.getGlobalElement(i);
166 if( nonOvlInput.isNodeGlobalElement(gcid) == false) {
167 ovlUnknownStatusGids.push_back(gcid);
168 }
169 }
170
171 // Communicate the number of DOFs on each processor
172 std::vector<int> myUnknownDofGIDs(comm->getSize(),0);
173 std::vector<int> numUnknownDofGIDs(comm->getSize(),0);
174 myUnknownDofGIDs[comm->getRank()] = ovlUnknownStatusGids.size();
175 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myUnknownDofGIDs[0],&numUnknownDofGIDs[0]);
176
177 // create array containing all DOF GIDs
178 size_t cntUnknownDofGIDs = 0;
179 for(int p = 0; p < comm->getSize(); p++) cntUnknownDofGIDs += numUnknownDofGIDs[p];
180 std::vector<GlobalOrdinal> lUnknownDofGIDs(cntUnknownDofGIDs,0); // local version to be filled
181 std::vector<GlobalOrdinal> gUnknownDofGIDs(cntUnknownDofGIDs,0); // global version after communication
182 // calculate the offset and fill chunk of memory with local data on each processor
183 size_t cntUnknownOffset = 0;
184 for(int p = 0; p < comm->getRank(); p++) cntUnknownOffset += numUnknownDofGIDs[p];
185 for(size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
186 lUnknownDofGIDs[k+cntUnknownOffset] = ovlUnknownStatusGids[k];
187 }
188 if(cntUnknownDofGIDs > 0) // only perform communication if there are unknown DOF GIDs
189 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lUnknownDofGIDs[0],&gUnknownDofGIDs[0]);
190 std::vector<GlobalOrdinal> lTranslatedDofGIDs(cntUnknownDofGIDs,0); // local version to be filled
191 std::vector<GlobalOrdinal> gTranslatedDofGIDs(cntUnknownDofGIDs,0); // global version after communication
192 // loop through all GIDs with unknown status
193 for(size_t k=0; k < gUnknownDofGIDs.size(); k++) {
194 GlobalOrdinal curgid = gUnknownDofGIDs[k];
195 if(nonOvlInput.isNodeGlobalElement(curgid)) {
196 lTranslatedDofGIDs[k] = origGID2newGID[curgid]; // curgid is in special map (on this processor)
197 }
198 }
199 if(cntUnknownDofGIDs > 0) // only perform communication if there are unknown DOF GIDs
200 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lTranslatedDofGIDs[0],&gTranslatedDofGIDs[0]);
201
202 for(size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
203 origGID2newGID[ovlUnknownStatusGids[k]] = gTranslatedDofGIDs[k+cntUnknownOffset];
204 }
205 Teuchos::Array<GlobalOrdinal> ovlDomainMapArray;
206 for(size_t i = 0; i<input.getLocalNumElements(); i++) {
207 GlobalOrdinal gcid = input.getGlobalElement(i);
208 ovlDomainMapArray.push_back(origGID2newGID[gcid]);
209 }
212 (nonOvlInput.lib(),Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),ovlDomainMapArray(),0,comm);
213 return ovlDomainMap;
214 }
215
235 const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& nonOvlReferenceInput) {
236 //TEUCHOS_TEST_FOR_EXCEPTION(nonOvlInput.getLocalNumElements() > input.getLocalNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::transformThyra2XpetraGIDs: the non-overlapping map must not have more local ids than the overlapping map.");
237 TEUCHOS_TEST_FOR_EXCEPTION(nonOvlInput.getLocalNumElements() != nonOvlReferenceInput.getLocalNumElements(), 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!");
238 //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());
239
240 RCP< const Teuchos::Comm<int> > comm = input.getComm();
241
242 // fill translation map as far as possible
243 std::map<const GlobalOrdinal, GlobalOrdinal> thyra2xpetraGID;
244 for(size_t i = 0; i < nonOvlInput.getLocalNumElements(); i++) {
245 thyra2xpetraGID[nonOvlInput.getGlobalElement(i)] =
246 nonOvlReferenceInput.getGlobalElement(i);
247 }
248
249 // find all GIDs of the overlapping Thyra map which are not owned by this proc
250 Teuchos::Array<GlobalOrdinal> ovlUnknownStatusGids;
251 // loop over global column map of A and find all GIDs where it is not sure, whether they are special or not
252 for(size_t i = 0; i<input.getLocalNumElements(); i++) {
253 GlobalOrdinal gcid = input.getGlobalElement(i);
254 if( nonOvlInput.isNodeGlobalElement(gcid) == false) {
255 ovlUnknownStatusGids.push_back(gcid);
256 }
257 }
258
259 // Communicate the number of DOFs on each processor
260 std::vector<int> myUnknownDofGIDs(comm->getSize(),0);
261 std::vector<int> numUnknownDofGIDs(comm->getSize(),0);
262 myUnknownDofGIDs[comm->getRank()] = ovlUnknownStatusGids.size();
263 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myUnknownDofGIDs[0],&numUnknownDofGIDs[0]);
264
265 // create array containing all DOF GIDs
266 size_t cntUnknownDofGIDs = 0;
267 for(int p = 0; p < comm->getSize(); p++) cntUnknownDofGIDs += numUnknownDofGIDs[p];
268 std::vector<GlobalOrdinal> lUnknownDofGIDs(cntUnknownDofGIDs,0); // local version to be filled
269 std::vector<GlobalOrdinal> gUnknownDofGIDs(cntUnknownDofGIDs,0); // global version after communication
270 // calculate the offset and fill chunk of memory with local data on each processor
271 size_t cntUnknownOffset = 0;
272 for(int p = 0; p < comm->getRank(); p++) cntUnknownOffset += numUnknownDofGIDs[p];
273 for(size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
274 lUnknownDofGIDs[k+cntUnknownOffset] = ovlUnknownStatusGids[k];
275 }
276 if(cntUnknownDofGIDs > 0) // only perform communication if there are unknown DOF GIDs
277 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lUnknownDofGIDs[0],&gUnknownDofGIDs[0]);
278 std::vector<GlobalOrdinal> lTranslatedDofGIDs(cntUnknownDofGIDs,0); // local version to be filled
279 std::vector<GlobalOrdinal> gTranslatedDofGIDs(cntUnknownDofGIDs,0); // global version after communication
280 // loop through all GIDs with unknown status
281 for(size_t k=0; k < gUnknownDofGIDs.size(); k++) {
282 GlobalOrdinal curgid = gUnknownDofGIDs[k];
283 if(nonOvlInput.isNodeGlobalElement(curgid)) {
284 lTranslatedDofGIDs[k] = thyra2xpetraGID[curgid];
285 }
286 }
287 if(cntUnknownDofGIDs > 0) // only perform communication if there are unknown DOF GIDs
288 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lTranslatedDofGIDs[0],&gTranslatedDofGIDs[0]);
289
290 for(size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
291 thyra2xpetraGID[ovlUnknownStatusGids[k]] = gTranslatedDofGIDs[k+cntUnknownOffset];
292 }
293 Teuchos::Array<GlobalOrdinal> ovlDomainMapArray;
294 for(size_t i = 0; i<input.getLocalNumElements(); i++) {
295 GlobalOrdinal gcid = input.getGlobalElement(i);
296 ovlDomainMapArray.push_back(thyra2xpetraGID[gcid]);
297 }
300 (nonOvlInput.lib(),Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),ovlDomainMapArray(),0,comm);
301
302 TEUCHOS_TEST_FOR_EXCEPTION(input.getLocalNumElements() != ovlDomainMap->getLocalNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::transformThyra2XpetraGIDs: the number of local Thyra reference GIDs (overlapping) and local Xpetra GIDs (overlapping) must be the same!");
303 //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.");
304
305 return ovlDomainMap;
306 }
307
315 const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& input, GlobalOrdinal offset) {
316
318 RCP< const Teuchos::Comm<int> > comm = input.getComm();
319
320 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rcpInput = Teuchos::rcpFromRef(input);
321
322 // check whether input map is a BlockedMap or a standard Map
323 RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> > rcpBlockedInput = Teuchos::rcp_dynamic_cast<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpInput);
324 if(rcpBlockedInput.is_null() == true) {
325
326 // create a new map with Xpetra GIDs (may start not from GID = 0)
327 std::vector<GlobalOrdinal> gids;
328 for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(rcpInput->getLocalNumElements()); ++l) {
329 GlobalOrdinal gid = rcpInput->getGlobalElement(l) + offset;
330 gids.push_back(gid);
331 }
332 Teuchos::ArrayView<GO> gidsView(&gids[0], gids.size());
333 RCP<Map> fullMap = MapFactory::Build(rcpInput->lib(), INVALID, gidsView, rcpInput->getIndexBase(), comm);
334 return fullMap;
335 }
336
337 // SPECIAL CASE: input is a blocked map
338 // we have to recursively call this routine to get a BlockedMap containing (unique) Xpetra style GIDs
339
340 size_t numMaps = rcpBlockedInput->getNumMaps();
341
342 // first calucale GID offsets in submaps
343 // we need that for generating Xpetra GIDs
344 std::vector<GlobalOrdinal> gidOffsets(numMaps,0);
345 for(size_t i = 1; i < numMaps; ++i) {
346 gidOffsets[i] = rcpBlockedInput->getMap(i-1,true)->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
347 }
348
349 std::vector<RCP<const Map> > mapsXpetra(rcpBlockedInput->getNumMaps(), Teuchos::null);
350 std::vector<RCP<const Map> > mapsThyra (rcpBlockedInput->getNumMaps(), Teuchos::null);
351 for (size_t b = 0; b < rcpBlockedInput->getNumMaps(); ++b){
352 // extract sub map with Thyra style gids
353 // this can be an underlying Map or BlockedMap object
354 RCP<const Map> subMapThyra = rcpBlockedInput->getMap(b,true);
355 RCP<const Map> subMapXpetra = MapUtils::transformThyra2XpetraGIDs(*subMapThyra, gidOffsets[b] + offset); // recursive call
356 mapsXpetra[b] = subMapXpetra; // map can be of type Map or BlockedMap
357 mapsThyra[b] = subMapThyra; // map can be of type Map or BlockedMap
358 }
359
361 return resultMap;
362 }
363
364};
365
366} // end namespace Xpetra
367
368#define XPETRA_MAPUTILS_SHORT
369
370#endif // PACKAGES_XPETRA_SUP_MAP_UTILS_HPP_
iterator end() const
iterator begin() const
size_type size() const
void push_back(const value_type &x)
bool is_null() const
Exception throws to report incompatible objects (like maps).
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)
Map constructor with Xpetra-defined contiguous uniform distribution.
Xpetra utility class for common map-related routines.
static Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > transformThyra2XpetraGIDs(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &input, GlobalOrdinal offset)
replace set of global ids by new global ids
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
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.
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.
virtual UnderlyingLib lib() const =0
Get the library used by this object (Tpetra or Epetra?)
virtual bool isNodeGlobalElement(GlobalOrdinal globalIndex) const =0
Whether the given global index is valid for this Map on this process.
virtual GlobalOrdinal getGlobalElement(LocalOrdinal localIndex) const =0
The global index corresponding to the given local index.
virtual GlobalOrdinal getMaxAllGlobalIndex() const =0
The maximum global index over all processes in the communicator.
virtual size_t getLocalNumElements() const =0
The number of elements belonging to the calling process.
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const =0
Get this Map's Comm object.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra namespace