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);
266 MDMap(
const Teuchos::RCP< const MDComm > mdComm,
267 const Teuchos::ArrayView< Slice > & myGlobalBounds,
268 const Teuchos::ArrayView< padding_type > & padding =
269 Teuchos::ArrayView< padding_type >(),
270 const Teuchos::ArrayView< const int > & replicatedBoundary =
271 Teuchos::ArrayView< const int >(),
272 const Layout layout = DEFAULT_ORDER);
337 const Teuchos::ArrayView< Slice > & slices,
338 const Teuchos::ArrayView< int > & bndryPad =
339 Teuchos::ArrayView< int >());
366 inline Teuchos::RCP< const MDComm >
getMDComm()
const;
382 inline Teuchos::RCP< const Teuchos::Comm< int > >
getTeuchosComm()
const;
489 bool withBndryPad=
false)
const;
500 bool withBndryPad=
false)
const;
516 bool withBndryPad=
false)
const;
531 bool withPad=
false)
const;
556 bool withPad=
false)
const;
665 bool isPad(
const Teuchos::ArrayView< dim_type > & index)
const;
673 bool isCommPad(
const Teuchos::ArrayView< dim_type > & index)
const;
681 bool isBndryPad(
const Teuchos::ArrayView< dim_type > & index)
const;
693 Teuchos::RCP< Node >
getNode()
const;
704 Teuchos::ArrayView< Teuchos::RCP< const Domi::MDMap< Node > > >
711 Teuchos::RCP< const Domi::MDMap< Node > >
getAxisMap(
int axis)
const;
729 Teuchos::RCP< const MDMap< Node > >
731 const dim_type trailingDim=0)
const;
743 Teuchos::RCP< const Epetra_Map >
744 getEpetraMap(
bool withCommPad=
true)
const;
756 Teuchos::RCP< const Epetra_Map >
757 getEpetraAxisMap(
int axis,
758 bool withCommPad=
true)
const;
773 template<
class LocalOrdinal,
774 class GlobalOrdinal = LocalOrdinal,
776 Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node2 > >
777 getTpetraMap(
bool withCommPad=
true)
const;
791 template<
class LocalOrdinal,
792 class GlobalOrdinal = LocalOrdinal,
794 Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node2 > >
795 getTpetraAxisMap(
int axis,
796 bool withCommPad=
true)
const;
809 Teuchos::Array< dim_type >
816 Teuchos::Array< dim_type >
830 getGlobalID(
const Teuchos::ArrayView< dim_type > & globalIndex)
const;
839 size_type
getLocalID(size_type globalID)
const;
846 getLocalID(
const Teuchos::ArrayView< dim_type > & localIndex)
const;
891 const int verbose = 0)
const;
911 void computeBounds();
914 Teuchos::RCP< const MDComm > _mdComm;
918 Teuchos::Array< dim_type > _globalDims;
922 Teuchos::Array< Slice > _globalBounds;
931 Teuchos::Array< Teuchos::Array< Slice > > _globalRankBounds;
936 Teuchos::Array< size_type > _globalStrides;
940 size_type _globalMin;
944 size_type _globalMax;
948 Teuchos::Array< dim_type > _localDims;
956 Teuchos::Array< Slice > _localBounds;
959 Teuchos::Array< size_type > _localStrides;
971 Teuchos::Array< int > _commPadSizes;
976 Teuchos::Array< padding_type > _pad;
979 Teuchos::Array< int > _bndryPadSizes;
987 Teuchos::Array< padding_type > _bndryPad;
994 Teuchos::Array< int > _replicatedBoundary;
1003 mutable Teuchos::Array< Teuchos::RCP< const MDMap< Node > > > _axisMaps;
1006 Teuchos::RCP< Node > _node;
1013 mutable Teuchos::RCP< const Epetra_Map > _epetraMap;
1020 mutable Teuchos::RCP< const Epetra_Map > _epetraOwnMap;
1026 mutable Teuchos::Array< Teuchos::RCP< const Epetra_Map > > _epetraAxisMaps;
1032 mutable Teuchos::Array< Teuchos::RCP< const Epetra_Map > > _epetraAxisOwnMaps;
1041 template<
class Node >
1043 MDMap(
const Teuchos::RCP< const MDComm > mdComm,
1044 const Teuchos::ArrayView< const dim_type > & dimensions,
1045 const Teuchos::ArrayView< const int > & commPad,
1046 const Teuchos::ArrayView< const int > & bndryPad,
1047 const Teuchos::ArrayView< const int > & replicatedBoundary,
1048 const Layout layout) :
1050 _globalDims(mdComm->numDims()),
1052 _globalRankBounds(mdComm->numDims()),
1056 _localDims(mdComm->numDims(), 0),
1061 _commPadSizes(mdComm->numDims(), 0),
1063 _bndryPadSizes(mdComm->numDims(), 0),
1065 _replicatedBoundary(createArrayOfInts(mdComm->numDims(),
1066 replicatedBoundary)),
1070 int numDims = mdComm->numDims();
1073 TEUCHOS_TEST_FOR_EXCEPTION(
1076 "Size of dimensions does not match MDComm number of dimensions");
1080 for (
int axis = 0; axis <
numDims; ++axis)
1082 if (axis < commPad.size() ) _commPadSizes[ axis] = commPad[ axis];
1083 if (axis < bndryPad.size()) _bndryPadSizes[axis] = bndryPad[axis];
1084 if (_mdComm->isPeriodic(axis))
1085 _bndryPad.push_back(Teuchos::tuple(_commPadSizes[axis],
1086 _commPadSizes[axis]));
1088 _bndryPad.push_back(Teuchos::tuple(_bndryPadSizes[axis],
1089 _bndryPadSizes[axis]));
1090 _globalDims[axis] = dimensions[axis] + _bndryPad[axis][0] +
1095 lower = _bndryPadSizes[axis];
1097 lower = _commPadSizes[axis];
1099 upper = _bndryPadSizes[axis];
1101 upper = _commPadSizes[axis];
1102 _pad.push_back(Teuchos::tuple(lower, upper));
1106 _globalMax = computeSize(_globalDims());
1111 _localMax = computeSize(_localDims());
1114 _globalStrides = computeStrides< size_type, dim_type >(_globalDims, _layout);
1115 _localStrides = computeStrides< size_type, dim_type >(_localDims , _layout);
1120 template<
class Node >
1123 _mdComm(Teuchos::rcp(new
MDComm(plist))),
1126 _globalRankBounds(),
1139 _replicatedBoundary(),
1147 int numDims = _mdComm->numDims();
1150 Teuchos::Array< dim_type > dimensions =
1151 plist.get(
"dimensions", Teuchos::Array< dim_type >());
1152 TEUCHOS_TEST_FOR_EXCEPTION(
1155 "Size of dimensions does not match MDComm number of dimensions");
1159 int commPad = plist.get(
"communication pad size",
int(0));
1160 int bndryPad = plist.get(
"boundary pad size" ,
int(0));
1161 Teuchos::Array< int > commPads =
1162 plist.get(
"communication pad sizes", Teuchos::Array< int >());
1163 Teuchos::Array< int > bndryPads =
1164 plist.get(
"boundary pad sizes" , Teuchos::Array< int >());
1165 _commPadSizes.resize(
numDims);
1166 _bndryPadSizes.resize(
numDims);
1171 for (
int axis = 0; axis <
numDims; ++axis)
1173 if (axis < commPads.size() ) _commPadSizes[ axis] = commPads[ axis];
1174 else _commPadSizes[ axis] = commPad;
1175 if (axis < bndryPads.size()) _bndryPadSizes[axis] = bndryPads[axis];
1176 else _bndryPadSizes[axis] = bndryPad;
1177 if (_mdComm->isPeriodic(axis))
1178 _bndryPad.push_back(Teuchos::tuple(_commPadSizes[axis],
1179 _commPadSizes[axis]));
1181 _bndryPad.push_back(Teuchos::tuple(_bndryPadSizes[axis],
1182 _bndryPadSizes[axis]));
1183 _globalDims[axis] = dimensions[axis] + _bndryPad[axis][0] +
1188 lower = _bndryPadSizes[axis];
1190 lower = _commPadSizes[axis];
1192 upper = _bndryPadSizes[axis];
1194 upper = _commPadSizes[axis];
1195 _pad.push_back(Teuchos::tuple(lower, upper));
1199 _globalMax = computeSize(_globalDims());
1203 _globalRankBounds.resize(
numDims);
1206 _localMax = computeSize(_localDims());
1209 Teuchos::Array< int > repBndry = plist.get(
"replicated boundary",
1210 Teuchos::Array< int >());
1211 _replicatedBoundary = createArrayOfInts(
numDims, repBndry);
1214 std::string layout = plist.get(
"layout",
"DEFAULT");
1215 std::transform(layout.begin(), layout.end(), layout.begin(), ::toupper);
1216 if (layout ==
"C ORDER")
1218 else if (layout ==
"FORTRAN ORDER")
1219 _layout = FORTRAN_ORDER;
1220 else if (layout ==
"ROW MAJOR")
1221 _layout = ROW_MAJOR;
1222 else if (layout ==
"COLUMN MAJOR")
1223 _layout = COLUMN_MAJOR;
1224 else if (layout ==
"LAST INDEX FASTEST")
1225 _layout = LAST_INDEX_FASTEST;
1226 else if (layout ==
"FIRST INDEX FASTEST")
1227 _layout = FIRST_INDEX_FASTEST;
1229 _layout = DEFAULT_ORDER;
1232 _globalStrides = computeStrides< size_type, dim_type >(_globalDims, _layout);
1233 _localStrides = computeStrides< size_type, dim_type >(_localDims , _layout);
1238 template<
class Node >
1240 MDMap(Teuchos::RCP<
const Teuchos::Comm< int > > teuchosComm,
1241 Teuchos::ParameterList & plist) :
1242 _mdComm(Teuchos::rcp(new
MDComm(teuchosComm, plist))),
1245 _globalRankBounds(),
1258 _replicatedBoundary(),
1266 int numDims = _mdComm->numDims();
1269 Teuchos::Array< dim_type > dimensions =
1270 plist.get(
"dimensions", Teuchos::Array< dim_type >());
1271 TEUCHOS_TEST_FOR_EXCEPTION(
1274 "Size of dimensions does not match MDComm number of dimensions");
1278 int commPad = plist.get(
"communication pad size",
int(0));
1279 int bndryPad = plist.get(
"boundary pad size" ,
int(0));
1280 Teuchos::Array< int > commPads =
1281 plist.get(
"communication pad sizes", Teuchos::Array< int >());
1282 Teuchos::Array< int > bndryPads =
1283 plist.get(
"boundary pad sizes" , Teuchos::Array< int >());
1284 _commPadSizes.resize(
numDims);
1285 _bndryPadSizes.resize(
numDims);
1290 for (
int axis = 0; axis <
numDims; ++axis)
1292 if (axis < commPads.size() ) _commPadSizes[ axis] = commPads[ axis];
1293 else _commPadSizes[ axis] = commPad;
1294 if (axis < bndryPads.size()) _bndryPadSizes[axis] = bndryPads[axis];
1295 else _bndryPadSizes[axis] = bndryPad;
1296 if (_mdComm->isPeriodic(axis))
1297 _bndryPad.push_back(Teuchos::tuple(_commPadSizes[axis],
1298 _commPadSizes[axis]));
1300 _bndryPad.push_back(Teuchos::tuple(_bndryPadSizes[axis],
1301 _bndryPadSizes[axis]));
1302 _globalDims[axis] = dimensions[axis] + _bndryPad[axis][0] +
1307 lower = _bndryPadSizes[axis];
1309 lower = _commPadSizes[axis];
1311 upper = _bndryPadSizes[axis];
1313 upper = _commPadSizes[axis];
1314 _pad.push_back(Teuchos::tuple(lower, upper));
1323 _globalMax = computeSize(_globalDims());
1327 _globalRankBounds.resize(
numDims);
1330 _localMax = computeSize(_localDims());
1333 Teuchos::Array< int > repBndry = plist.get(
"replicated boundary",
1334 Teuchos::Array< int >());
1335 _replicatedBoundary = createArrayOfInts(
numDims, repBndry);
1338 std::string layout = plist.get(
"layout",
"DEFAULT");
1339 std::transform(layout.begin(), layout.end(), layout.begin(), ::toupper);
1340 if (layout ==
"C ORDER")
1342 else if (layout ==
"FORTRAN ORDER")
1343 _layout = FORTRAN_ORDER;
1344 else if (layout ==
"ROW MAJOR")
1345 _layout = ROW_MAJOR;
1346 else if (layout ==
"COLUMN MAJOR")
1347 _layout = COLUMN_MAJOR;
1348 else if (layout ==
"LAST INDEX FASTEST")
1349 _layout = LAST_INDEX_FASTEST;
1350 else if (layout ==
"FIRST INDEX FASTEST")
1351 _layout = FIRST_INDEX_FASTEST;
1353 _layout = DEFAULT_ORDER;
1356 _globalStrides = computeStrides< size_type, dim_type >(_globalDims, _layout);
1357 _localStrides = computeStrides< size_type, dim_type >(_localDims , _layout);
1362 template<
class Node >
1364 MDMap(Teuchos::RCP< const MDComm > mdComm,
1365 Teuchos::ParameterList & plist) :
1367 _globalDims(mdComm->numDims()),
1369 _globalRankBounds(mdComm->numDims()),
1373 _localDims(mdComm->numDims(), 0),
1378 _commPadSizes(mdComm->numDims(), 0),
1380 _bndryPadSizes(mdComm->numDims(), 0),
1382 _replicatedBoundary(),
1386 plist.validateParameters(*getValidParameters());
1389 int numDims = _mdComm->numDims();
1392 Teuchos::Array< dim_type > dimensions =
1393 plist.get(
"dimensions", Teuchos::Array< dim_type >());
1394 TEUCHOS_TEST_FOR_EXCEPTION(
1397 "Number of dimensions does not match MDComm number of dimensions");
1401 int commPad = plist.get(
"communication pad size",
int(0));
1402 int bndryPad = plist.get(
"boundary pad size" ,
int(0));
1403 Teuchos::Array< int > commPads =
1404 plist.get(
"communication pad sizes", Teuchos::Array< int >());
1405 Teuchos::Array< int > bndryPads =
1406 plist.get(
"boundary pad sizes" , Teuchos::Array< int >());
1407 _commPadSizes.resize(
numDims);
1408 _bndryPadSizes.resize(
numDims);
1413 for (
int axis = 0; axis <
numDims; ++axis)
1415 if (axis < commPads.size() ) _commPadSizes[ axis] = commPads[ axis];
1416 else _commPadSizes[ axis] = commPad;
1417 if (axis < bndryPads.size()) _bndryPadSizes[axis] = bndryPads[axis];
1418 else _bndryPadSizes[axis] = bndryPad;
1419 if (_mdComm->isPeriodic(axis))
1420 _bndryPad.push_back(Teuchos::tuple(_commPadSizes[axis],
1421 _commPadSizes[axis]));
1423 _bndryPad.push_back(Teuchos::tuple(_bndryPadSizes[axis],
1424 _bndryPadSizes[axis]));
1425 _globalDims[axis] = dimensions[axis] + _bndryPad[axis][0] +
1430 lower = _bndryPadSizes[axis];
1432 lower = _commPadSizes[axis];
1434 upper = _bndryPadSizes[axis];
1436 upper = _commPadSizes[axis];
1437 _pad.push_back(Teuchos::tuple(lower, upper));
1441 _globalMax = computeSize(_globalDims());
1445 _globalRankBounds.resize(
numDims);
1448 _localMax = computeSize(_localDims());
1451 Teuchos::Array< int > repBndry = plist.get(
"replicated boundary",
1452 Teuchos::Array< int >());
1453 _replicatedBoundary = createArrayOfInts(
numDims, repBndry);
1456 std::string layout = plist.get(
"layout",
"DEFAULT");
1457 std::transform(layout.begin(), layout.end(), layout.begin(), ::toupper);
1458 if (layout ==
"C ORDER")
1460 else if (layout ==
"FORTRAN ORDER")
1461 _layout = FORTRAN_ORDER;
1462 else if (layout ==
"ROW MAJOR")
1463 _layout = ROW_MAJOR;
1464 else if (layout ==
"COLUMN MAJOR")
1465 _layout = COLUMN_MAJOR;
1466 else if (layout ==
"LAST INDEX FASTEST")
1467 _layout = LAST_INDEX_FASTEST;
1468 else if (layout ==
"FIRST INDEX FASTEST")
1469 _layout = FIRST_INDEX_FASTEST;
1471 _layout = DEFAULT_ORDER;
1474 _globalStrides = computeStrides< size_type, dim_type >(_globalDims, _layout);
1475 _localStrides = computeStrides< size_type, dim_type >(_localDims , _layout);
1480 template<
class Node >
1482 MDMap(
const Teuchos::RCP< const MDComm > mdComm,
1483 const Teuchos::ArrayView< Slice > & myGlobalBounds,
1484 const Teuchos::ArrayView< padding_type > & padding,
1485 const Teuchos::ArrayView< const int > & replicatedBoundary,
1486 const Layout layout) :
1488 _globalDims(mdComm->numDims()),
1489 _globalBounds(mdComm->numDims()),
1490 _globalRankBounds(mdComm->numDims()),
1491 _globalStrides(mdComm->numDims()),
1494 _localDims(mdComm->numDims(), 0),
1495 _localBounds(mdComm->numDims()),
1496 _localStrides(mdComm->numDims()),
1499 _commPadSizes(mdComm->numDims(), 0),
1500 _pad(mdComm->numDims(), Teuchos::tuple(0,0)),
1501 _bndryPadSizes(mdComm->numDims(), 0),
1502 _bndryPad(mdComm->numDims()),
1503 _replicatedBoundary(createArrayOfInts(mdComm->numDims(),
1504 replicatedBoundary)),
1508 int numDims = _mdComm->numDims();
1509 TEUCHOS_TEST_FOR_EXCEPTION(
1510 myGlobalBounds.size() <
numDims,
1512 "MDMap: myGlobalBounds too small");
1515 int maxAxis = std::min(
numDims, (
int)padding.size());
1516 for (
int axis = 0; axis < maxAxis; ++axis)
1517 _pad[axis] = padding[axis];
1525 for (
int axis = 0; axis <
numDims; ++axis)
1526 _globalRankBounds[axis].resize(_mdComm->getCommDim(axis));
1529 int numProc = _mdComm->getTeuchosComm()->getSize();
1531 FIRST_INDEX_FASTEST);
1533 FIRST_INDEX_FASTEST);
1534 for (
int axis = 0; axis <
numDims; ++axis)
1536 sendBuffer(0,axis) = mdComm->getCommIndex(axis);
1537 sendBuffer(1,axis) = myGlobalBounds[axis].start();
1538 sendBuffer(2,axis) = myGlobalBounds[axis].stop();
1539 sendBuffer(3,axis) = _pad[axis][0];
1540 sendBuffer(4,axis) = _pad[axis][1];
1544 Teuchos::gatherAll(*(_mdComm->getTeuchosComm()),
1545 (
int)sendBuffer.
size(),
1547 (int)recvBuffer.
size(),
1553 for (
int axis = 0; axis <
numDims; ++axis)
1555 for (
int commIndex = 0; commIndex < _mdComm->getCommDim(axis); ++commIndex)
1559 for (
int rank = 0; rank < numProc; ++rank)
1561 if (recvBuffer(0,axis,rank) == commIndex)
1563 dim_type start = recvBuffer(1,axis,rank);
1564 dim_type stop = recvBuffer(2,axis,rank);
1565 int loPad = recvBuffer(3,axis,rank);
1566 int hiPad = recvBuffer(4,axis,rank);
1569 bounds =
Slice(start, stop);
1575 TEUCHOS_TEST_FOR_EXCEPTION(
1576 (bounds.
start() != start) || (bounds.
stop() != stop),
1578 "Global rank bounds mismatch: bounds = " << bounds <<
1579 ", (start,stop) = (" << start <<
"," << stop <<
")");
1580 TEUCHOS_TEST_FOR_EXCEPTION(
1581 (pad[0] != loPad) || (pad[1] != hiPad),
1583 "Padding value mismatch: pad = " << pad <<
", (loPad,hiPad) = (" 1584 << loPad <<
"," << hiPad <<
")");
1590 if (commIndex == 0 ) _bndryPad[axis][0] = pad[0];
1591 if (commIndex == mdComm->getCommDim(axis)-1) _bndryPad[axis][1] = pad[1];
1594 _globalRankBounds[axis][commIndex] = bounds;
1599 for (
int axis = 0; axis <
numDims; ++axis)
1601 for (
int commIndex = 1; commIndex < _mdComm->getCommDim(axis); ++commIndex)
1603 TEUCHOS_TEST_FOR_EXCEPTION(
1604 _globalRankBounds[axis][commIndex-1].stop() !=
1605 _globalRankBounds[axis][commIndex ].start(),
1607 "Global rank bounds are not contiguous");
1612 for (
int axis = 0; axis <
numDims; ++axis)
1614 int commSize = _mdComm->getCommDim(axis);
1616 _globalRankBounds[axis][0 ].start() - _pad[axis][0];
1618 _globalRankBounds[axis][commSize-1].stop() + _pad[axis][1];
1619 _globalDims[axis] = stop - start;
1620 _globalBounds[axis] =
Slice(start, stop);
1622 _globalStrides = computeStrides< size_type, dim_type >(_globalDims, _layout);
1625 for (
int axis = 0; axis <
numDims; ++axis)
1627 _globalMin += _globalBounds[axis].start() * _globalStrides[axis];
1628 _globalMax += _globalBounds[axis].stop() * _globalStrides[axis];
1632 for (
int axis = 0; axis <
numDims; ++axis)
1634 int commIndex = _mdComm->getCommIndex(axis);
1636 _globalRankBounds[axis][commIndex].start() - _pad[axis][0];
1638 _globalRankBounds[axis][commIndex].stop() + _pad[axis][1];
1639 _localDims[axis] = stop - start;
1640 _localBounds[axis] =
Slice(stop - start);
1642 _localStrides = computeStrides< size_type, dim_type >(_localDims, _layout);
1645 for (
int axis = 0; axis <
numDims; ++axis)
1646 _localMax += (_localDims[axis] - 1) * _localStrides[axis];
1651 template<
class Node >
1654 _mdComm(source._mdComm),
1655 _globalDims(source._globalDims),
1656 _globalBounds(source._globalBounds),
1657 _globalRankBounds(source._globalRankBounds),
1658 _globalStrides(source._globalStrides),
1659 _globalMin(source._globalMin),
1660 _globalMax(source._globalMax),
1661 _localDims(source._localDims),
1662 _localBounds(source._localBounds),
1663 _localStrides(source._localStrides),
1664 _localMin(source._localMin),
1665 _localMax(source._localMax),
1666 _commPadSizes(source._commPadSizes),
1668 _bndryPadSizes(source._bndryPadSizes),
1669 _bndryPad(source._bndryPad),
1670 _replicatedBoundary(source._replicatedBoundary),
1671 _layout(source._layout)
1677 template<
class Node >
1682 _mdComm(parent._mdComm),
1685 _globalRankBounds(),
1698 _replicatedBoundary(),
1699 _layout(parent._layout)
1704 TEUCHOS_TEST_FOR_EXCEPTION(
1705 ((axis < 0) || (axis >=
numDims)),
1707 "axis = " << axis <<
" is invalid for communicator with " <<
1711 TEUCHOS_TEST_FOR_EXCEPTION(
1712 ((index < globalBounds.
start()) || (index >= globalBounds.
stop())),
1714 "index = " << index <<
" is invalid for MDMap axis " <<
1715 axis <<
" with bounds " << globalBounds);
1719 int thisAxisRank = -1;
1720 for (
int axisRank = 0; axisRank < parent.
getCommDim(axis); ++axisRank)
1721 if (index >= parent._globalRankBounds[axis][axisRank].start() &&
1722 index < parent._globalRankBounds[axis][axisRank].stop())
1723 thisAxisRank = axisRank;
1724 TEUCHOS_TEST_FOR_EXCEPTION(
1725 (thisAxisRank == -1),
1727 "error computing axis rank for sub-communicator");
1728 _mdComm = Teuchos::rcp(
new MDComm(*(parent._mdComm), axis, thisAxisRank));
1735 if (_mdComm->onSubcommunicator())
1740 _globalDims.push_back(1);
1742 Teuchos::Array< Slice > bounds(1);
1744 _globalRankBounds.push_back(bounds);
1745 _globalStrides.push_back(1);
1746 _globalMin = index * parent._globalStrides[axis];
1747 _globalMax = _globalMin;
1748 _localDims.push_back(1);
1750 _localStrides.push_back(1);
1751 _localMin = parent._localMin +
1752 (index - parent._globalRankBounds[axis][0].start()) *
1753 parent._localStrides[axis];
1754 _localMax = _localMin + 1;
1755 _commPadSizes.push_back(0);
1756 _pad.push_back(Teuchos::tuple(0,0));
1757 _bndryPadSizes.push_back(0);
1758 _bndryPad.push_back(Teuchos::tuple(0,0));
1759 _replicatedBoundary.push_back(0);
1763 _globalMin = parent._globalMin;
1764 _globalMax = parent._globalMax;
1765 _localMin = parent._localMin;
1766 _localMax = parent._localMax;
1767 for (
int myAxis = 0; myAxis <
numDims; ++myAxis)
1771 _globalDims.push_back(parent._globalDims[myAxis]);
1772 _globalBounds.push_back(parent._globalBounds[myAxis]);
1773 _globalRankBounds.push_back(parent._globalRankBounds[myAxis]);
1774 _globalStrides.push_back(parent._globalStrides[myAxis]);
1775 _localDims.push_back(parent._localDims[myAxis]);
1776 _localBounds.push_back(parent._localBounds[myAxis]);
1777 _localStrides.push_back(parent._localStrides[myAxis]);
1778 _commPadSizes.push_back(parent._commPadSizes[myAxis]);
1779 _pad.push_back(parent._pad[myAxis]);
1780 _bndryPadSizes.push_back(parent._bndryPadSizes[myAxis]);
1781 _bndryPad.push_back(parent._bndryPad[myAxis]);
1782 _replicatedBoundary.push_back(parent._replicatedBoundary[myAxis]);
1787 _globalMin += index * parent._globalStrides[axis];
1788 _globalMax -= (parent._globalBounds[axis].stop() - index) *
1789 parent._globalStrides[axis];
1790 _localMin += (index-parent._globalRankBounds[axis][axisRank].start())
1791 * parent._localStrides[axis];
1792 _localMax -= (parent._globalRankBounds[axis][axisRank].stop()-index-1)
1793 * parent._localStrides[axis];
1802 template<
class Node >
1806 const Slice & slice,
1808 _mdComm(parent._mdComm),
1809 _globalDims(parent._globalDims),
1810 _globalBounds(parent._globalBounds),
1811 _globalRankBounds(parent._globalRankBounds),
1812 _globalStrides(parent._globalStrides),
1813 _globalMin(parent._globalMin),
1814 _globalMax(parent._globalMax),
1815 _localDims(parent._localDims),
1816 _localBounds(parent._localBounds),
1817 _localStrides(parent._localStrides),
1818 _localMin(parent._localMin),
1819 _localMax(parent._localMax),
1820 _commPadSizes(parent._commPadSizes),
1822 _bndryPadSizes(parent._bndryPadSizes),
1823 _bndryPad(parent._bndryPad),
1824 _replicatedBoundary(parent._replicatedBoundary),
1825 _layout(parent._layout)
1833 TEUCHOS_TEST_FOR_EXCEPTION(
1834 ((axis < 0) || (axis >=
numDims)),
1836 "axis = " << axis <<
" is invalid for MDMap with " <<
1842 TEUCHOS_TEST_FOR_EXCEPTION(
1846 "Slice along axis " << axis <<
" is " << bounds <<
" but must be within " 1848 TEUCHOS_TEST_FOR_EXCEPTION(
1851 "Slice along axis " << axis <<
" has length zero");
1855 _bndryPadSizes[axis] = bndryPad;
1856 _bndryPad[axis][0] = bndryPad;
1857 _bndryPad[axis][1] = bndryPad;
1860 for (
int axisRank = 0; axisRank < parent.
getCommDim(axis);
1863 dim_type start = _globalRankBounds[axis][axisRank].start();
1864 dim_type stop = _globalRankBounds[axis][axisRank].stop();
1865 if (start < bounds.
start()) start = bounds.
start();
1866 if (stop > bounds.
stop() ) stop = bounds.
stop();
1867 _globalRankBounds[axis][axisRank] =
ConcreteSlice(start, stop);
1871 dim_type start = bounds.
start() - _bndryPadSizes[axis];
1874 _bndryPad[axis][0] = bounds.
start();
1877 dim_type stop = bounds.
stop() + _bndryPadSizes[axis];
1887 _globalDims[axis] = stop - start;
1888 _globalMin += start * _globalStrides[axis];
1889 _globalMax -= (parent.
getGlobalDim(axis,
true) - stop) *
1890 _globalStrides[axis];
1893 if ((parent.
getGlobalBounds(axis,
true).start() < _globalBounds[axis].start())
1894 || (parent.
getGlobalBounds(axis,
true).stop() > _globalBounds[axis].stop()))
1895 _replicatedBoundary[axis] = 0;
1900 for (
int axisRank = 0; axisRank < parent.
getCommDim(axis);
1903 if ((_globalRankBounds[axis][axisRank].start() - _bndryPad[axis][0]
1904 <= _globalBounds[axis].start()) &&
1905 (_globalBounds[axis].start() <
1906 _globalRankBounds[axis][axisRank].stop() + _bndryPad[axis][1]))
1907 if (pStart == -1) pStart = axisRank;
1908 if ((_globalRankBounds[axis][axisRank].start() - _bndryPad[axis][0]
1909 < _globalBounds[axis].stop()) &&
1910 (_globalBounds[axis].stop() <=
1911 _globalRankBounds[axis][axisRank].stop() + _bndryPad[axis][1]))
1914 TEUCHOS_TEST_FOR_EXCEPTION(
1915 (pStart == -1 || pStop == -1),
1917 "error computing axis rank slice");
1921 _mdComm = Teuchos::rcp(
new MDComm(*(parent._mdComm), axis, axisRankSlice));
1926 if (_mdComm->onSubcommunicator())
1930 int myAxisRank = _mdComm->getCommIndex(axis);
1931 if (myAxisRank == 0)
1932 _pad[axis][0] = _bndryPad[axis][0];
1933 if (myAxisRank == _mdComm->getCommDim(axis)-1)
1934 _pad[axis][1] = _bndryPad[axis][1];
1940 dim_type start = (_globalRankBounds[axis][parentAxisRank].start() -
1942 (parent._globalRankBounds[axis][parentAxisRank].start() -
1943 parent._pad[axis][0]);
1944 dim_type stop = (_globalRankBounds[axis][parentAxisRank].stop() +
1946 (parent._globalRankBounds[axis][parentAxisRank].start() -
1947 parent._pad[axis][0]);
1951 _localDims[axis] = stop - start;
1952 _localMin += start * _localStrides[axis];
1954 _localStrides[axis];
1958 Teuchos::Array< Slice > newRankBounds;
1959 for (
int axisRank = 0; axisRank < parent.
getCommDim(axis);
1961 if ((axisRank >= axisRankSlice.
start()) &&
1962 (axisRank < axisRankSlice.
stop() ) )
1963 newRankBounds.push_back(_globalRankBounds[axis][axisRank]);
1964 _globalRankBounds[axis] = newRankBounds;
1971 _localBounds.clear();
1972 _commPadSizes.clear();
1974 _bndryPadSizes.clear();
1976 _localStrides.clear();
1983 template<
class Node >
1986 const Teuchos::ArrayView< Slice > & slices,
1987 const Teuchos::ArrayView< int > & bndryPad)
1990 int numDims = parent.
numDims();
1993 TEUCHOS_TEST_FOR_EXCEPTION(
1994 (slices.size() != numDims),
1996 "number of slices = " << slices.size() <<
" != parent MDMap number of " 1997 "dimensions = " << numDims);
2001 for (
int axis = 0; axis < numDims; ++axis)
2003 int bndryPadding = (axis < bndryPad.size()) ? bndryPad[axis] : 0;
2008 tempMDMap1 = tempMDMap2;
2015 template<
class Node >
2022 template<
class Node >
2026 _mdComm = source._mdComm;
2027 _globalDims = source._globalDims;
2028 _globalBounds = source._globalBounds;
2029 _globalRankBounds = source._globalRankBounds;
2030 _globalStrides = source._globalStrides;
2031 _globalMin = source._globalMin;
2032 _globalMax = source._globalMax;
2033 _localDims = source._localDims;
2034 _localBounds = source._localBounds;
2035 _localStrides = source._localStrides;
2036 _localMin = source._localMin;
2037 _localMax = source._localMax;
2038 _commPadSizes = source._commPadSizes;
2040 _bndryPadSizes = source._bndryPadSizes;
2041 _bndryPad = source._bndryPad;
2042 _layout = source._layout;
2048 template<
class Node >
2049 Teuchos::RCP< const MDComm >
2057 template<
class Node >
2061 return _mdComm->onSubcommunicator();
2066 template<
class Node >
2067 Teuchos::RCP< const Teuchos::Comm< int > >
2070 return _mdComm->getTeuchosComm();
2075 template<
class Node >
2079 return _mdComm->numDims();
2084 template<
class Node >
2085 Teuchos::Array< int >
2088 return _mdComm->getCommDims();
2093 template<
class Node >
2097 return _mdComm->getCommDim(axis);
2102 template<
class Node >
2106 return _mdComm->isPeriodic(axis);
2111 template<
class Node >
2115 return _mdComm->getCommIndex(axis);
2120 template<
class Node >
2124 return _mdComm->getLowerNeighbor(axis);
2129 template<
class Node >
2133 return _mdComm->getUpperNeighbor(axis);
2138 template<
class Node >
2139 Teuchos::Array< dim_type >
2148 template<
class Node >
2152 bool withBndryPad)
const 2154 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2155 TEUCHOS_TEST_FOR_EXCEPTION(
2156 ((axis < 0) || (axis >= numDims())),
2158 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2162 return _globalDims[axis];
2164 return _globalDims[axis] - _bndryPad[axis][0] - _bndryPad[axis][1];
2169 template<
class Node >
2173 bool withBndryPad)
const 2175 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2176 TEUCHOS_TEST_FOR_EXCEPTION(
2177 ((axis < 0) || (axis >= numDims())),
2179 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2183 return _globalBounds[axis];
2186 dim_type start = _globalBounds[axis].start() + _bndryPad[axis][0];
2187 dim_type stop = _globalBounds[axis].stop() - _bndryPad[axis][1];
2194 template<
class Node >
2200 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2201 TEUCHOS_TEST_FOR_EXCEPTION(
2202 ((axis < 0) || (axis >= numDims())),
2204 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2208 return _localDims[axis];
2210 return _localDims[axis] - _pad[axis][0] - _pad[axis][1];
2215 template<
class Node >
2216 Teuchos::Array< dim_type >
2225 template<
class Node >
2229 bool withBndryPad)
const 2231 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2232 TEUCHOS_TEST_FOR_EXCEPTION(
2233 ((axis < 0) || (axis >= numDims())),
2235 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2238 int axisRank = getCommIndex(axis);
2241 dim_type start = _globalRankBounds[axis][axisRank].start();
2242 dim_type stop = _globalRankBounds[axis][axisRank].stop();
2243 if (getCommIndex(axis) == 0)
2244 start -= _bndryPad[axis][0];
2245 if (getCommIndex(axis) == getCommDim(axis)-1)
2246 stop += _bndryPad[axis][1];
2250 return _globalRankBounds[axis][axisRank];
2255 template<
class Node >
2256 Teuchos::ArrayView< const Slice >
2260 return _localBounds();
2265 template<
class Node >
2271 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2272 TEUCHOS_TEST_FOR_EXCEPTION(
2273 ((axis < 0) || (axis >= numDims())),
2275 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2279 return _localBounds[axis];
2282 dim_type start = _localBounds[axis].start() + _pad[axis][0];
2283 dim_type stop = _localBounds[axis].stop() - _pad[axis][1];
2290 template<
class Node >
2295 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2296 TEUCHOS_TEST_FOR_EXCEPTION(
2297 ((axis < 0) || (axis >= numDims())),
2299 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2302 dim_type start = _localBounds[axis].start() + _pad[axis][0];
2303 dim_type stop = _localBounds[axis].stop() - _pad[axis][1];
2304 if (_mdComm->getLowerNeighbor(axis) == -1) ++start;
2305 if (_mdComm->getUpperNeighbor(axis) == -1) --stop;
2311 template<
class Node >
2315 bool result =
false;
2316 for (
int axis = 0; axis < numDims(); ++axis)
2317 if (_pad[axis][0] + _pad[axis][1]) result =
true;
2323 template<
class Node >
2327 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2328 TEUCHOS_TEST_FOR_EXCEPTION(
2329 ((axis < 0) || (axis >= numDims())),
2331 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2334 return _pad[axis][0];
2339 template<
class Node >
2343 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2344 TEUCHOS_TEST_FOR_EXCEPTION(
2345 ((axis < 0) || (axis >= numDims())),
2347 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2350 return _pad[axis][1];
2355 template<
class Node >
2359 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2360 TEUCHOS_TEST_FOR_EXCEPTION(
2361 ((axis < 0) || (axis >= numDims())),
2363 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2366 return _commPadSizes[axis];
2371 template<
class Node >
2375 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2376 TEUCHOS_TEST_FOR_EXCEPTION(
2377 ((axis < 0) || (axis >= numDims())),
2379 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2382 return _bndryPad[axis][0];
2387 template<
class Node >
2391 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2392 TEUCHOS_TEST_FOR_EXCEPTION(
2393 ((axis < 0) || (axis >= numDims())),
2395 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2398 return _bndryPad[axis][1];
2403 template<
class Node >
2404 Teuchos::Array< int >
2407 return _bndryPadSizes;
2412 template<
class Node >
2416 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2417 TEUCHOS_TEST_FOR_EXCEPTION(
2418 ((axis < 0) || (axis >= numDims())),
2420 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2423 return _bndryPadSizes[axis];
2428 template<
class Node >
2431 isPad(
const Teuchos::ArrayView< dim_type > & index)
const 2433 bool result =
false;
2434 for (
int axis = 0; axis < numDims(); ++axis)
2436 if (index[axis] < getLowerPadSize(axis))
2438 if (index[axis] >= getLocalDim(axis,
true) - getUpperPadSize(axis))
2446 template<
class Node >
2449 isCommPad(
const Teuchos::ArrayView< dim_type > & index)
const 2451 bool result =
false;
2452 for (
int axis = 0; axis < numDims(); ++axis)
2457 if (getLowerNeighbor(axis) >= 0)
2459 if (index[axis] < getLowerPadSize(axis))
2462 if (getUpperNeighbor(axis) >= 0)
2464 if (index[axis] >= getLocalDim(axis,
true) - getUpperPadSize(axis))
2473 template<
class Node >
2476 isBndryPad(
const Teuchos::ArrayView< dim_type > & index)
const 2478 bool result =
false;
2479 for (
int axis = 0; axis < numDims(); ++axis)
2484 if (getLowerNeighbor(axis) == -1)
2486 if (index[axis] < getLowerPadSize(axis))
2489 if (getUpperNeighbor(axis) == -1)
2491 if (index[axis] >= getLocalDim(axis,
true) - getUpperPadSize(axis))
2500 template<
class Node >
2504 return _mdComm->isPeriodic(axis) && bool(_replicatedBoundary[axis]);
2509 template<
class Node >
2518 template<
class Node >
2519 Teuchos::RCP< Node >
2527 template<
class Node >
2528 Teuchos::ArrayView< Teuchos::RCP< const MDMap< Node > > >
2531 if (_axisMaps.size() == 0) _axisMaps.resize(numDims());
2532 for (
int axis = 0; axis < numDims(); ++axis)
2533 if (_axisMaps[axis].is_null()) getAxisMap(axis);
2539 template<
class Node >
2540 Teuchos::RCP< const MDMap< Node > >
2543 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2544 TEUCHOS_TEST_FOR_EXCEPTION(
2545 ((axis < 0) || (axis >= numDims())),
2547 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2550 if (_axisMaps.size() == 0) _axisMaps.resize(numDims());
2551 if (_axisMaps[axis].is_null())
2553 Teuchos::RCP< const MDComm > axisComm = _mdComm->getAxisComm(axis);
2554 Domi::dim_type axisDim = _globalDims[axis] - 2*_bndryPadSizes[axis];
2556 Teuchos::rcp(
new MDMap(axisComm,
2557 Teuchos::tuple(axisDim),
2558 Teuchos::tuple(_commPadSizes[axis]),
2559 Teuchos::tuple(_bndryPadSizes[axis]),
2560 Teuchos::tuple(_replicatedBoundary[axis]),
2563 return _axisMaps[axis];
2568 template<
class Node >
2569 Teuchos::RCP< const MDMap< Node > >
2571 const dim_type trailingDim)
const 2578 if (leadingDim > 0) ++newNumDims;
2579 if (trailingDim > 0) ++newNumDims;
2582 if (newNumDims == 0)
return Teuchos::rcp(newMdMap);
2585 int oldNumDims = numDims();
2586 Teuchos::Array< int > newCommDims(oldNumDims);
2587 Teuchos::Array< int > newPeriodic(oldNumDims);
2588 Teuchos::Array< int > newReplicatedBndry(oldNumDims);
2590 for (
int axis = 0; axis < oldNumDims; ++axis)
2592 newCommDims[axis] = getCommDim(axis);
2593 newPeriodic[axis] = int(isPeriodic(axis));
2594 newReplicatedBndry[axis] = int(isReplicatedBoundary(axis));
2598 newCommDims.insert(newCommDims.begin(),1);
2599 newPeriodic.insert(newPeriodic.begin(),0);
2600 newReplicatedBndry.insert(newReplicatedBndry.begin(),0);
2602 if (trailingDim > 0)
2604 newCommDims.push_back(1);
2605 newPeriodic.push_back(0);
2606 newReplicatedBndry.push_back(0);
2608 newMdMap->_mdComm = Teuchos::rcp(
new MDComm(getTeuchosComm(),
2611 newMdMap->_replicatedBoundary = newReplicatedBndry;
2615 padding_type pad(Teuchos::tuple(0,0));
2618 newMdMap->_globalDims.insert(newMdMap->_globalDims.begin(), leadingDim);
2619 newMdMap->_globalBounds.insert(newMdMap->_globalBounds.begin(), slice);
2620 newMdMap->_globalRankBounds.insert(newMdMap->_globalRankBounds.begin(),
2621 Teuchos::Array< Slice >(1,slice));
2622 newMdMap->_localDims.insert(newMdMap->_localDims.begin(), leadingDim);
2623 newMdMap->_localBounds.insert(newMdMap->_localBounds.begin(), slice);
2624 newMdMap->_commPadSizes.insert(newMdMap->_commPadSizes.begin(),0);
2625 newMdMap->_pad.insert(newMdMap->_pad.begin(), pad);
2626 newMdMap->_bndryPadSizes.insert(newMdMap->_bndryPadSizes.begin(),0);
2627 newMdMap->_bndryPad.insert(newMdMap->_bndryPad.begin(), pad);
2631 slice =
Slice(trailingDim);
2632 if (trailingDim > 0)
2634 newMdMap->_globalDims.push_back(trailingDim);
2635 newMdMap->_globalBounds.push_back(slice);
2636 newMdMap->_globalRankBounds.push_back(Teuchos::Array< Slice >(1,slice));
2637 newMdMap->_localDims.push_back(trailingDim);
2638 newMdMap->_localBounds.push_back(slice);
2639 newMdMap->_commPadSizes.push_back(0);
2640 newMdMap->_pad.push_back(pad);
2641 newMdMap->_bndryPadSizes.push_back(0);
2642 newMdMap->_bndryPad.push_back(pad);
2646 newMdMap->_globalStrides =
2647 computeStrides< size_type, dim_type >(newMdMap->_globalDims,
2649 newMdMap->_localStrides =
2650 computeStrides< size_type, dim_type >(newMdMap->_localDims,
2652 newMdMap->_globalMin = 0;
2653 newMdMap->_globalMax = 0;
2654 newMdMap->_localMin = 0;
2655 newMdMap->_localMax = 0;
2656 for (
int axis = 0; axis < oldNumDims + newNumDims; ++axis)
2658 newMdMap->_globalMin += newMdMap->_globalBounds[axis].start() *
2659 newMdMap->_globalStrides[axis];
2660 newMdMap->_globalMax += newMdMap->_globalBounds[axis].stop() *
2661 newMdMap->_globalStrides[axis];
2662 newMdMap->_localMin += newMdMap->_localBounds[axis].start() *
2663 newMdMap->_localStrides[axis];
2664 newMdMap->_localMax += newMdMap->_localBounds[axis].stop() *
2665 newMdMap->_localStrides[axis];
2669 return Teuchos::rcp(newMdMap);
2676 template<
class Node >
2677 Teuchos::RCP< const Epetra_Map >
2682 if (_epetraMap.is_null())
2686 TEUCHOS_TEST_FOR_EXCEPTION(
2687 computeSize(_globalDims) - 1 > std::numeric_limits< int >::max(),
2689 "The maximum global ID of this MDMap is too large for an Epetra_Map");
2692 int num_dims = numDims();
2693 Teuchos::Array<dim_type> localDims(num_dims);
2694 for (
int axis = 0; axis < num_dims; ++axis)
2695 localDims[axis] = _localDims[axis];
2697 Teuchos::Array<int> index(num_dims);
2701 it != myElements.end(); ++it)
2704 for (
int axis = 0; axis < num_dims; ++axis)
2706 int axisRank = getCommIndex(axis);
2707 int start = _globalRankBounds[axis][axisRank].start() -
2709 globalID += (start + it.index(axis)) * _globalStrides[axis];
2715 Teuchos::RCP< const Epetra_Comm > epetraComm = _mdComm->getEpetraComm();
2716 _epetraMap = Teuchos::rcp(
new Epetra_Map(-1,
2718 myElements.getRawPtr(),
2726 if (_epetraOwnMap.is_null())
2730 if (computeSize(_globalDims) - 1 > std::numeric_limits< int >::max())
2731 throw MapOrdinalError(
"The maximum global ID of this MDMap is too " 2732 "large for an Epetra_Map");
2735 int num_dims = numDims();
2736 Teuchos::Array<int> index(num_dims);
2737 Teuchos::Array<dim_type> myDims(num_dims);
2738 for (
int axis = 0; axis < num_dims; ++axis)
2740 myDims[axis] = _localDims[axis] - _pad[axis][0] - _pad[axis][1];
2741 int axisRank = getCommIndex(axis);
2743 myDims[axis] += _bndryPad[axis][0];
2744 if (axisRank == getCommDim(axis)-1)
2745 myDims[axis] += _bndryPad[axis][1];
2747 MDArray<int> myElements(myDims());
2750 for (MDArray<int>::iterator it = myElements.begin();
2751 it != myElements.end(); ++it)
2754 for (
int axis = 0; axis < num_dims; ++axis)
2756 int axisRank = getCommIndex(axis);
2757 int start = _globalRankBounds[axis][axisRank].start();
2759 start -= _bndryPad[axis][0];
2760 if (axisRank == getCommDim(axis)-1)
2761 start += _bndryPad[axis][1];
2762 globalID += (start + it.index(axis)) * _globalStrides[axis];
2767 Teuchos::RCP< const Epetra_Comm > epetraComm = _mdComm->getEpetraComm();
2768 _epetraOwnMap = Teuchos::rcp(
new Epetra_Map(-1,
2770 myElements.getRawPtr(),
2774 return _epetraOwnMap;
2784 template<
class Node >
2785 Teuchos::RCP< const Epetra_Map >
2787 getEpetraAxisMap(
int axis,
2788 bool withCommPad)
const 2790 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2791 TEUCHOS_TEST_FOR_EXCEPTION(
2792 ((axis < 0) || (axis >= numDims())),
2794 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2797 if ((withCommPad && (_epetraAxisMaps.size() == 0)) ||
2798 (not withCommPad && (_epetraAxisOwnMaps.size() == 0)))
2800 int num_dims = numDims();
2801 Teuchos::RCP< const Epetra_Comm > epetraComm = _mdComm->getEpetraComm();
2802 for (
int axis=0; axis < num_dims; ++axis)
2804 Teuchos::Array<int> elements(getLocalDim(axis, withCommPad));
2805 int start = getGlobalRankBounds(axis,
true).start();
2806 if (withCommPad && (getCommIndex(axis) != 0)) start -= _pad[axis][0];
2807 for (
int i = 0; i < elements.size(); ++i)
2808 elements[i] = i + start;
2811 _epetraAxisMaps.push_back(
2812 Teuchos::rcp(
new Epetra_Map(-1,
2814 elements.getRawPtr(),
2820 _epetraAxisOwnMaps.push_back(
2821 Teuchos::rcp(
new Epetra_Map(-1,
2823 elements.getRawPtr(),
2831 return _epetraAxisMaps[axis];
2833 return _epetraAxisOwnMaps[axis];
2842 template<
class Node >
2843 template<
class LocalOrdinal,
2844 class GlobalOrdinal,
2846 Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node2 > >
2847 MDMap< Node >::getTpetraMap(
bool withCommPad)
const 2852 int num_dims = numDims();
2853 Teuchos::Array<dim_type> localDims(num_dims);
2854 for (
int axis = 0; axis < num_dims; ++axis)
2855 localDims[axis] = _localDims[axis];
2856 MDArray< GlobalOrdinal > elementMDArray(localDims);
2857 Teuchos::Array< LocalOrdinal > index(num_dims);
2860 for (
typename MDArray< GlobalOrdinal >::iterator it = elementMDArray.begin();
2861 it != elementMDArray.end(); ++it)
2863 GlobalOrdinal globalID = 0;
2864 for (
int axis = 0; axis < num_dims; ++axis)
2866 int axisRank = getCommIndex(axis);
2867 GlobalOrdinal start = _globalRankBounds[axis][axisRank].start() -
2869 globalID += (start + it.index(axis)) * _globalStrides[axis];
2875 const Teuchos::Array< GlobalOrdinal > & myElements =
2876 elementMDArray.array();
2877 Teuchos::RCP< const Teuchos::Comm< int > > teuchosComm =
2878 _mdComm->getTeuchosComm();
2880 Teuchos::rcp(
new Tpetra::Map< LocalOrdinal,
2882 Node2 >(Teuchos::OrdinalTraits<
2883 Tpetra::global_size_t>::invalid(),
2891 int num_dims = numDims();
2892 Teuchos::Array< LocalOrdinal > index(num_dims);
2893 Teuchos::Array< dim_type > myDims(num_dims);
2894 for (
int axis = 0; axis < num_dims; ++axis)
2897 _localDims[axis] - _pad[axis][0] - _pad[axis][1];
2898 int axisRank = getCommIndex(axis);
2900 myDims[axis] += _bndryPad[axis][0];
2901 if (axisRank == getCommDim(axis)-1)
2902 myDims[axis] += _bndryPad[axis][1];
2904 MDArray< GlobalOrdinal > elementMDArray(myDims());
2907 for (
typename MDArray< GlobalOrdinal >::iterator it = elementMDArray.begin();
2908 it != elementMDArray.end(); ++it)
2910 GlobalOrdinal globalID = 0;
2911 for (
int axis = 0; axis < num_dims; ++axis)
2913 int axisRank = getCommIndex(axis);
2914 GlobalOrdinal start = _globalRankBounds[axis][axisRank].start();
2916 start -= _bndryPad[axis][0];
2917 if (axisRank == getCommDim(axis)-1)
2918 start += _bndryPad[axis][1];
2919 globalID += (start + it.index(axis)) * _globalStrides[axis];
2924 const Teuchos::Array< GlobalOrdinal> & myElements =
2925 elementMDArray.array();
2926 Teuchos::RCP< const Teuchos::Comm< int > > teuchosComm =
2927 _mdComm->getTeuchosComm();
2929 Teuchos::rcp(
new Tpetra::Map< LocalOrdinal,
2931 Node >(Teuchos::OrdinalTraits<
2932 Tpetra::global_size_t>::invalid(),
2944 template<
class Node >
2945 template<
class LocalOrdinal,
2946 class GlobalOrdinal,
2948 Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node2 > >
2950 getTpetraAxisMap(
int axis,
2951 bool withCommPad)
const 2953 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2954 TEUCHOS_TEST_FOR_EXCEPTION(
2955 ((axis < 0) || (axis >= numDims())),
2957 "invalid axis index = " << axis <<
" (number of dimensions = " <<
2960 int num_dims = numDims();
2961 Teuchos::RCP< const Teuchos::Comm< int > > teuchosComm =
2962 _mdComm->getTeuchosComm();
2963 Teuchos::Array< GlobalOrdinal > elements(getLocalDim(axis,withCommPad));
2964 GlobalOrdinal start = getGlobalRankBounds(axis,
true).start();
2965 if (withCommPad && (getCommIndex(axis) != 0)) start -= _pad[axis][0];
2966 for (LocalOrdinal i = 0; i < elements.size(); ++i)
2967 elements[i] = i + start;
2968 return Teuchos::rcp(
new Tpetra::Map< LocalOrdinal,
2970 Node >(Teuchos::OrdinalTraits<
2971 Tpetra::global_size_t>::invalid(),
2980 template<
class Node >
2981 Teuchos::Array< dim_type >
2985 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 2986 TEUCHOS_TEST_FOR_EXCEPTION(
2987 ((globalID < _globalMin) || (globalID >= _globalMax)),
2989 "invalid global index = " << globalID <<
" (should be between " <<
2990 _globalMin <<
" and " << _globalMax <<
")");
2992 int num_dims = numDims();
2993 Teuchos::Array< dim_type > result(num_dims);
2994 size_type index = globalID;
2995 if (_layout == LAST_INDEX_FASTEST)
2997 for (
int axis = 0; axis < num_dims-1; ++axis)
2999 result[axis] = index / _globalStrides[axis];
3000 index = index % _globalStrides[axis];
3002 result[num_dims-1] = index;
3006 for (
int axis = num_dims-1; axis > 0; --axis)
3008 result[axis] = index / _globalStrides[axis];
3009 index = index % _globalStrides[axis];
3018 template<
class Node >
3019 Teuchos::Array< dim_type >
3023 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 3024 TEUCHOS_TEST_FOR_EXCEPTION(
3025 ((localID < _localMin) || (localID >= _localMax)),
3027 "invalid local index = " << localID <<
" (should be between " <<
3028 _localMin <<
" and " << _localMax <<
")");
3030 int num_dims = numDims();
3031 Teuchos::Array< dim_type > result(num_dims);
3032 size_type index = localID;
3033 if (_layout == LAST_INDEX_FASTEST)
3035 for (
int axis = 0; axis < num_dims-1; ++axis)
3037 result[axis] = index / _localStrides[axis];
3038 index = index % _localStrides[axis];
3040 result[num_dims-1] = index;
3044 for (
int axis = num_dims-1; axis > 0; --axis)
3046 result[axis] = index / _localStrides[axis];
3047 index = index % _localStrides[axis];
3056 template<
class Node >
3061 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 3062 TEUCHOS_TEST_FOR_EXCEPTION(
3063 ((localID < 0) || (localID >= _localMax)),
3065 "invalid local index = " << localID <<
" (local size = " <<
3068 Teuchos::Array< dim_type > localIndex = getLocalIndex(localID);
3069 size_type result = 0;
3070 for (
int axis = 0; axis < numDims(); ++axis)
3072 dim_type globalIndex = localIndex[axis] +
3073 _globalRankBounds[axis][getCommIndex(axis)].start() - _pad[axis][0];
3074 result += globalIndex * _globalStrides[axis];
3081 template<
class Node >
3086 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 3087 TEUCHOS_TEST_FOR_EXCEPTION(
3088 (globalIndex.size() != numDims()),
3090 "globalIndex has " << globalIndex.size() <<
" entries; expecting " 3092 for (
int axis = 0; axis < numDims(); ++axis)
3094 TEUCHOS_TEST_FOR_EXCEPTION(
3095 ((globalIndex[axis] < 0) ||
3096 (globalIndex[axis] >= _globalDims[axis])),
3098 "invalid globalIndex[" << axis <<
"] = " << globalIndex[axis] <<
3099 " (global dimension = " << _globalDims[axis] <<
")");
3102 size_type result = 0;
3103 for (
int axis = 0; axis < numDims(); ++axis)
3104 result += globalIndex[axis] * _globalStrides[axis];
3110 template<
class Node >
3115 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 3116 TEUCHOS_TEST_FOR_EXCEPTION(
3117 ((globalID < _globalMin) || (globalID >= _globalMax)),
3119 "invalid global index = " << globalID <<
" (should be between " <<
3120 _globalMin <<
" and " << _globalMax <<
")");
3122 Teuchos::Array< dim_type > globalIndex =
3123 getGlobalIndex(globalID);
3124 size_type result = 0;
3125 for (
int axis = 0; axis < numDims(); ++axis)
3127 dim_type localIndex = globalIndex[axis] -
3128 _globalRankBounds[axis][getCommIndex(axis)].start() + _pad[axis][0];
3129 TEUCHOS_TEST_FOR_EXCEPTION(
3130 (localIndex < 0 || localIndex >= _localDims[axis]),
3132 "global index not on local processor")
3133 result += localIndex * _localStrides[axis];
3140 template<
class Node >
3143 getLocalID(
const Teuchos::ArrayView< dim_type > & localIndex)
const 3145 #ifdef HAVE_DOMI_ARRAY_BOUNDSCHECK 3146 TEUCHOS_TEST_FOR_EXCEPTION(
3147 (localIndex.size() != numDims()),
3149 "localIndex has " << localIndex.size() <<
" entries; expecting " 3151 for (
int axis = 0; axis < numDims(); ++axis)
3153 TEUCHOS_TEST_FOR_EXCEPTION(
3154 ((localIndex[axis] < 0) ||
3155 (localIndex[axis] >= _localDims[axis])),
3157 "invalid localIndex[" << axis <<
"] = " << localIndex[axis] <<
3158 " (local dimension = " << _localDims[axis] <<
")");
3161 size_type result = 0;
3162 for (
int axis = 0; axis < numDims(); ++axis)
3163 result += localIndex[axis] * _localStrides[axis];
3169 template<
class Node >
3175 if (
this == &mdMap)
return true;
3179 int num_dims = numDims();
3180 if (num_dims != mdMap.
numDims())
return false;
3184 for (
int axis = 0; axis < num_dims; ++axis)
3185 if (getCommDim(axis) != mdMap.
getCommDim(axis))
return false;
3189 if (_globalDims != mdMap._globalDims)
return false;
3194 int localResult = 1;
3195 int globalResult = 1;
3196 for (
int axis = 0; axis < num_dims; ++axis)
3197 if (getLocalDim(axis,
false) != mdMap.
getLocalDim(axis,
false))
3199 Teuchos::reduceAll(*(getTeuchosComm()),
3200 Teuchos::REDUCE_MIN,
3206 return bool(globalResult);
3211 template<
class Node >
3214 const int verbose)
const 3218 if (
this == &mdMap)
return true;
3224 int localResult = 1;
3225 Teuchos::RCP< const Teuchos::Comm< int > > comm = getTeuchosComm();
3226 int rank = comm->getRank();
3229 if (! isCompatible(mdMap))
3233 std::cout << rank <<
": MDMaps are incompatible" << std::endl;
3241 std::cout << rank <<
": this Comm size = " << comm->getSize() <<
" != " 3250 std::cout << rank <<
": this Comm rank = " << rank <<
" != " 3255 if (_globalBounds != mdMap._globalBounds)
3259 std::cout << rank <<
": global bounds " << _globalBounds <<
" != " 3260 << mdMap._globalBounds << std::endl;
3264 if (_localDims != mdMap._localDims)
3268 std::cout << rank <<
": local dimensions " << _localDims <<
" != " 3269 << mdMap._localDims << std::endl;
3273 int globalResult = 1;
3274 Teuchos::reduceAll(*(getTeuchosComm()),
3275 Teuchos::REDUCE_MIN,
3281 return bool(globalResult);
3286 template<
class Node >
3291 Teuchos::Array< size_type > contiguousStrides =
3292 computeStrides< size_type, dim_type >(_localDims, _layout);
3295 int localResult = int(_localStrides != contiguousStrides);
3298 int globalResult = 0;
3299 Teuchos::reduceAll(*(_mdComm->getTeuchosComm()),
3300 Teuchos::REDUCE_SUM,
3304 return (globalResult == 0);
3309 template<
class Node >
3314 int num_dims = numDims();
3317 for (
int axis = 0; axis < num_dims; ++axis)
3320 int commDim = getCommDim(axis);
3321 for (
int axisRank = 0; axisRank < commDim; ++axisRank)
3326 dim_type localDim = (_globalDims[axis] - _bndryPad[axis][0] -
3327 _bndryPad[axis][1]) / commDim;
3328 dim_type axisStart = axisRank * localDim;
3337 dim_type remainder = (_globalDims[axis] - _bndryPad[axis][0] -
3338 _bndryPad[axis][1]) % commDim;
3339 if (commDim - axisRank - 1 < remainder)
3342 axisStart += (remainder - commDim + axisRank);
3346 axisStart += _bndryPad[axis][0];
3349 _globalRankBounds[axis].push_back(
3350 ConcreteSlice(axisStart, axisStart + localDim));
3354 if (axisRank == getCommIndex(axis))
3359 _localDims[axis] = localDim + _pad[axis][0] + _pad[axis][1];
3362 _localBounds.push_back(ConcreteSlice(_localDims[axis]));
Teuchos::RCP< const MDComm > getMDComm() const
Access the underlying MDComm object.
Definition: Domi_MDMap.hpp:2050
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:1043
Teuchos::Array< int > getCommDims() const
Get the communicator sizes along each axis.
Definition: Domi_MDMap.hpp:2086
size_type getGlobalID(size_type localID) const
Convert a local ID to a global ID.
Definition: Domi_MDMap.hpp:3059
Teuchos::RCP< Node > getNode() const
Get the Kokkos node.
Definition: Domi_MDMap.hpp:2520
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:2570
int getCommDim(int axis) const
Get the communicator size along the given axis.
Definition: Domi_MDMap.hpp:2095
MDMap< Node > & operator=(const MDMap< Node > &source)
Assignment operator.
Definition: Domi_MDMap.hpp:2024
int getLowerPadSize(int axis) const
Get the size of the lower padding along the given axis.
Definition: Domi_MDMap.hpp:2325
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:2983
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:2293
bool isReplicatedBoundary(int axis) const
Return whether the given axis supports a replicated boundary.
Definition: Domi_MDMap.hpp:2502
int numDims() const
Get the number of dimensions.
Definition: Domi_MDMap.hpp:2077
bool onSubcommunicator() const
Query whether this processor is on the sub-communicator.
Definition: Domi_MDMap.hpp:2059
dim_type getLocalDim(int axis, bool withPad=false) const
Get the local dimension along the specified axis.
Definition: Domi_MDMap.hpp:2197
bool isCompatible(const MDMap< Node > &mdMap) const
True if two MDMaps are "compatible".
Definition: Domi_MDMap.hpp:3171
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:2228
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:2104
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:2529
Teuchos::RCP< const Teuchos::Comm< int > > getTeuchosComm() const
Get the Teuchos communicator.
Definition: Domi_MDMap.hpp:2068
Teuchos::RCP< const Domi::MDMap< Node > > getAxisMap(int axis) const
Return an axis map for the given axis.
Definition: Domi_MDMap.hpp:2541
Teuchos::Array< int > getBndryPadSizes() const
Get the boundary padding sizes along each axis.
Definition: Domi_MDMap.hpp:2405
Teuchos::ArrayView< const Slice > getLocalBounds() const
Get the local loop bounds along every axis.
Definition: Domi_MDMap.hpp:2258
int getUpperNeighbor(int axis) const
Get the rank of the upper neighbor.
Definition: Domi_MDMap.hpp:2131
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:2357
bool isSameAs(const MDMap< Node > &mdMap, const int verbose=0) const
True if two MDMaps are "identical".
Definition: Domi_MDMap.hpp:3213
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:2341
bool isContiguous() const
True if there are no stride gaps on any processor.
Definition: Domi_MDMap.hpp:3288
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:2373
Slice getGlobalBounds(int axis, bool withBndryPad=false) const
Get the bounds of the global problem.
Definition: Domi_MDMap.hpp:2172
int getUpperBndryPad(int axis) const
Get the size of the upper boundary padding along the given axis.
Definition: Domi_MDMap.hpp:2389
~MDMap()
Destructor.
Definition: Domi_MDMap.hpp:2016
dim_type getGlobalDim(int axis, bool withBndryPad=false) const
Get the global dimension along the specified axis.
Definition: Domi_MDMap.hpp:2151
bool hasPadding() const
Return true if there is any padding stored locally.
Definition: Domi_MDMap.hpp:2313
int getBndryPadSize(int axis) const
Get the boundary padding size along the given axis.
Definition: Domi_MDMap.hpp:2414
size_type getLocalID(size_type globalID) const
Convert a global ID to a local ID.
Definition: Domi_MDMap.hpp:3113
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:2122
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:2141
Teuchos::Array< dim_type > getLocalIndex(size_type localID) const
Convert a local ID to a local index.
Definition: Domi_MDMap.hpp:3021
Teuchos::Array< dim_type > getLocalDims() const
Get an array of the local dimensions, including padding.
Definition: Domi_MDMap.hpp:2218
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:2113
Map Ordinal Error exception type.
Definition: Domi_Exceptions.hpp:90
Layout getLayout() const
Get the storage order.
Definition: Domi_MDMap.hpp:2511