43 #ifndef DOMI_MDMAP_HPP 44 #define DOMI_MDMAP_HPP 52 #include "Domi_ConfigDefs.hpp" 53 #include "Domi_DefaultNode.hpp" 54 #include "Domi_Utils.hpp" 55 #include "Domi_getValidParameters.hpp" 57 #include "Domi_MDComm.hpp" 58 #include "Domi_MDArray.hpp" 61 #include "Teuchos_Comm.hpp" 62 #include "Teuchos_DefaultComm.hpp" 63 #include "Teuchos_CommHelpers.hpp" 64 #include "Teuchos_Tuple.hpp" 67 #include "Kokkos_Core.hpp" 70 #include "Epetra_Map.h" 74 #include "Tpetra_Map.hpp" 75 #include "Tpetra_Util.hpp" 143 template<
class Node = DefaultNode::DefaultNodeType >
184 MDMap(
const Teuchos::RCP< const MDComm > mdComm,
185 const Teuchos::ArrayView< const dim_type > & dimensions,
186 const Teuchos::ArrayView< const int > & commPad =
187 Teuchos::ArrayView< const int >(),
188 const Teuchos::ArrayView< const int > & bndryPad =
189 Teuchos::ArrayView< const int >(),
190 const Teuchos::ArrayView< const int > & replicatedBoundary =
191 Teuchos::ArrayView< const int >(),
192 const Layout layout = DEFAULT_ORDER);
205 MDMap(Teuchos::ParameterList & plist);
219 MDMap(
const Teuchos::RCP<
const Teuchos::Comm< int > > teuchosComm,
220 Teuchos::ParameterList & plist);
234 MDMap(
const Teuchos::RCP< const MDComm > mdComm,
235 Teuchos::ParameterList & plist);
260 MDMap(
const Teuchos::RCP< const MDComm > mdComm,
261 const Teuchos::ArrayView< Slice > & myGlobalBounds,
262 const Teuchos::ArrayView< padding_type > & padding =
263 Teuchos::ArrayView< padding_type >(),
264 const Teuchos::ArrayView< const int > & replicatedBoundary =
265 Teuchos::ArrayView< const int >(),
266 const Layout layout = DEFAULT_ORDER);
331 const Teuchos::ArrayView< Slice > & slices,
332 const Teuchos::ArrayView< int > & bndryPad =
333 Teuchos::ArrayView< int >());
360 inline Teuchos::RCP< const MDComm >
getMDComm()
const;
376 inline Teuchos::RCP< const Teuchos::Comm< int > >
getTeuchosComm()
const;
483 bool withBndryPad=
false)
const;
494 bool withBndryPad=
false)
const;
510 bool withBndryPad=
false)
const;
525 bool withPad=
false)
const;
550 bool withPad=
false)
const;
659 bool isPad(
const Teuchos::ArrayView< dim_type > & index)
const;
667 bool isCommPad(
const Teuchos::ArrayView< dim_type > & index)
const;
675 bool isBndryPad(
const Teuchos::ArrayView< dim_type > & index)
const;
687 Teuchos::RCP< Node >
getNode()
const;
698 Teuchos::ArrayView< Teuchos::RCP< const Domi::MDMap< Node > > >
705 Teuchos::RCP< const Domi::MDMap< Node > >
getAxisMap(
int axis)
const;
723 Teuchos::RCP< const MDMap< Node > >
725 const dim_type trailingDim=0)
const;
737 Teuchos::RCP< const Epetra_Map >
738 getEpetraMap(
bool withCommPad=
true)
const;
750 Teuchos::RCP< const Epetra_Map >
751 getEpetraAxisMap(
int axis,
752 bool withCommPad=
true)
const;
767 template<
class LocalOrdinal,
768 class GlobalOrdinal = LocalOrdinal,
770 Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node2 > >
771 getTpetraMap(
bool withCommPad=
true)
const;
785 template<
class LocalOrdinal,
786 class GlobalOrdinal = LocalOrdinal,
788 Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node2 > >
789 getTpetraAxisMap(
int axis,
790 bool withCommPad=
true)
const;
803 Teuchos::Array< dim_type >
810 Teuchos::Array< dim_type >
824 getGlobalID(
const Teuchos::ArrayView< dim_type > & globalIndex)
const;
833 size_type
getLocalID(size_type globalID)
const;
840 getLocalID(
const Teuchos::ArrayView< dim_type > & localIndex)
const;
885 const int verbose = 0)
const;
905 void computeBounds();
908 Teuchos::RCP< const MDComm > _mdComm;
912 Teuchos::Array< dim_type > _globalDims;
916 Teuchos::Array< Slice > _globalBounds;
925 Teuchos::Array< Teuchos::Array< Slice > > _globalRankBounds;
930 Teuchos::Array< size_type > _globalStrides;
934 size_type _globalMin;
938 size_type _globalMax;
942 Teuchos::Array< dim_type > _localDims;
950 Teuchos::Array< Slice > _localBounds;
953 Teuchos::Array< size_type > _localStrides;
965 Teuchos::Array< int > _commPadSizes;
970 Teuchos::Array< padding_type > _pad;
973 Teuchos::Array< int > _bndryPadSizes;
981 Teuchos::Array< padding_type > _bndryPad;
988 Teuchos::Array< int > _replicatedBoundary;
997 mutable Teuchos::Array< Teuchos::RCP< const MDMap< Node > > > _axisMaps;
1000 Teuchos::RCP< Node > _node;
1007 mutable Teuchos::RCP< const Epetra_Map > _epetraMap;
1014 mutable Teuchos::RCP< const Epetra_Map > _epetraOwnMap;
1020 mutable Teuchos::Array< Teuchos::RCP< const Epetra_Map > > _epetraAxisMaps;
1026 mutable Teuchos::Array< Teuchos::RCP< const Epetra_Map > > _epetraAxisOwnMaps;
1035 template<
class Node >
1037 MDMap(
const Teuchos::RCP< const MDComm > mdComm,
1038 const Teuchos::ArrayView< const dim_type > & dimensions,
1039 const Teuchos::ArrayView< const int > & commPad,
1040 const Teuchos::ArrayView< const int > & bndryPad,
1041 const Teuchos::ArrayView< const int > & replicatedBoundary,
1042 const Layout layout) :
1044 _globalDims(mdComm->numDims()),
1046 _globalRankBounds(mdComm->numDims()),
1050 _localDims(mdComm->numDims(), 0),
1055 _commPadSizes(mdComm->numDims(), 0),
1057 _bndryPadSizes(mdComm->numDims(), 0),
1059 _replicatedBoundary(createArrayOfInts(mdComm->numDims(),
1060 replicatedBoundary)),
1064 int numDims = mdComm->numDims();
1067 TEUCHOS_TEST_FOR_EXCEPTION(
1070 "Size of dimensions does not match MDComm number of dimensions");
1074 for (
int axis = 0; axis <
numDims; ++axis)
1076 if (axis < commPad.size() ) _commPadSizes[ axis] = commPad[ axis];
1077 if (axis < bndryPad.size()) _bndryPadSizes[axis] = bndryPad[axis];
1078 if (_mdComm->isPeriodic(axis))
1079 _bndryPad.push_back(Teuchos::tuple(_commPadSizes[axis],
1080 _commPadSizes[axis]));
1082 _bndryPad.push_back(Teuchos::tuple(_bndryPadSizes[axis],
1083 _bndryPadSizes[axis]));
1084 _globalDims[axis] = dimensions[axis] + _bndryPad[axis][0] +
1089 lower = _bndryPadSizes[axis];
1091 lower = _commPadSizes[axis];
1093 upper = _bndryPadSizes[axis];
1095 upper = _commPadSizes[axis];
1096 _pad.push_back(Teuchos::tuple(lower, upper));
1100 _globalMax = computeSize(_globalDims());
1105 _localMax = computeSize(_localDims());
1108 _globalStrides = computeStrides< size_type, dim_type >(_globalDims, _layout);
1109 _localStrides = computeStrides< size_type, dim_type >(_localDims , _layout);
1114 template<
class Node >
1117 _mdComm(Teuchos::rcp(new
MDComm(plist))),
1120 _globalRankBounds(),
1133 _replicatedBoundary(),
1141 int numDims = _mdComm->numDims();
1144 Teuchos::Array< dim_type > dimensions =
1145 plist.get(
"dimensions", Teuchos::Array< dim_type >());
1146 TEUCHOS_TEST_FOR_EXCEPTION(
1149 "Size of dimensions does not match MDComm number of dimensions");
1153 int commPad = plist.get(
"communication pad size",
int(0));
1154 int bndryPad = plist.get(
"boundary pad size" ,
int(0));
1155 Teuchos::Array< int > commPads =
1156 plist.get(
"communication pad sizes", Teuchos::Array< int >());
1157 Teuchos::Array< int > bndryPads =
1158 plist.get(
"boundary pad sizes" , Teuchos::Array< int >());
1159 _commPadSizes.resize(
numDims);
1160 _bndryPadSizes.resize(
numDims);
1165 for (
int axis = 0; axis <
numDims; ++axis)
1167 if (axis < commPads.size() ) _commPadSizes[ axis] = commPads[ axis];
1168 else _commPadSizes[ axis] = commPad;
1169 if (axis < bndryPads.size()) _bndryPadSizes[axis] = bndryPads[axis];
1170 else _bndryPadSizes[axis] = bndryPad;
1171 if (_mdComm->isPeriodic(axis))
1172 _bndryPad.push_back(Teuchos::tuple(_commPadSizes[axis],
1173 _commPadSizes[axis]));
1175 _bndryPad.push_back(Teuchos::tuple(_bndryPadSizes[axis],
1176 _bndryPadSizes[axis]));
1177 _globalDims[axis] = dimensions[axis] + _bndryPad[axis][0] +
1182 lower = _bndryPadSizes[axis];
1184 lower = _commPadSizes[axis];
1186 upper = _bndryPadSizes[axis];
1188 upper = _commPadSizes[axis];
1189 _pad.push_back(Teuchos::tuple(lower, upper));
1193 _globalMax = computeSize(_globalDims());
1197 _globalRankBounds.resize(
numDims);
1200 _localMax = computeSize(_localDims());
1203 Teuchos::Array< int > repBndry = plist.get(
"replicated boundary",
1204 Teuchos::Array< int >());
1205 _replicatedBoundary = createArrayOfInts(
numDims, repBndry);
1208 std::string layout = plist.get(
"layout",
"DEFAULT");
1209 std::transform(layout.begin(), layout.end(), layout.begin(), ::toupper);
1210 if (layout ==
"C ORDER")
1212 else if (layout ==
"FORTRAN ORDER")
1213 _layout = FORTRAN_ORDER;
1214 else if (layout ==
"ROW MAJOR")
1215 _layout = ROW_MAJOR;
1216 else if (layout ==
"COLUMN MAJOR")
1217 _layout = COLUMN_MAJOR;
1218 else if (layout ==
"LAST INDEX FASTEST")
1219 _layout = LAST_INDEX_FASTEST;
1220 else if (layout ==
"FIRST INDEX FASTEST")
1221 _layout = FIRST_INDEX_FASTEST;
1223 _layout = DEFAULT_ORDER;
1226 _globalStrides = computeStrides< size_type, dim_type >(_globalDims, _layout);
1227 _localStrides = computeStrides< size_type, dim_type >(_localDims , _layout);
1232 template<
class Node >
1234 MDMap(Teuchos::RCP<
const Teuchos::Comm< int > > teuchosComm,
1235 Teuchos::ParameterList & plist) :
1236 _mdComm(Teuchos::rcp(new
MDComm(teuchosComm, plist))),
1239 _globalRankBounds(),
1252 _replicatedBoundary(),
1260 int numDims = _mdComm->numDims();
1263 Teuchos::Array< dim_type > dimensions =
1264 plist.get(
"dimensions", Teuchos::Array< dim_type >());
1265 TEUCHOS_TEST_FOR_EXCEPTION(
1268 "Size of dimensions does not match MDComm number of dimensions");
1272 int commPad = plist.get(
"communication pad size",
int(0));
1273 int bndryPad = plist.get(
"boundary pad size" ,
int(0));
1274 Teuchos::Array< int > commPads =
1275 plist.get(
"communication pad sizes", Teuchos::Array< int >());
1276 Teuchos::Array< int > bndryPads =
1277 plist.get(
"boundary pad sizes" , Teuchos::Array< int >());
1278 _commPadSizes.resize(
numDims);
1279 _bndryPadSizes.resize(
numDims);
1284 for (
int axis = 0; axis <
numDims; ++axis)
1286 if (axis < commPads.size() ) _commPadSizes[ axis] = commPads[ axis];
1287 else _commPadSizes[ axis] = commPad;
1288 if (axis < bndryPads.size()) _bndryPadSizes[axis] = bndryPads[axis];
1289 else _bndryPadSizes[axis] = bndryPad;
1290 if (_mdComm->isPeriodic(axis))
1291 _bndryPad.push_back(Teuchos::tuple(_commPadSizes[axis],
1292 _commPadSizes[axis]));
1294 _bndryPad.push_back(Teuchos::tuple(_bndryPadSizes[axis],
1295 _bndryPadSizes[axis]));
1296 _globalDims[axis] = dimensions[axis] + _bndryPad[axis][0] +
1301 lower = _bndryPadSizes[axis];
1303 lower = _commPadSizes[axis];
1305 upper = _bndryPadSizes[axis];
1307 upper = _commPadSizes[axis];
1308 _pad.push_back(Teuchos::tuple(lower, upper));
1317 _globalMax = computeSize(_globalDims());
1321 _globalRankBounds.resize(
numDims);
1324 _localMax = computeSize(_localDims());
1327 Teuchos::Array< int > repBndry = plist.get(
"replicated boundary",
1328 Teuchos::Array< int >());
1329 _replicatedBoundary = createArrayOfInts(
numDims, repBndry);
1332 std::string layout = plist.get(
"layout",
"DEFAULT");
1333 std::transform(layout.begin(), layout.end(), layout.begin(), ::toupper);
1334 if (layout ==
"C ORDER")
1336 else if (layout ==
"FORTRAN ORDER")
1337 _layout = FORTRAN_ORDER;
1338 else if (layout ==
"ROW MAJOR")
1339 _layout = ROW_MAJOR;
1340 else if (layout ==
"COLUMN MAJOR")
1341 _layout = COLUMN_MAJOR;
1342 else if (layout ==
"LAST INDEX FASTEST")
1343 _layout = LAST_INDEX_FASTEST;
1344 else if (layout ==
"FIRST INDEX FASTEST")
1345 _layout = FIRST_INDEX_FASTEST;
1347 _layout = DEFAULT_ORDER;
1350 _globalStrides = computeStrides< size_type, dim_type >(_globalDims, _layout);
1351 _localStrides = computeStrides< size_type, dim_type >(_localDims , _layout);
1356 template<
class Node >
1358 MDMap(Teuchos::RCP< const MDComm > mdComm,
1359 Teuchos::ParameterList & plist) :
1361 _globalDims(mdComm->numDims()),
1363 _globalRankBounds(mdComm->numDims()),
1367 _localDims(mdComm->numDims(), 0),
1372 _commPadSizes(mdComm->numDims(), 0),
1374 _bndryPadSizes(mdComm->numDims(), 0),
1376 _replicatedBoundary(),
1380 plist.validateParameters(*getValidParameters());
1383 int numDims = _mdComm->numDims();
1386 Teuchos::Array< dim_type > dimensions =
1387 plist.get(
"dimensions", Teuchos::Array< dim_type >());
1388 TEUCHOS_TEST_FOR_EXCEPTION(
1391 "Number of dimensions does not match MDComm number of dimensions");
1395 int commPad = plist.get(
"communication pad size",
int(0));
1396 int bndryPad = plist.get(
"boundary pad size" ,
int(0));
1397 Teuchos::Array< int > commPads =
1398 plist.get(
"communication pad sizes", Teuchos::Array< int >());
1399 Teuchos::Array< int > bndryPads =
1400 plist.get(
"boundary pad sizes" , Teuchos::Array< int >());
1401 _commPadSizes.resize(
numDims);
1402 _bndryPadSizes.resize(
numDims);
1407 for (
int axis = 0; axis <
numDims; ++axis)
1409 if (axis < commPads.size() ) _commPadSizes[ axis] = commPads[ axis];
1410 else _commPadSizes[ axis] = commPad;
1411 if (axis < bndryPads.size()) _bndryPadSizes[axis] = bndryPads[axis];
1412 else _bndryPadSizes[axis] = bndryPad;
1413 if (_mdComm->isPeriodic(axis))
1414 _bndryPad.push_back(Teuchos::tuple(_commPadSizes[axis],
1415 _commPadSizes[axis]));
1417 _bndryPad.push_back(Teuchos::tuple(_bndryPadSizes[axis],
1418 _bndryPadSizes[axis]));
1419 _globalDims[axis] = dimensions[axis] + _bndryPad[axis][0] +
1424 lower = _bndryPadSizes[axis];
1426 lower = _commPadSizes[axis];
1428 upper = _bndryPadSizes[axis];
1430 upper = _commPadSizes[axis];
1431 _pad.push_back(Teuchos::tuple(lower, upper));
1435 _globalMax = computeSize(_globalDims());
1439 _globalRankBounds.resize(
numDims);
1442 _localMax = computeSize(_localDims());
1445 Teuchos::Array< int > repBndry = plist.get(
"replicated boundary",
1446 Teuchos::Array< int >());
1447 _replicatedBoundary = createArrayOfInts(
numDims, repBndry);
1450 std::string layout = plist.get(
"layout",
"DEFAULT");
1451 std::transform(layout.begin(), layout.end(), layout.begin(), ::toupper);
1452 if (layout ==
"C ORDER")
1454 else if (layout ==
"FORTRAN ORDER")
1455 _layout = FORTRAN_ORDER;
1456 else if (layout ==
"ROW MAJOR")
1457 _layout = ROW_MAJOR;
1458 else if (layout ==
"COLUMN MAJOR")
1459 _layout = COLUMN_MAJOR;
1460 else if (layout ==
"LAST INDEX FASTEST")
1461 _layout = LAST_INDEX_FASTEST;
1462 else if (layout ==
"FIRST INDEX FASTEST")
1463 _layout = FIRST_INDEX_FASTEST;
1465 _layout = DEFAULT_ORDER;
1468 _globalStrides = computeStrides< size_type, dim_type >(_globalDims, _layout);
1469 _localStrides = computeStrides< size_type, dim_type >(_localDims , _layout);
1474 template<
class Node >
1476 MDMap(
const Teuchos::RCP< const MDComm > mdComm,
1477 const Teuchos::ArrayView< Slice > & myGlobalBounds,
1478 const Teuchos::ArrayView< padding_type > & padding,
1479 const Teuchos::ArrayView< const int > & replicatedBoundary,
1480 const Layout layout) :
1482 _globalDims(mdComm->numDims()),
1483 _globalBounds(mdComm->numDims()),
1484 _globalRankBounds(mdComm->numDims()),
1485 _globalStrides(mdComm->numDims()),
1488 _localDims(mdComm->numDims(), 0),
1489 _localBounds(mdComm->numDims()),
1490 _localStrides(mdComm->numDims()),
1493 _commPadSizes(mdComm->numDims(), 0),
1494 _pad(mdComm->numDims(), Teuchos::tuple(0,0)),
1495 _bndryPadSizes(mdComm->numDims(), 0),
1496 _bndryPad(mdComm->numDims()),
1497 _replicatedBoundary(createArrayOfInts(mdComm->numDims(),
1498 replicatedBoundary)),
1502 int numDims = _mdComm->numDims();
1503 TEUCHOS_TEST_FOR_EXCEPTION(
1504 myGlobalBounds.size() <
numDims,
1506 "MDMap: myGlobalBounds too small");
1509 int maxAxis = std::min(
numDims, (
int)padding.size());
1510 for (
int axis = 0; axis < maxAxis; ++axis)
1511 _pad[axis] = padding[axis];
1519 for (
int axis = 0; axis <
numDims; ++axis)
1520 _globalRankBounds[axis].resize(_mdComm->getCommDim(axis));
1523 int numProc = _mdComm->getTeuchosComm()->getSize();
1525 FIRST_INDEX_FASTEST);
1527 FIRST_INDEX_FASTEST);
1528 for (
int axis = 0; axis <
numDims; ++axis)
1530 sendBuffer(0,axis) = mdComm->getCommIndex(axis);
1531 sendBuffer(1,axis) = myGlobalBounds[axis].start();
1532 sendBuffer(2,axis) = myGlobalBounds[axis].stop();
1533 sendBuffer(3,axis) = _pad[axis][0];
1534 sendBuffer(4,axis) = _pad[axis][1];
1538 Teuchos::gatherAll(*(_mdComm->getTeuchosComm()),
1539 (
int)sendBuffer.
size(),
1541 (int)recvBuffer.
size(),
1547 for (
int axis = 0; axis <
numDims; ++axis)
1549 for (
int commIndex = 0; commIndex < _mdComm->getCommDim(axis); ++commIndex)
1553 for (
int rank = 0; rank < numProc; ++rank)
1555 if (recvBuffer(0,axis,rank) == commIndex)
1557 dim_type start = recvBuffer(1,axis,rank);
1558 dim_type stop = recvBuffer(2,axis,rank);
1559 int loPad = recvBuffer(3,axis,rank);
1560 int hiPad = recvBuffer(4,axis,rank);
1563 bounds =
Slice(start, stop);
1569 TEUCHOS_TEST_FOR_EXCEPTION(
1570 (bounds.
start() != start) || (bounds.
stop() != stop),
1572 "Global rank bounds mismatch: bounds = " << bounds <<
1573 ", (start,stop) = (" << start <<
"," << stop <<
")");
1574 TEUCHOS_TEST_FOR_EXCEPTION(
1575 (pad[0] != loPad) || (pad[1] != hiPad),
1577 "Padding value mismatch: pad = " << pad <<
", (loPad,hiPad) = (" 1578 << loPad <<
"," << hiPad <<
")");
1584 if (commIndex == 0 ) _bndryPad[axis][0] = pad[0];
1585 if (commIndex == mdComm->getCommDim(axis)-1) _bndryPad[axis][1] = pad[1];
1588 _globalRankBounds[axis][commIndex] = bounds;
1593 for (
int axis = 0; axis <
numDims; ++axis)
1595 for (
int commIndex = 1; commIndex < _mdComm->getCommDim(axis); ++commIndex)
1597 TEUCHOS_TEST_FOR_EXCEPTION(
1598 _globalRankBounds[axis][commIndex-1].stop() !=
1599 _globalRankBounds[axis][commIndex ].start(),
1601 "Global rank bounds are not contiguous");
1606 for (
int axis = 0; axis <
numDims; ++axis)
1608 int commSize = _mdComm->getCommDim(axis);
1610 _globalRankBounds[axis][0 ].start() - _pad[axis][0];
1612 _globalRankBounds[axis][commSize-1].stop() + _pad[axis][1];
1613 _globalDims[axis] = stop - start;
1614 _globalBounds[axis] =
Slice(start, stop);
1616 _globalStrides = computeStrides< size_type, dim_type >(_globalDims, _layout);
1619 for (
int axis = 0; axis <
numDims; ++axis)
1621 _globalMin += _globalBounds[axis].start() * _globalStrides[axis];
1622 _globalMax += _globalBounds[axis].stop() * _globalStrides[axis];
1626 for (
int axis = 0; axis <
numDims; ++axis)
1628 int commIndex = _mdComm->getCommIndex(axis);
1630 _globalRankBounds[axis][commIndex].start() - _pad[axis][0];
1632 _globalRankBounds[axis][commIndex].stop() + _pad[axis][1];
1633 _localDims[axis] = stop - start;
1634 _localBounds[axis] =
Slice(stop - start);
1636 _localStrides = computeStrides< size_type, dim_type >(_localDims, _layout);
1639 for (
int axis = 0; axis <
numDims; ++axis)
1640 _localMax += (_localDims[axis] - 1) * _localStrides[axis];
1645 template<
class Node >
1648 _mdComm(source._mdComm),
1649 _globalDims(source._globalDims),
1650 _globalBounds(source._globalBounds),
1651 _globalRankBounds(source._globalRankBounds),
1652 _globalStrides(source._globalStrides),
1653 _globalMin(source._globalMin),
1654 _globalMax(source._globalMax),
1655 _localDims(source._localDims),
1656 _localBounds(source._localBounds),
1657 _localStrides(source._localStrides),
1658 _localMin(source._localMin),
1659 _localMax(source._localMax),
1660 _commPadSizes(source._commPadSizes),
1662 _bndryPadSizes(source._bndryPadSizes),
1663 _bndryPad(source._bndryPad),
1664 _replicatedBoundary(source._replicatedBoundary),
1665 _layout(source._layout)
1671 template<
class Node >
1676 _mdComm(parent._mdComm),
1679 _globalRankBounds(),
1692 _replicatedBoundary(),
1693 _layout(parent._layout)
1698 TEUCHOS_TEST_FOR_EXCEPTION(
1699 ((axis < 0) || (axis >=
numDims)),
1701 "axis = " << axis <<
" is invalid for communicator with " <<
1705 TEUCHOS_TEST_FOR_EXCEPTION(
1706 ((index < globalBounds.
start()) || (index >= globalBounds.
stop())),
1708 "index = " << index <<
" is invalid for MDMap axis " <<
1709 axis <<
" with bounds " << globalBounds);
1713 int thisAxisRank = -1;
1714 for (
int axisRank = 0; axisRank < parent.
getCommDim(axis); ++axisRank)
1715 if (index >= parent._globalRankBounds[axis][axisRank].start() &&
1716 index < parent._globalRankBounds[axis][axisRank].stop())
1717 thisAxisRank = axisRank;
1718 TEUCHOS_TEST_FOR_EXCEPTION(
1719 (thisAxisRank == -1),
1721 "error computing axis rank for sub-communicator");
1722 _mdComm = Teuchos::rcp(
new MDComm(*(parent._mdComm), axis, thisAxisRank));
1729 if (_mdComm->onSubcommunicator())
1734 _globalDims.push_back(1);
1736 Teuchos::Array< Slice > bounds(1);
1738 _globalRankBounds.push_back(bounds);
1739 _globalStrides.push_back(1);
1740 _globalMin = index * parent._globalStrides[axis];
1741 _globalMax = _globalMin;
1742 _localDims.push_back(1);
1744 _localStrides.push_back(1);
1745 _localMin = parent._localMin +
1746 (index - parent._globalRankBounds[axis][0].start()) *
1747 parent._localStrides[axis];
1748 _localMax = _localMin + 1;
1749 _commPadSizes.push_back(0);
1750 _pad.push_back(Teuchos::tuple(0,0));
1751 _bndryPadSizes.push_back(0);
1752 _bndryPad.push_back(Teuchos::tuple(0,0));
1753 _replicatedBoundary.push_back(0);
1757 _globalMin = parent._globalMin;
1758 _globalMax = parent._globalMax;
1759 _localMin = parent._localMin;
1760 _localMax = parent._localMax;
1761 for (
int myAxis = 0; myAxis <
numDims; ++myAxis)
1765 _globalDims.push_back(parent._globalDims[myAxis]);
1766 _globalBounds.push_back(parent._globalBounds[myAxis]);
1767 _globalRankBounds.push_back(parent._globalRankBounds[myAxis]);
1768 _globalStrides.push_back(parent._globalStrides[myAxis]);
1769 _localDims.push_back(parent._localDims[myAxis]);
1770 _localBounds.push_back(parent._localBounds[myAxis]);
1771 _localStrides.push_back(parent._localStrides[myAxis]);
1772 _commPadSizes.push_back(parent._commPadSizes[myAxis]);
1773 _pad.push_back(parent._pad[myAxis]);
1774 _bndryPadSizes.push_back(parent._bndryPadSizes[myAxis]);
1775 _bndryPad.push_back(parent._bndryPad[myAxis]);
1776 _replicatedBoundary.push_back(parent._replicatedBoundary[myAxis]);
1781 _globalMin += index * parent._globalStrides[axis];
1782 _globalMax -= (parent._globalBounds[axis].stop() - index) *
1783 parent._globalStrides[axis];
1784 _localMin += (index-parent._globalRankBounds[axis][axisRank].start())
1785 * parent._localStrides[axis];
1786 _localMax -= (parent._globalRankBounds[axis][axisRank].stop()-index-1)
1787 * parent._localStrides[axis];
1796 template<
class Node >
1800 const Slice & slice,
1802 _mdComm(parent._mdComm),
1803 _globalDims(parent._globalDims),
1804 _globalBounds(parent._globalBounds),
1805 _globalRankBounds(parent._globalRankBounds),
1806 _globalStrides(parent._globalStrides),
1807 _globalMin(parent._globalMin),
1808 _globalMax(parent._globalMax),
1809 _localDims(parent._localDims),
1810 _localBounds(parent._localBounds),
1811 _localStrides(parent._localStrides),
1812 _localMin(parent._localMin),
1813 _localMax(parent._localMax),
1814 _commPadSizes(parent._commPadSizes),
1816 _bndryPadSizes(parent._bndryPadSizes),
1817 _bndryPad(parent._bndryPad),
1818 _replicatedBoundary(parent._replicatedBoundary),
1819 _layout(parent._layout)
1827 TEUCHOS_TEST_FOR_EXCEPTION(
1828 ((axis < 0) || (axis >=
numDims)),
1830 "axis = " << axis <<
" is invalid for MDMap with " <<
1836 TEUCHOS_TEST_FOR_EXCEPTION(
1840 "Slice along axis " << axis <<
" is " << bounds <<
" but must be within " 1842 TEUCHOS_TEST_FOR_EXCEPTION(
1845 "Slice along axis " << axis <<
" has length zero");
1849 _bndryPadSizes[axis] = bndryPad;
1850 _bndryPad[axis][0] = bndryPad;
1851 _bndryPad[axis][1] = bndryPad;
1854 for (
int axisRank = 0; axisRank < parent.
getCommDim(axis);
1857 dim_type start = _globalRankBounds[axis][axisRank].start();
1858 dim_type stop = _globalRankBounds[axis][axisRank].stop();
1859 if (start < bounds.
start()) start = bounds.
start();
1860 if (stop > bounds.
stop() ) stop = bounds.
stop();
1861 _globalRankBounds[axis][axisRank] =
ConcreteSlice(start, stop);
1865 dim_type start = bounds.
start() - _bndryPadSizes[axis];
1868 _bndryPad[axis][0] = bounds.
start();
1871 dim_type stop = bounds.
stop() + _bndryPadSizes[axis];
1881 _globalDims[axis] = stop - start;
1882 _globalMin += start * _globalStrides[axis];
1883 _globalMax -= (parent.
getGlobalDim(axis,
true) - stop) *
1884 _globalStrides[axis];
1887 if ((parent.
getGlobalBounds(axis,
true).start() < _globalBounds[axis].start())
1888 || (parent.
getGlobalBounds(axis,
true).stop() > _globalBounds[axis].stop()))
1889 _replicatedBoundary[axis] = 0;
1894 for (
int axisRank = 0; axisRank < parent.
getCommDim(axis);
1897 if ((_globalRankBounds[axis][axisRank].start() - _bndryPad[axis][0]
1898 <= _globalBounds[axis].start()) &&
1899 (_globalBounds[axis].start() <
1900 _globalRankBounds[axis][axisRank].stop() + _bndryPad[axis][1]))
1901 if (pStart == -1) pStart = axisRank;
1902 if ((_globalRankBounds[axis][axisRank].start() - _bndryPad[axis][0]
1903 < _globalBounds[axis].stop()) &&
1904 (_globalBounds[axis].stop() <=
1905 _globalRankBounds[axis][axisRank].stop() + _bndryPad[axis][1]))
1908 TEUCHOS_TEST_FOR_EXCEPTION(
1909 (pStart == -1 || pStop == -1),
1911 "error computing axis rank slice");
1915 _mdComm = Teuchos::rcp(
new MDComm(*(parent._mdComm), axis, axisRankSlice));
1920 if (_mdComm->onSubcommunicator())
1924 int myAxisRank = _mdComm->getCommIndex(axis);
1925 if (myAxisRank == 0)
1926 _pad[axis][0] = _bndryPad[axis][0];
1927 if (myAxisRank == _mdComm->getCommDim(axis)-1)
1928 _pad[axis][1] = _bndryPad[axis][1];
1934 dim_type start = (_globalRankBounds[axis][parentAxisRank].start() -
1936 (parent._globalRankBounds[axis][parentAxisRank].start() -
1937 parent._pad[axis][0]);
1938 dim_type stop = (_globalRankBounds[axis][parentAxisRank].stop() +
1940 (parent._globalRankBounds[axis][parentAxisRank].start() -
1941 parent._pad[axis][0]);
1945 _localDims[axis] = stop - start;
1946 _localMin += start * _localStrides[axis];
1948 _localStrides[axis];
1952 Teuchos::Array< Slice > newRankBounds;
1953 for (
int axisRank = 0; axisRank < parent.
getCommDim(axis);
1955 if ((axisRank >= axisRankSlice.
start()) &&
1956 (axisRank < axisRankSlice.
stop() ) )
1957 newRankBounds.push_back(_globalRankBounds[axis][axisRank]);
1958 _globalRankBounds[axis] = newRankBounds;
1965 _localBounds.clear();
1966 _commPadSizes.clear();
1968 _bndryPadSizes.clear();
1970 _localStrides.clear();
1977 template<
class Node >
1980 const Teuchos::ArrayView< Slice > & slices,
1981 const Teuchos::ArrayView< int > & bndryPad)
1984 int numDims = parent.
numDims();
1987 TEUCHOS_TEST_FOR_EXCEPTION(
1988 (slices.size() != numDims),
1990 "number of slices = " << slices.size() <<
" != parent MDMap number of " 1991 "dimensions = " << numDims);
1995 for (
int axis = 0; axis < numDims; ++axis)
1997 int bndryPadding = (axis < bndryPad.size()) ? bndryPad[axis] : 0;
2002 tempMDMap1 = tempMDMap2;
2009 template<
class Node >
2016 template<
class Node >
2020 _mdComm = source._mdComm;
2021 _globalDims = source._globalDims;
2022 _globalBounds = source._globalBounds;
2023 _globalRankBounds = source._globalRankBounds;
2024 _globalStrides = source._globalStrides;
2025 _globalMin = source._globalMin;
2026 _globalMax = source._globalMax;
2027 _localDims = source._localDims;
2028 _localBounds = source._localBounds;
2029 _localStrides = source._localStrides;
2030 _localMin = source._localMin;
2031 _localMax = source._localMax;
2032 _commPadSizes = source._commPadSizes;
2034 _bndryPadSizes = source._bndryPadSizes;
2035 _bndryPad = source._bndryPad;
2036 _layout = source._layout;
2042 template<
class Node >
2043 Teuchos::RCP< const MDComm >
2051 template<
class Node >
2055 return _mdComm->onSubcommunicator();
2060 template<
class Node >
2061 Teuchos::RCP< const Teuchos::Comm< int > >
2064 return _mdComm->getTeuchosComm();
2069 template<
class Node >
2073 return _mdComm->numDims();
2078 template<
class Node >
2079 Teuchos::Array< int >
2082 return _mdComm->getCommDims();
2087 template<
class Node >
2091 return _mdComm->getCommDim(axis);
2096 template<
class Node >
2100 return _mdComm->isPeriodic(axis);
2105 template<
class Node >
2109 return _mdComm->getCommIndex(axis);
2114 template<
class Node >
2118 return _mdComm->getLowerNeighbor(axis);
2123 template<
class Node >
2127 return _mdComm->getUpperNeighbor(axis);
2132 template<
class Node >
2133 Teuchos::Array< dim_type >
2142 template<
class Node >
2146 bool withBndryPad)
const 2148 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2149 TEUCHOS_TEST_FOR_EXCEPTION(
2150 ((axis < 0) || (axis >= numDims())),
2152 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2156 return _globalDims[axis];
2158 return _globalDims[axis] - _bndryPad[axis][0] - _bndryPad[axis][1];
2163 template<
class Node >
2167 bool withBndryPad)
const 2169 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2170 TEUCHOS_TEST_FOR_EXCEPTION(
2171 ((axis < 0) || (axis >= numDims())),
2173 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2177 return _globalBounds[axis];
2180 dim_type start = _globalBounds[axis].start() + _bndryPad[axis][0];
2181 dim_type stop = _globalBounds[axis].stop() - _bndryPad[axis][1];
2188 template<
class Node >
2194 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2195 TEUCHOS_TEST_FOR_EXCEPTION(
2196 ((axis < 0) || (axis >= numDims())),
2198 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2202 return _localDims[axis];
2204 return _localDims[axis] - _pad[axis][0] - _pad[axis][1];
2209 template<
class Node >
2210 Teuchos::Array< dim_type >
2219 template<
class Node >
2223 bool withBndryPad)
const 2225 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2226 TEUCHOS_TEST_FOR_EXCEPTION(
2227 ((axis < 0) || (axis >= numDims())),
2229 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2232 int axisRank = getCommIndex(axis);
2235 dim_type start = _globalRankBounds[axis][axisRank].start();
2236 dim_type stop = _globalRankBounds[axis][axisRank].stop();
2237 if (getCommIndex(axis) == 0)
2238 start -= _bndryPad[axis][0];
2239 if (getCommIndex(axis) == getCommDim(axis)-1)
2240 stop += _bndryPad[axis][1];
2244 return _globalRankBounds[axis][axisRank];
2249 template<
class Node >
2250 Teuchos::ArrayView< const Slice >
2254 return _localBounds();
2259 template<
class Node >
2265 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2266 TEUCHOS_TEST_FOR_EXCEPTION(
2267 ((axis < 0) || (axis >= numDims())),
2269 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2273 return _localBounds[axis];
2276 dim_type start = _localBounds[axis].start() + _pad[axis][0];
2277 dim_type stop = _localBounds[axis].stop() - _pad[axis][1];
2284 template<
class Node >
2289 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2290 TEUCHOS_TEST_FOR_EXCEPTION(
2291 ((axis < 0) || (axis >= numDims())),
2293 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2296 dim_type start = _localBounds[axis].start() + _pad[axis][0];
2297 dim_type stop = _localBounds[axis].stop() - _pad[axis][1];
2298 if (_mdComm->getLowerNeighbor(axis) == -1) ++start;
2299 if (_mdComm->getUpperNeighbor(axis) == -1) --stop;
2305 template<
class Node >
2309 bool result =
false;
2310 for (
int axis = 0; axis < numDims(); ++axis)
2311 if (_pad[axis][0] + _pad[axis][1]) result =
true;
2317 template<
class Node >
2321 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2322 TEUCHOS_TEST_FOR_EXCEPTION(
2323 ((axis < 0) || (axis >= numDims())),
2325 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2328 return _pad[axis][0];
2333 template<
class Node >
2337 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2338 TEUCHOS_TEST_FOR_EXCEPTION(
2339 ((axis < 0) || (axis >= numDims())),
2341 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2344 return _pad[axis][1];
2349 template<
class Node >
2353 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2354 TEUCHOS_TEST_FOR_EXCEPTION(
2355 ((axis < 0) || (axis >= numDims())),
2357 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2360 return _commPadSizes[axis];
2365 template<
class Node >
2369 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2370 TEUCHOS_TEST_FOR_EXCEPTION(
2371 ((axis < 0) || (axis >= numDims())),
2373 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2376 return _bndryPad[axis][0];
2381 template<
class Node >
2385 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2386 TEUCHOS_TEST_FOR_EXCEPTION(
2387 ((axis < 0) || (axis >= numDims())),
2389 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2392 return _bndryPad[axis][1];
2397 template<
class Node >
2398 Teuchos::Array< int >
2401 return _bndryPadSizes;
2406 template<
class Node >
2410 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2411 TEUCHOS_TEST_FOR_EXCEPTION(
2412 ((axis < 0) || (axis >= numDims())),
2414 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2417 return _bndryPadSizes[axis];
2422 template<
class Node >
2425 isPad(
const Teuchos::ArrayView< dim_type > & index)
const 2427 bool result =
false;
2428 for (
int axis = 0; axis < numDims(); ++axis)
2430 if (index[axis] < getLowerPadSize(axis))
2432 if (index[axis] >= getLocalDim(axis,
true) - getUpperPadSize(axis))
2440 template<
class Node >
2443 isCommPad(
const Teuchos::ArrayView< dim_type > & index)
const 2445 bool result =
false;
2446 for (
int axis = 0; axis < numDims(); ++axis)
2451 if (getLowerNeighbor(axis) >= 0)
2453 if (index[axis] < getLowerPadSize(axis))
2456 if (getUpperNeighbor(axis) >= 0)
2458 if (index[axis] >= getLocalDim(axis,
true) - getUpperPadSize(axis))
2467 template<
class Node >
2470 isBndryPad(
const Teuchos::ArrayView< dim_type > & index)
const 2472 bool result =
false;
2473 for (
int axis = 0; axis < numDims(); ++axis)
2478 if (getLowerNeighbor(axis) == -1)
2480 if (index[axis] < getLowerPadSize(axis))
2483 if (getUpperNeighbor(axis) == -1)
2485 if (index[axis] >= getLocalDim(axis,
true) - getUpperPadSize(axis))
2494 template<
class Node >
2498 return _mdComm->isPeriodic(axis) && bool(_replicatedBoundary[axis]);
2503 template<
class Node >
2512 template<
class Node >
2513 Teuchos::RCP< Node >
2521 template<
class Node >
2522 Teuchos::ArrayView< Teuchos::RCP< const MDMap< Node > > >
2525 if (_axisMaps.size() == 0) _axisMaps.resize(numDims());
2526 for (
int axis = 0; axis < numDims(); ++axis)
2527 if (_axisMaps[axis].is_null()) getAxisMap(axis);
2533 template<
class Node >
2534 Teuchos::RCP< const MDMap< Node > >
2537 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2538 TEUCHOS_TEST_FOR_EXCEPTION(
2539 ((axis < 0) || (axis >= numDims())),
2541 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2544 if (_axisMaps.size() == 0) _axisMaps.resize(numDims());
2545 if (_axisMaps[axis].is_null())
2547 Teuchos::RCP< const MDComm > axisComm = _mdComm->getAxisComm(axis);
2548 Domi::dim_type axisDim = _globalDims[axis] - 2*_bndryPadSizes[axis];
2550 Teuchos::rcp(
new MDMap(axisComm,
2551 Teuchos::tuple(axisDim),
2552 Teuchos::tuple(_commPadSizes[axis]),
2553 Teuchos::tuple(_bndryPadSizes[axis]),
2554 Teuchos::tuple(_replicatedBoundary[axis]),
2557 return _axisMaps[axis];
2562 template<
class Node >
2563 Teuchos::RCP< const MDMap< Node > >
2565 const dim_type trailingDim)
const 2572 if (leadingDim > 0) ++newNumDims;
2573 if (trailingDim > 0) ++newNumDims;
2576 if (newNumDims == 0)
return Teuchos::rcp(newMdMap);
2579 int oldNumDims = numDims();
2580 Teuchos::Array< int > newCommDims(oldNumDims);
2581 Teuchos::Array< int > newPeriodic(oldNumDims);
2582 for (
int axis = 0; axis < oldNumDims; ++axis)
2584 newCommDims[axis] = getCommDim(axis);
2585 newPeriodic[axis] = int(isPeriodic(axis));
2589 newCommDims.insert(newCommDims.begin(),1);
2590 newPeriodic.insert(newPeriodic.begin(),0);
2592 if (trailingDim > 0)
2594 newCommDims.push_back(1);
2595 newPeriodic.push_back(0);
2597 newMdMap->_mdComm = Teuchos::rcp(
new MDComm(getTeuchosComm(),
2603 padding_type pad(Teuchos::tuple(0,0));
2606 newMdMap->_globalDims.insert(newMdMap->_globalDims.begin(), leadingDim);
2607 newMdMap->_globalBounds.insert(newMdMap->_globalBounds.begin(), slice);
2608 newMdMap->_globalRankBounds.insert(newMdMap->_globalRankBounds.begin(),
2609 Teuchos::Array< Slice >(1,slice));
2610 newMdMap->_localDims.insert(newMdMap->_localDims.begin(), leadingDim);
2611 newMdMap->_localBounds.insert(newMdMap->_localBounds.begin(), slice);
2612 newMdMap->_commPadSizes.insert(newMdMap->_commPadSizes.begin(),0);
2613 newMdMap->_pad.insert(newMdMap->_pad.begin(), pad);
2614 newMdMap->_bndryPadSizes.insert(newMdMap->_bndryPadSizes.begin(),0);
2615 newMdMap->_bndryPad.insert(newMdMap->_bndryPad.begin(), pad);
2619 slice =
Slice(trailingDim);
2620 if (trailingDim > 0)
2622 newMdMap->_globalDims.push_back(trailingDim);
2623 newMdMap->_globalBounds.push_back(slice);
2624 newMdMap->_globalRankBounds.push_back(Teuchos::Array< Slice >(1,slice));
2625 newMdMap->_localDims.push_back(trailingDim);
2626 newMdMap->_localBounds.push_back(slice);
2627 newMdMap->_commPadSizes.push_back(0);
2628 newMdMap->_pad.push_back(pad);
2629 newMdMap->_bndryPadSizes.push_back(0);
2630 newMdMap->_bndryPad.push_back(pad);
2634 newMdMap->_globalStrides =
2635 computeStrides< size_type, dim_type >(newMdMap->_globalDims,
2637 newMdMap->_localStrides =
2638 computeStrides< size_type, dim_type >(newMdMap->_localDims,
2640 newMdMap->_globalMin = 0;
2641 newMdMap->_globalMax = 0;
2642 newMdMap->_localMin = 0;
2643 newMdMap->_localMax = 0;
2644 for (
int axis = 0; axis < oldNumDims + newNumDims; ++axis)
2646 newMdMap->_globalMin += newMdMap->_globalBounds[axis].start() *
2647 newMdMap->_globalStrides[axis];
2648 newMdMap->_globalMax += newMdMap->_globalBounds[axis].stop() *
2649 newMdMap->_globalStrides[axis];
2650 newMdMap->_localMin += newMdMap->_localBounds[axis].start() *
2651 newMdMap->_localStrides[axis];
2652 newMdMap->_localMax += newMdMap->_localBounds[axis].stop() *
2653 newMdMap->_localStrides[axis];
2657 return Teuchos::rcp(newMdMap);
2664 template<
class Node >
2665 Teuchos::RCP< const Epetra_Map >
2670 if (_epetraMap.is_null())
2674 TEUCHOS_TEST_FOR_EXCEPTION(
2675 computeSize(_globalDims) - 1 > std::numeric_limits< int >::max(),
2677 "The maximum global ID of this MDMap is too large for an Epetra_Map");
2680 int num_dims = numDims();
2681 Teuchos::Array<dim_type> localDims(num_dims);
2682 for (
int axis = 0; axis < num_dims; ++axis)
2683 localDims[axis] = _localDims[axis];
2685 Teuchos::Array<int> index(num_dims);
2689 it != myElements.end(); ++it)
2692 for (
int axis = 0; axis < num_dims; ++axis)
2694 int axisRank = getCommIndex(axis);
2695 int start = _globalRankBounds[axis][axisRank].start() -
2697 globalID += (start + it.index(axis)) * _globalStrides[axis];
2703 Teuchos::RCP< const Epetra_Comm > epetraComm = _mdComm->getEpetraComm();
2704 _epetraMap = Teuchos::rcp(
new Epetra_Map(-1,
2706 myElements.getRawPtr(),
2714 if (_epetraOwnMap.is_null())
2718 if (computeSize(_globalDims) - 1 > std::numeric_limits< int >::max())
2719 throw MapOrdinalError(
"The maximum global ID of this MDMap is too " 2720 "large for an Epetra_Map");
2723 int num_dims = numDims();
2724 Teuchos::Array<int> index(num_dims);
2725 Teuchos::Array<dim_type> myDims(num_dims);
2726 for (
int axis = 0; axis < num_dims; ++axis)
2728 myDims[axis] = _localDims[axis] - _pad[axis][0] - _pad[axis][1];
2729 int axisRank = getCommIndex(axis);
2731 myDims[axis] += _bndryPad[axis][0];
2732 if (axisRank == getCommDim(axis)-1)
2733 myDims[axis] += _bndryPad[axis][1];
2735 MDArray<int> myElements(myDims());
2738 for (MDArray<int>::iterator it = myElements.begin();
2739 it != myElements.end(); ++it)
2742 for (
int axis = 0; axis < num_dims; ++axis)
2744 int axisRank = getCommIndex(axis);
2745 int start = _globalRankBounds[axis][axisRank].start();
2747 start -= _bndryPad[axis][0];
2748 if (axisRank == getCommDim(axis)-1)
2749 start += _bndryPad[axis][1];
2750 globalID += (start + it.index(axis)) * _globalStrides[axis];
2755 Teuchos::RCP< const Epetra_Comm > epetraComm = _mdComm->getEpetraComm();
2756 _epetraOwnMap = Teuchos::rcp(
new Epetra_Map(-1,
2758 myElements.getRawPtr(),
2762 return _epetraOwnMap;
2772 template<
class Node >
2773 Teuchos::RCP< const Epetra_Map >
2775 getEpetraAxisMap(
int axis,
2776 bool withCommPad)
const 2778 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2779 TEUCHOS_TEST_FOR_EXCEPTION(
2780 ((axis < 0) || (axis >= numDims())),
2782 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2785 if ((withCommPad && (_epetraAxisMaps.size() == 0)) ||
2786 (not withCommPad && (_epetraAxisOwnMaps.size() == 0)))
2788 int num_dims = numDims();
2789 Teuchos::RCP< const Epetra_Comm > epetraComm = _mdComm->getEpetraComm();
2790 for (
int axis=0; axis < num_dims; ++axis)
2792 Teuchos::Array<int> elements(getLocalDim(axis, withCommPad));
2793 int start = getGlobalRankBounds(axis,
true).start();
2794 if (withCommPad && (getCommIndex(axis) != 0)) start -= _pad[axis][0];
2795 for (
int i = 0; i < elements.size(); ++i)
2796 elements[i] = i + start;
2799 _epetraAxisMaps.push_back(
2800 Teuchos::rcp(
new Epetra_Map(-1,
2802 elements.getRawPtr(),
2808 _epetraAxisOwnMaps.push_back(
2809 Teuchos::rcp(
new Epetra_Map(-1,
2811 elements.getRawPtr(),
2819 return _epetraAxisMaps[axis];
2821 return _epetraAxisOwnMaps[axis];
2830 template<
class Node >
2831 template<
class LocalOrdinal,
2832 class GlobalOrdinal,
2834 Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node2 > >
2835 MDMap< Node >::getTpetraMap(
bool withCommPad)
const 2840 int num_dims = numDims();
2841 Teuchos::Array<dim_type> localDims(num_dims);
2842 for (
int axis = 0; axis < num_dims; ++axis)
2843 localDims[axis] = _localDims[axis];
2844 MDArray< GlobalOrdinal > elementMDArray(localDims);
2845 Teuchos::Array< LocalOrdinal > index(num_dims);
2848 for (
typename MDArray< GlobalOrdinal >::iterator it = elementMDArray.begin();
2849 it != elementMDArray.end(); ++it)
2851 GlobalOrdinal globalID = 0;
2852 for (
int axis = 0; axis < num_dims; ++axis)
2854 int axisRank = getCommIndex(axis);
2855 GlobalOrdinal start = _globalRankBounds[axis][axisRank].start() -
2857 globalID += (start + it.index(axis)) * _globalStrides[axis];
2863 const Teuchos::Array< GlobalOrdinal > & myElements =
2864 elementMDArray.array();
2865 Teuchos::RCP< const Teuchos::Comm< int > > teuchosComm =
2866 _mdComm->getTeuchosComm();
2868 Teuchos::rcp(
new Tpetra::Map< LocalOrdinal,
2870 Node2 >(Teuchos::OrdinalTraits<
2871 Tpetra::global_size_t>::invalid(),
2879 int num_dims = numDims();
2880 Teuchos::Array< LocalOrdinal > index(num_dims);
2881 Teuchos::Array< dim_type > myDims(num_dims);
2882 for (
int axis = 0; axis < num_dims; ++axis)
2885 _localDims[axis] - _pad[axis][0] - _pad[axis][1];
2886 int axisRank = getCommIndex(axis);
2888 myDims[axis] += _bndryPad[axis][0];
2889 if (axisRank == getCommDim(axis)-1)
2890 myDims[axis] += _bndryPad[axis][1];
2892 MDArray< GlobalOrdinal > elementMDArray(myDims());
2895 for (
typename MDArray< GlobalOrdinal >::iterator it = elementMDArray.begin();
2896 it != elementMDArray.end(); ++it)
2898 GlobalOrdinal globalID = 0;
2899 for (
int axis = 0; axis < num_dims; ++axis)
2901 int axisRank = getCommIndex(axis);
2902 GlobalOrdinal start = _globalRankBounds[axis][axisRank].start();
2904 start -= _bndryPad[axis][0];
2905 if (axisRank == getCommDim(axis)-1)
2906 start += _bndryPad[axis][1];
2907 globalID += (start + it.index(axis)) * _globalStrides[axis];
2912 const Teuchos::Array< GlobalOrdinal> & myElements =
2913 elementMDArray.array();
2914 Teuchos::RCP< const Teuchos::Comm< int > > teuchosComm =
2915 _mdComm->getTeuchosComm();
2917 Teuchos::rcp(
new Tpetra::Map< LocalOrdinal,
2919 Node >(Teuchos::OrdinalTraits<
2920 Tpetra::global_size_t>::invalid(),
2932 template<
class Node >
2933 template<
class LocalOrdinal,
2934 class GlobalOrdinal,
2936 Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node2 > >
2938 getTpetraAxisMap(
int axis,
2939 bool withCommPad)
const 2941 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2942 TEUCHOS_TEST_FOR_EXCEPTION(
2943 ((axis < 0) || (axis >= numDims())),
2945 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2948 int num_dims = numDims();
2949 Teuchos::RCP< const Teuchos::Comm< int > > teuchosComm =
2950 _mdComm->getTeuchosComm();
2951 Teuchos::Array< GlobalOrdinal > elements(getLocalDim(axis,withCommPad));
2952 GlobalOrdinal start = getGlobalRankBounds(axis,
true).start();
2953 if (withCommPad && (getCommIndex(axis) != 0)) start -= _pad[axis][0];
2954 for (LocalOrdinal i = 0; i < elements.size(); ++i)
2955 elements[i] = i + start;
2956 return Teuchos::rcp(
new Tpetra::Map< LocalOrdinal,
2958 Node >(Teuchos::OrdinalTraits<
2959 Tpetra::global_size_t>::invalid(),
2968 template<
class Node >
2969 Teuchos::Array< dim_type >
2973 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2974 TEUCHOS_TEST_FOR_EXCEPTION(
2975 ((globalID < _globalMin) || (globalID >= _globalMax)),
2977 "invalid global index = " << globalID <<
" (should be between " <<
2978 _globalMin <<
" and " << _globalMax <<
")");
2980 int num_dims = numDims();
2981 Teuchos::Array< dim_type > result(num_dims);
2982 size_type index = globalID;
2983 if (_layout == LAST_INDEX_FASTEST)
2985 for (
int axis = 0; axis < num_dims-1; ++axis)
2987 result[axis] = index / _globalStrides[axis];
2988 index = index % _globalStrides[axis];
2990 result[num_dims-1] = index;
2994 for (
int axis = num_dims-1; axis > 0; --axis)
2996 result[axis] = index / _globalStrides[axis];
2997 index = index % _globalStrides[axis];
3006 template<
class Node >
3007 Teuchos::Array< dim_type >
3011 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 3012 TEUCHOS_TEST_FOR_EXCEPTION(
3013 ((localID < _localMin) || (localID >= _localMax)),
3015 "invalid local index = " << localID <<
" (should be between " <<
3016 _localMin <<
" and " << _localMax <<
")");
3018 int num_dims = numDims();
3019 Teuchos::Array< dim_type > result(num_dims);
3020 size_type index = localID;
3021 if (_layout == LAST_INDEX_FASTEST)
3023 for (
int axis = 0; axis < num_dims-1; ++axis)
3025 result[axis] = index / _localStrides[axis];
3026 index = index % _localStrides[axis];
3028 result[num_dims-1] = index;
3032 for (
int axis = num_dims-1; axis > 0; --axis)
3034 result[axis] = index / _localStrides[axis];
3035 index = index % _localStrides[axis];
3044 template<
class Node >
3049 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 3050 TEUCHOS_TEST_FOR_EXCEPTION(
3051 ((localID < 0) || (localID >= _localMax)),
3053 "invalid local index = " << localID <<
" (local size = " <<
3056 Teuchos::Array< dim_type > localIndex = getLocalIndex(localID);
3057 size_type result = 0;
3058 for (
int axis = 0; axis < numDims(); ++axis)
3060 dim_type globalIndex = localIndex[axis] +
3061 _globalRankBounds[axis][getCommIndex(axis)].start() - _pad[axis][0];
3062 result += globalIndex * _globalStrides[axis];
3069 template<
class Node >
3074 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 3075 TEUCHOS_TEST_FOR_EXCEPTION(
3076 (globalIndex.size() != numDims()),
3078 "globalIndex has " << globalIndex.size() <<
" entries; expecting " 3080 for (
int axis = 0; axis < numDims(); ++axis)
3082 TEUCHOS_TEST_FOR_EXCEPTION(
3083 ((globalIndex[axis] < 0) ||
3084 (globalIndex[axis] >= _globalDims[axis])),
3086 "invalid globalIndex[" << axis <<
"] = " << globalIndex[axis] <<
3087 " (global dimension = " << _globalDims[axis] <<
")");
3090 size_type result = 0;
3091 for (
int axis = 0; axis < numDims(); ++axis)
3092 result += globalIndex[axis] * _globalStrides[axis];
3098 template<
class Node >
3103 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 3104 TEUCHOS_TEST_FOR_EXCEPTION(
3105 ((globalID < _globalMin) || (globalID >= _globalMax)),
3107 "invalid global index = " << globalID <<
" (should be between " <<
3108 _globalMin <<
" and " << _globalMax <<
")");
3110 Teuchos::Array< dim_type > globalIndex =
3111 getGlobalIndex(globalID);
3112 size_type result = 0;
3113 for (
int axis = 0; axis < numDims(); ++axis)
3115 dim_type localIndex = globalIndex[axis] -
3116 _globalRankBounds[axis][getCommIndex(axis)].start() + _pad[axis][0];
3117 TEUCHOS_TEST_FOR_EXCEPTION(
3118 (localIndex < 0 || localIndex >= _localDims[axis]),
3120 "global index not on local processor")
3121 result += localIndex * _localStrides[axis];
3128 template<
class Node >
3131 getLocalID(
const Teuchos::ArrayView< dim_type > & localIndex)
const 3133 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 3134 TEUCHOS_TEST_FOR_EXCEPTION(
3135 (localIndex.size() != numDims()),
3137 "localIndex has " << localIndex.size() <<
" entries; expecting " 3139 for (
int axis = 0; axis < numDims(); ++axis)
3141 TEUCHOS_TEST_FOR_EXCEPTION(
3142 ((localIndex[axis] < 0) ||
3143 (localIndex[axis] >= _localDims[axis])),
3145 "invalid localIndex[" << axis <<
"] = " << localIndex[axis] <<
3146 " (local dimension = " << _localDims[axis] <<
")");
3149 size_type result = 0;
3150 for (
int axis = 0; axis < numDims(); ++axis)
3151 result += localIndex[axis] * _localStrides[axis];
3157 template<
class Node >
3163 if (
this == &mdMap)
return true;
3167 int num_dims = numDims();
3168 if (num_dims != mdMap.
numDims())
return false;
3172 for (
int axis = 0; axis < num_dims; ++axis)
3173 if (getCommDim(axis) != mdMap.
getCommDim(axis))
return false;
3177 if (_globalDims != mdMap._globalDims)
return false;
3182 int localResult = 1;
3183 int globalResult = 1;
3184 for (
int axis = 0; axis < num_dims; ++axis)
3185 if (getLocalDim(axis,
false) != mdMap.
getLocalDim(axis,
false))
3187 Teuchos::reduceAll(*(getTeuchosComm()),
3188 Teuchos::REDUCE_MIN,
3194 return bool(globalResult);
3199 template<
class Node >
3202 const int verbose)
const 3206 if (
this == &mdMap)
return true;
3212 int localResult = 1;
3213 Teuchos::RCP< const Teuchos::Comm< int > > comm = getTeuchosComm();
3214 int rank = comm->getRank();
3217 if (! isCompatible(mdMap))
3221 std::cout << rank <<
": MDMaps are incompatible" << std::endl;
3229 std::cout << rank <<
": this Comm size = " << comm->getSize() <<
" != " 3238 std::cout << rank <<
": this Comm rank = " << rank <<
" != " 3243 if (_globalBounds != mdMap._globalBounds)
3247 std::cout << rank <<
": global bounds " << _globalBounds <<
" != " 3248 << mdMap._globalBounds << std::endl;
3252 if (_localDims != mdMap._localDims)
3256 std::cout << rank <<
": local dimensions " << _localDims <<
" != " 3257 << mdMap._localDims << std::endl;
3261 int globalResult = 1;
3262 Teuchos::reduceAll(*(getTeuchosComm()),
3263 Teuchos::REDUCE_MIN,
3269 return bool(globalResult);
3274 template<
class Node >
3279 Teuchos::Array< size_type > contiguousStrides =
3280 computeStrides< size_type, dim_type >(_localDims, _layout);
3283 int localResult = int(_localStrides != contiguousStrides);
3286 int globalResult = 0;
3287 Teuchos::reduceAll(*(_mdComm->getTeuchosComm()),
3288 Teuchos::REDUCE_SUM,
3292 return (globalResult == 0);
3297 template<
class Node >
3302 int num_dims = numDims();
3305 for (
int axis = 0; axis < num_dims; ++axis)
3308 int commDim = getCommDim(axis);
3309 for (
int axisRank = 0; axisRank < commDim; ++axisRank)
3314 dim_type localDim = (_globalDims[axis] - _bndryPad[axis][0] -
3315 _bndryPad[axis][1]) / commDim;
3316 dim_type axisStart = axisRank * localDim;
3325 dim_type remainder = (_globalDims[axis] - _bndryPad[axis][0] -
3326 _bndryPad[axis][1]) % commDim;
3327 if (commDim - axisRank - 1 < remainder)
3330 axisStart += (remainder - commDim + axisRank);
3334 axisStart += _bndryPad[axis][0];
3337 _globalRankBounds[axis].push_back(
3338 ConcreteSlice(axisStart, axisStart + localDim));
3342 if (axisRank == getCommIndex(axis))
3347 _localDims[axis] = localDim + _pad[axis][0] + _pad[axis][1];
3350 _localBounds.push_back(ConcreteSlice(_localDims[axis]));
Teuchos::RCP< const MDComm > getMDComm() const
Access the underlying MDComm object.
Definition: Domi_MDMap.hpp:2044
MDMap(const Teuchos::RCP< const MDComm > mdComm, const Teuchos::ArrayView< const dim_type > &dimensions, const Teuchos::ArrayView< const int > &commPad=Teuchos::ArrayView< const int >(), const Teuchos::ArrayView< const int > &bndryPad=Teuchos::ArrayView< const int >(), const Teuchos::ArrayView< const int > &replicatedBoundary=Teuchos::ArrayView< const int >(), const Layout layout=DEFAULT_ORDER)
Main constructor.
Definition: Domi_MDMap.hpp:1037
Teuchos::Array< int > getCommDims() const
Get the communicator sizes along each axis.
Definition: Domi_MDMap.hpp:2080
size_type getGlobalID(size_type localID) const
Convert a local ID to a global ID.
Definition: Domi_MDMap.hpp:3047
Teuchos::RCP< Node > getNode() const
Get the Kokkos node.
Definition: Domi_MDMap.hpp:2514
Teuchos::RCP< const MDMap< Node > > getAugmentedMDMap(const dim_type leadingDim, const dim_type trailingDim=0) const
Return an RCP to a new MDMap that is a simple augmentation of this MDMap.
Definition: Domi_MDMap.hpp:2564
int getCommDim(int axis) const
Get the communicator size along the given axis.
Definition: Domi_MDMap.hpp:2089
MDMap< Node > & operator=(const MDMap< Node > &source)
Assignment operator.
Definition: Domi_MDMap.hpp:2018
int getLowerPadSize(int axis) const
Get the size of the lower padding along the given axis.
Definition: Domi_MDMap.hpp:2319
Iterator class suitable for multi-dimensional arrays.
Definition: Domi_MDIterator.hpp:100
Teuchos::Array< dim_type > getGlobalIndex(size_type globalID) const
Convert a global ID to a global index.
Definition: Domi_MDMap.hpp:2971
Range Error exception type.
Definition: Domi_Exceptions.hpp:65
Slice getLocalInteriorBounds(int axis) const
Get the local interior loop bounds along the specified axis.
Definition: Domi_MDMap.hpp:2287
bool isReplicatedBoundary(int axis) const
Return whether the given axis supports a replicated boundary.
Definition: Domi_MDMap.hpp:2496
int numDims() const
Get the number of dimensions.
Definition: Domi_MDMap.hpp:2071
bool onSubcommunicator() const
Query whether this processor is on the sub-communicator.
Definition: Domi_MDMap.hpp:2053
dim_type getLocalDim(int axis, bool withPad=false) const
Get the local dimension along the specified axis.
Definition: Domi_MDMap.hpp:2191
bool isCompatible(const MDMap< Node > &mdMap) const
True if two MDMaps are "compatible".
Definition: Domi_MDMap.hpp:3159
A Slice defines a subset of a container.
Multi-dimensional communicator object.
Definition: Domi_MDComm.hpp:108
Slice getGlobalRankBounds(int axis, bool withBndryPad=false) const
Get the global loop bounds along the specified axis.
Definition: Domi_MDMap.hpp:2222
Bounds Error exception type.
Definition: Domi_Exceptions.hpp:138
const dim_type stop() const
Stop index.
Definition: Domi_Slice.hpp:438
bool isPeriodic(int axis) const
Return the periodic flag for the given axis.
Definition: Domi_MDMap.hpp:2098
A ConcreteSlice is a Slice that does not accept Default or negative start or stop values...
Definition: Domi_Slice.hpp:340
A Slice contains a start, stop, and step index, describing a subset of an ordered container...
Definition: Domi_Slice.hpp:137
Teuchos::ArrayView< Teuchos::RCP< const Domi::MDMap< Node > > > getAxisMaps() const
Return an array of axis maps.
Definition: Domi_MDMap.hpp:2523
Teuchos::RCP< const Teuchos::Comm< int > > getTeuchosComm() const
Get the Teuchos communicator.
Definition: Domi_MDMap.hpp:2062
Teuchos::RCP< const Domi::MDMap< Node > > getAxisMap(int axis) const
Return an axis map for the given axis.
Definition: Domi_MDMap.hpp:2535
Teuchos::Array< int > getBndryPadSizes() const
Get the boundary padding sizes along each axis.
Definition: Domi_MDMap.hpp:2399
Teuchos::ArrayView< const Slice > getLocalBounds() const
Get the local loop bounds along every axis.
Definition: Domi_MDMap.hpp:2252
int getUpperNeighbor(int axis) const
Get the rank of the upper neighbor.
Definition: Domi_MDMap.hpp:2125
MDMap Error exception type.
Definition: Domi_Exceptions.hpp:114
Definition: Domi_ConfigDefs.hpp:97
int getCommPadSize(int axis) const
Get the communication padding size along the given axis.
Definition: Domi_MDMap.hpp:2351
bool isSameAs(const MDMap< Node > &mdMap, const int verbose=0) const
True if two MDMaps are "identical".
Definition: Domi_MDMap.hpp:3201
virtual Slice bounds(dim_type size) const
Return a Slice with concrete start and stop values.
int getUpperPadSize(int axis) const
Get the size of the upper padding along the given axis.
Definition: Domi_MDMap.hpp:2335
bool isContiguous() const
True if there are no stride gaps on any processor.
Definition: Domi_MDMap.hpp:3276
static const dim_type Default
Default value for Slice constructors.
Definition: Domi_Slice.hpp:155
int getLowerBndryPad(int axis) const
Get the size of the lower boundary padding along the given axis.
Definition: Domi_MDMap.hpp:2367
Slice getGlobalBounds(int axis, bool withBndryPad=false) const
Get the bounds of the global problem.
Definition: Domi_MDMap.hpp:2166
int getUpperBndryPad(int axis) const
Get the size of the upper boundary padding along the given axis.
Definition: Domi_MDMap.hpp:2383
~MDMap()
Destructor.
Definition: Domi_MDMap.hpp:2010
dim_type getGlobalDim(int axis, bool withBndryPad=false) const
Get the global dimension along the specified axis.
Definition: Domi_MDMap.hpp:2145
bool hasPadding() const
Return true if there is any padding stored locally.
Definition: Domi_MDMap.hpp:2307
int getBndryPadSize(int axis) const
Get the boundary padding size along the given axis.
Definition: Domi_MDMap.hpp:2408
size_type getLocalID(size_type globalID) const
Convert a global ID to a local ID.
Definition: Domi_MDMap.hpp:3101
size_type size() const
Return the total size of the MDArray
Definition: Domi_MDArray.hpp:1037
int getLowerNeighbor(int axis) const
Get the rank of the lower neighbor.
Definition: Domi_MDMap.hpp:2116
Multi-dimensional map.
Definition: Domi_MDMap.hpp:144
const T * getRawPtr() const
Return a const raw pointer to the beginning of the MDArray or NULL if unsized.
Definition: Domi_MDArray.hpp:1714
const dim_type start() const
Start index.
Definition: Domi_Slice.hpp:431
Memory-safe templated multi-dimensional array class.
Definition: Domi_MDArray.hpp:65
Teuchos::Array< dim_type > getGlobalDims() const
Get an array of the the global dimensions, including boundary padding.
Definition: Domi_MDMap.hpp:2135
Teuchos::Array< dim_type > getLocalIndex(size_type localID) const
Convert a local ID to a local index.
Definition: Domi_MDMap.hpp:3009
Teuchos::Array< dim_type > getLocalDims() const
Get an array of the local dimensions, including padding.
Definition: Domi_MDMap.hpp:2212
Invalid argument exception type.
Definition: Domi_Exceptions.hpp:53
int getCommIndex(int axis) const
Get the axis rank of this processor.
Definition: Domi_MDMap.hpp:2107
Map Ordinal Error exception type.
Definition: Domi_Exceptions.hpp:90
Layout getLayout() const
Get the storage order.
Definition: Domi_MDMap.hpp:2505