43 #ifndef DOMI_MDVECTOR_HPP 44 #define DOMI_MDVECTOR_HPP 52 #include "Domi_ConfigDefs.hpp" 53 #include "Domi_DefaultNode.hpp" 54 #include "Domi_MDMap.hpp" 55 #include "Domi_MDArrayRCP.hpp" 58 #include "Teuchos_DataAccess.hpp" 59 #include "Teuchos_Describable.hpp" 60 #include "Teuchos_ScalarTraitsDecl.hpp" 61 #include "Teuchos_Comm.hpp" 62 #include "Teuchos_CommHelpers.hpp" 65 #include "Epetra_IntVector.h" 66 #include "Epetra_Vector.h" 70 #include "Tpetra_Vector.hpp" 81 struct ompi_datatype_t {};
173 template<
class Scalar,
174 class Node = DefaultNode::DefaultNodeType >
190 bool zeroOut =
true);
214 const dim_type leadingDim,
215 const dim_type trailingDim = 0,
216 bool zeroOut =
true);
251 Teuchos::DataAccess access = Teuchos::View);
266 MDVector(
const Teuchos::RCP<
const Teuchos::Comm< int > > teuchosComm,
267 Teuchos::ParameterList & plist);
281 MDVector(
const Teuchos::RCP< const MDComm > mdComm,
282 Teuchos::ParameterList & plist);
331 const Teuchos::ArrayView< Slice > & slices,
332 const Teuchos::ArrayView< int > & bndryPad =
333 Teuchos::ArrayView< int >());
353 const Teuchos::RCP< const MDMap< Node > >
370 Teuchos::RCP< const Teuchos::Comm< int > >
getTeuchosComm()
const;
458 dim_type
getGlobalDim(
int axis,
bool withBndryPad=
false)
const;
494 dim_type
getLocalDim(
int axis,
bool withCommPad=
false)
const;
668 Teuchos::RCP< Epetra_IntVector > getEpetraIntVectorView()
const;
685 Teuchos::RCP< Epetra_Vector > getEpetraVectorView()
const;
707 Teuchos::RCP< Epetra_MultiVector > getEpetraMultiVectorView()
const;
718 Teuchos::RCP< Epetra_IntVector > getEpetraIntVectorCopy()
const;
729 Teuchos::RCP< Epetra_Vector > getEpetraVectorCopy()
const;
747 Teuchos::RCP< Epetra_MultiVector > getEpetraMultiVectorCopy()
const;
763 template<
class LocalOrdinal,
764 class GlobalOrdinal = LocalOrdinal,
766 Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node2 > >
767 getTpetraVectorView()
const;
786 template <
class LocalOrdinal,
787 class GlobalOrdinal = LocalOrdinal,
789 Teuchos::RCP< Tpetra::MultiVector< Scalar,
793 getTpetraMultiVectorView()
const;
800 template<
class LocalOrdinal,
801 class GlobalOrdinal = LocalOrdinal,
803 Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node2 > >
804 getTpetraVectorCopy()
const;
818 template <
class LocalOrdinal,
819 class GlobalOrdinal = LocalOrdinal,
821 Teuchos::RCP< Tpetra::MultiVector< Scalar,
825 getTpetraMultiVectorCopy()
const;
890 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
norm1()
const;
894 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
norm2()
const;
898 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
normInf()
const;
904 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
928 describe(Teuchos::FancyOStream &out,
929 const Teuchos::EVerbosityLevel verbLevel =
930 Teuchos::Describable::verbLevel_default)
const;
945 bool includePadding =
true);
1045 bool includeBndryPad =
false)
const;
1054 void readBinary(
const std::string & filename,
1055 bool includeBndryPad =
false);
1064 Teuchos::RCP< const Teuchos::Comm< int > > _teuchosComm;
1068 Teuchos::RCP< const MDMap< Node > > _mdMap;
1089 Teuchos::Array< MPI_Request > _requests;
1102 Teuchos::RCP< MPI_Datatype > datatype;
1116 Teuchos::Array< Teuchos::Tuple< MessageInfo, 2 > > _sendMessages;
1121 Teuchos::Array< Teuchos::Tuple< MessageInfo, 2 > > _recvMessages;
1125 void initializeMessages();
1137 Teuchos::Array< dim_type > fileShape;
1138 Teuchos::Array< dim_type > bufferShape;
1139 Teuchos::Array< dim_type > dataShape;
1140 Teuchos::Array< dim_type > fileStart;
1141 Teuchos::Array< dim_type > dataStart;
1143 Teuchos::RCP< MPI_Datatype > filetype;
1144 Teuchos::RCP< MPI_Datatype > datatype;
1152 mutable Teuchos::RCP< FileInfo > _fileInfo;
1158 mutable Teuchos::RCP< FileInfo > _fileInfoWithBndry;
1166 Teuchos::RCP< FileInfo > & computeFileInfo(
bool includeBndryPad)
const;
1176 template<
class Scalar,
1181 _teuchosComm(mdMap->getTeuchosComm()),
1192 setObjectLabel(
"Domi::MDVector");
1195 int numDims = _mdMap->numDims();
1196 Teuchos::Array< dim_type > dims(
numDims);
1197 for (
int axis = 0; axis <
numDims; ++axis)
1198 dims[axis] = _mdMap->getLocalDim(axis,
true);
1201 _mdArrayRcp.
resize(dims);
1202 _mdArrayView = _mdArrayRcp();
1207 template<
class Scalar,
1211 const dim_type leadingDim,
1212 const dim_type trailingDim,
1214 _teuchosComm(mdMap->getTeuchosComm()),
1215 _mdMap(mdMap->getAugmentedMDMap(leadingDim, trailingDim)),
1225 setObjectLabel(
"Domi::MDVector");
1228 int numDims = _mdMap->numDims();
1229 Teuchos::Array< dim_type > dims(
numDims);
1230 for (
int axis = 0; axis <
numDims; ++axis)
1231 dims[axis] = _mdMap->getLocalDim(axis,
true);
1234 _mdArrayRcp.
resize(dims);
1235 _mdArrayView = _mdArrayRcp();
1240 template<
class Scalar,
1246 _mdArrayRcp(source),
1247 _mdArrayView(_mdArrayRcp()),
1255 setObjectLabel(
"Domi::MDVector");
1256 int numDims = _mdMap->numDims();
1257 TEUCHOS_TEST_FOR_EXCEPTION(
1260 "MDMap and source array do not have the same number of dimensions");
1262 for (
int axis = 0; axis <
numDims; ++axis)
1264 TEUCHOS_TEST_FOR_EXCEPTION(
1265 _mdMap->getLocalDim(axis) != _mdArrayRcp.
dimension(axis),
1267 "Axis " << axis <<
": MDMap dimension = " << _mdMap->getLocalDim(axis)
1268 <<
", MDArray dimension = " << _mdArrayRcp.
dimension(axis));
1274 template<
class Scalar,
1280 _mdArrayRcp(mdArrayRcp),
1281 _mdArrayView(_mdArrayRcp()),
1289 #ifdef DOMI_MDVECTOR_VERBOSE 1290 std::cout <<
"_mdArrayRcp = " << _mdArrayRcp << std::endl;
1291 std::cout <<
"_mdArrayView.getRawPtr() = " << _mdArrayView.
getRawPtr()
1292 <<
" (constructor)" << std::endl;
1293 std::cout <<
"_mdArrayView = " << _mdArrayView << std::endl;
1295 setObjectLabel(
"Domi::MDVector");
1296 int numDims = _mdMap->numDims();
1297 TEUCHOS_TEST_FOR_EXCEPTION(
1300 "MDMap and source array do not have the same number of dimensions");
1302 for (
int axis = 0; axis <
numDims; ++axis)
1304 TEUCHOS_TEST_FOR_EXCEPTION(
1305 _mdMap->getLocalDim(axis) != _mdArrayRcp.
dimension(axis),
1307 "Axis " << axis <<
": MDMap dimension = " << _mdMap->getLocalDim(axis)
1308 <<
", MDArray dimension = " << _mdArrayRcp.
dimension(axis));
1314 template<
class Scalar,
1318 Teuchos::DataAccess access) :
1319 _teuchosComm(source.getMDMap()->getTeuchosComm()),
1320 _mdMap(source.getMDMap()),
1321 _mdArrayRcp(source._mdArrayRcp),
1322 _mdArrayView(source._mdArrayView),
1330 setObjectLabel(
"Domi::MDVector");
1332 if (access == Teuchos::Copy)
1334 #ifdef DOMI_MDVECTOR_VERBOSE 1335 std::cout <<
"Inside MDVector copy constructor with copy access" << std::endl;
1338 int numDims = _mdMap->numDims();
1339 Teuchos::Array< dim_type > dims(
numDims);
1340 for (
int axis = 0; axis <
numDims; ++axis)
1341 dims[axis] = _mdMap->getLocalDim(axis,
true);
1345 _mdArrayView = _mdArrayRcp();
1350 const_iterator src = source.
getData().begin();
1351 for (iterator trg = _mdArrayView.
begin(); trg != _mdArrayView.
end(); ++trg)
1357 #ifdef DOMI_MDVECTOR_VERBOSE 1360 std::cout <<
"Inside MDVector copy constructor with view access" 1368 template<
class Scalar,
1371 MDVector(
const Teuchos::RCP<
const Teuchos::Comm< int > > teuchosComm,
1372 Teuchos::ParameterList & plist) :
1373 _teuchosComm(teuchosComm),
1384 setObjectLabel(
"Domi::MDVector");
1388 dim_type leadingDim = plist.get(
"leading dimension" , 0);
1389 dim_type trailingDim = plist.get(
"trailing dimension", 0);
1390 if (leadingDim + trailingDim > 0)
1396 _mdMap = Teuchos::rcp(myMdMap);
1399 int numDims = _mdMap->numDims();
1400 Teuchos::Array< dim_type > dims(
numDims);
1401 for (
int axis = 0; axis <
numDims; ++axis)
1402 dims[axis] = _mdMap->getLocalDim(axis,
true);
1405 _mdArrayRcp.
resize(dims);
1406 _mdArrayView = _mdArrayRcp();
1411 template<
class Scalar,
1415 Teuchos::ParameterList & plist) :
1416 _teuchosComm(mdComm->getTeuchosComm()),
1417 _mdMap(Teuchos::rcp(new
MDMap< Node >(mdComm, plist))),
1427 setObjectLabel(
"Domi::MDVector");
1431 dim_type leadingDim = plist.get(
"leading dimension" , 0);
1432 dim_type trailingDim = plist.get(
"trailing dimension", 0);
1433 if (leadingDim + trailingDim > 0)
1439 _mdMap = Teuchos::rcp(myMdMap);
1442 int numDims = _mdMap->numDims();
1443 Teuchos::Array< dim_type > dims(
numDims);
1444 for (
int axis = 0; axis <
numDims; ++axis)
1445 dims[axis] = _mdMap->getLocalDim(axis,
true);
1448 _mdArrayRcp.
resize(dims);
1449 _mdArrayView = _mdArrayRcp();
1454 template<
class Scalar,
1459 dim_type globalIndex) :
1460 _teuchosComm(parent._teuchosComm),
1462 _mdArrayRcp(parent._mdArrayRcp),
1463 _mdArrayView(parent._mdArrayView),
1471 setObjectLabel(
"Domi::MDVector");
1474 Teuchos::RCP< const MDMap< Node > > parentMdMap = parent.
getMDMap();
1482 if (_mdMap->onSubcommunicator())
1489 dim_type origin = parentMdMap->getGlobalRankBounds(axis,
false).start() -
1490 parentMdMap->getLowerPadSize(axis);
1495 dim_type localIndex = globalIndex - origin;
1499 _mdArrayView = newView;
1505 _mdArrayRcp.
clear();
1512 template<
class Scalar,
1517 const Slice & slice,
1521 _mdArrayRcp(parent._mdArrayRcp),
1522 _mdArrayView(parent._mdArrayView),
1530 #ifdef DOMI_MDVECTOR_VERBOSE 1531 std::cout <<
"slice axis " << axis << std::endl;
1532 std::cout <<
" _mdArrayRcp @ " << _mdArrayRcp.
getRawPtr() << std::endl;
1533 std::cout <<
" _mdArrayView @ " << _mdArrayView.
getRawPtr() << std::endl;
1536 setObjectLabel(
"Domi::MDVector");
1539 Teuchos::RCP< const MDMap< Node > > parentMdMap = parent.
getMDMap();
1546 _teuchosComm = _mdMap->getTeuchosComm();
1549 if (_mdMap->onSubcommunicator())
1552 Slice bounds = slice.
bounds(parentMdMap->getGlobalDim(axis,
true));
1559 dim_type origin = parentMdMap->getGlobalRankBounds(axis,
false).
start() -
1560 parentMdMap->getLowerPadSize(axis);
1567 dim_type start = std::max(0, bounds.
start() - origin - bndryPad);
1574 dim_type stop = std::min(bounds.
stop() - origin + bndryPad,
1575 parentMdMap->getLocalDim(axis,
true));
1579 _mdArrayView = newView;
1585 _mdArrayRcp.
clear();
1588 #ifdef DOMI_MDVECTOR_VERBOSE 1589 std::cout <<
" _mdArrayView @ " << _mdArrayView.
getRawPtr() << std::endl;
1590 std::cout <<
" offset = " << int(_mdArrayView.
getRawPtr() -
1597 template<
class Scalar,
1601 const Teuchos::ArrayView< Slice > & slices,
1602 const Teuchos::ArrayView< int > & bndryPad)
1604 setObjectLabel(
"Domi::MDVector");
1607 int numDims = parent.
numDims();
1610 TEUCHOS_TEST_FOR_EXCEPTION(
1611 (slices.size() != numDims),
1613 "number of slices = " << slices.size() <<
" != parent MDVector number of " 1614 "dimensions = " << numDims);
1618 for (
int axis = 0; axis < numDims; ++axis)
1620 int bndryPadding = (axis < bndryPad.size()) ? bndryPad[axis] : 0;
1625 tempMDVector1 = tempMDVector2;
1627 *
this = tempMDVector1;
1632 template<
class Scalar,
1638 _teuchosComm = source._teuchosComm;
1639 _mdMap = source._mdMap;
1640 _mdArrayRcp = source._mdArrayRcp;
1641 _mdArrayView = source._mdArrayView;
1642 _nextAxis = source._nextAxis;
1644 _requests = source._requests;
1646 _sendMessages = source._sendMessages;
1647 _recvMessages = source._recvMessages;
1653 template<
class Scalar,
1661 template<
class Scalar,
1663 const Teuchos::RCP< const MDMap< Node > >
1672 template<
class Scalar,
1678 return _mdMap->onSubcommunicator();
1683 template<
class Scalar,
1685 Teuchos::RCP< const Teuchos::Comm< int > >
1689 return _mdMap->getTeuchosComm();
1694 template<
class Scalar,
1700 return _mdMap->numDims();
1705 template<
class Scalar,
1711 return _mdMap->getCommDim(axis);
1716 template<
class Scalar,
1722 return _mdMap->isPeriodic(axis);
1727 template<
class Scalar,
1733 return _mdMap->getCommIndex(axis);
1738 template<
class Scalar,
1744 return _mdMap->getLowerNeighbor(axis);
1749 template<
class Scalar,
1755 return _mdMap->getUpperNeighbor(axis);
1760 template<
class Scalar,
1766 return _mdMap->getGlobalDim(axis, withBndryPad);
1771 template<
class Scalar,
1777 return _mdMap->getGlobalBounds(axis, withBndryPad);
1782 template<
class Scalar,
1788 return _mdMap->getGlobalRankBounds(axis, withBndryPad);
1793 template<
class Scalar,
1799 return _mdMap->getLocalDim(axis, withCommPad);
1804 template<
class Scalar,
1806 Teuchos::ArrayView< const Slice >
1810 return _mdMap->getLocalBounds();
1815 template<
class Scalar,
1821 return _mdMap->getLocalBounds(axis, withCommPad);
1826 template<
class Scalar,
1832 return _mdMap->getLocalInteriorBounds(axis);
1837 template<
class Scalar,
1843 return _mdMap->hasPadding();
1848 template<
class Scalar,
1854 return _mdMap->getLowerPadSize(axis);
1859 template<
class Scalar,
1865 return _mdMap->getUpperPadSize(axis);
1870 template<
class Scalar,
1876 return _mdMap->getCommPadSize(axis);
1881 template<
class Scalar,
1887 return _mdMap->getLowerBndryPad(axis);
1892 template<
class Scalar,
1898 return _mdMap->getUpperBndryPad(axis);
1903 template<
class Scalar,
1909 return _mdMap->getBndryPadSize(axis);
1914 template<
class Scalar,
1927 template<
class Scalar,
1940 template<
class Scalar,
1946 return _mdMap->isReplicatedBoundary(axis);
1951 template<
class Scalar,
1957 return _mdMap->getLayout();
1962 template<
class Scalar,
1968 return _mdMap->isContiguous();
1982 template<
class Scalar,
1984 Teuchos::RCP< Epetra_IntVector >
1989 TEUCHOS_TEST_FOR_EXCEPTION(
1990 typeid(Scalar) !=
typeid(
int),
1992 "MDVector is of scalar type '" <<
typeid(Scalar).name() <<
"', but " 1993 "Epetra_IntVector requires scalar type 'int'");
1996 TEUCHOS_TEST_FOR_EXCEPTION(
1999 "This MDVector's MDMap is non-contiguous. This can happen when you take " 2000 "a slice of a parent MDVector.");
2003 Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(
true);
2008 void * buffer = (
void*) _mdArrayView.getRawPtr();
2009 return Teuchos::rcp(
new Epetra_IntVector(View,
2023 template<
class Scalar,
2025 Teuchos::RCP< Epetra_Vector >
2026 MDVector< Scalar, Node >::
2027 getEpetraVectorView()
const 2030 TEUCHOS_TEST_FOR_EXCEPTION(
2031 typeid(Scalar) !=
typeid(
double),
2033 "MDVector is of scalar type '" <<
typeid(Scalar).name() <<
"', but " 2034 "Epetra_Vector requires scalar type 'double'");
2037 TEUCHOS_TEST_FOR_EXCEPTION(
2039 MDMapNoncontiguousError,
2040 "This MDVector's MDMap is non-contiguous. This can happen when you take " 2041 "a slice of a parent MDVector.");
2044 Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(
true);
2049 void * buffer = (
void*) _mdArrayView.getRawPtr();
2050 return Teuchos::rcp(
new Epetra_Vector(View,
2064 template<
class Scalar,
2066 Teuchos::RCP< Epetra_MultiVector >
2067 MDVector< Scalar, Node >::
2068 getEpetraMultiVectorView()
const 2071 TEUCHOS_TEST_FOR_EXCEPTION(
2072 typeid(Scalar) !=
typeid(
double),
2074 "MDVector is of scalar type '" <<
typeid(Scalar).name() <<
"', but " 2075 "Epetra_Vector requires scalar type 'double'");
2078 int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
2079 int padding = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
2080 int commDim = getCommDim(vectorAxis);
2081 int numVectors = getGlobalDim(vectorAxis);
2084 Teuchos::RCP< const MDMap< Node > > newMdMap;
2085 if (padding == 0 && commDim == 1)
2086 newMdMap = Teuchos::rcp(
new MDMap< Node >(*_mdMap, vectorAxis, 0));
2092 TEUCHOS_TEST_FOR_EXCEPTION(
2093 ! newMdMap->isContiguous(),
2094 MDMapNoncontiguousError,
2095 "This MDVector's MDMap is non-contiguous. This can happen when you take " 2096 "a slice of a parent MDVector.");
2101 size_type stride = newMdMap->getLocalDim(0,
true);
2102 for (
int axis = 0; axis < newMdMap->numDims(); ++axis)
2103 stride *= newMdMap->getLocalDim(axis,
true);
2104 TEUCHOS_TEST_FOR_EXCEPTION(
2105 stride*numVectors > Teuchos::OrdinalTraits<int>::max(),
2107 "Buffer size " << stride*numVectors <<
" is too large for Epetra int " 2109 int lda = (int)stride;
2112 Teuchos::RCP< const Epetra_Map > epetraMap = newMdMap->getEpetraMap(
true);
2117 void * buffer = (
void*) _mdArrayView.getRawPtr();
2118 return Teuchos::rcp(
new Epetra_MultiVector(View,
2127 template<
class Scalar,
2129 Teuchos::RCP< Epetra_IntVector >
2130 MDVector< Scalar, Node >::
2131 getEpetraIntVectorCopy()
const 2133 typedef typename MDArrayView< const Scalar >::iterator iterator;
2136 Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(
true);
2139 Teuchos::RCP< Epetra_IntVector > result =
2140 Teuchos::rcp(
new Epetra_IntVector(*epetraMap));
2145 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2146 result->operator[](ii++) = (int) *it;
2154 template<
class Scalar,
2156 Teuchos::RCP< Epetra_Vector >
2157 MDVector< Scalar, Node >::
2158 getEpetraVectorCopy()
const 2160 typedef typename MDArrayView< const Scalar >::iterator iterator;
2163 Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(
true);
2166 Teuchos::RCP< Epetra_Vector > result =
2167 Teuchos::rcp(
new Epetra_Vector(*epetraMap));
2172 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2173 result->operator[](ii++) = (double) *it;
2181 template<
class Scalar,
2183 Teuchos::RCP< Epetra_MultiVector >
2184 MDVector< Scalar, Node >::
2185 getEpetraMultiVectorCopy()
const 2187 typedef typename MDArrayView< Scalar >::iterator iterator;
2188 typedef typename MDArrayView< const Scalar >::iterator citerator;
2191 int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
2192 int padding = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
2193 int commDim = getCommDim(vectorAxis);
2194 int numVectors = getGlobalDim(vectorAxis);
2197 Teuchos::RCP< const MDMap< Node > > newMdMap;
2198 if (padding == 0 && commDim == 1)
2199 newMdMap = Teuchos::rcp(
new MDMap< Node >(*_mdMap, vectorAxis, 0));
2207 Teuchos::RCP< const Epetra_Map > epetraMap = newMdMap->getEpetraMap(
true);
2210 Teuchos::RCP< Epetra_MultiVector > result =
2211 Teuchos::rcp(
new Epetra_MultiVector(*epetraMap, numVectors));
2216 if (numVectors == 1)
2218 for (citerator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2219 result->operator[](0)[ii++] = (double) *it;
2223 for (
int iv = 0; iv < numVectors; ++iv)
2227 for (iterator it = subArray.begin(); it != subArray.end(); ++it)
2228 result->operator[](iv)[ii++] = (double) *it;
2242 template<
class Scalar,
class Node >
2243 template<
class LocalOrdinal,
2244 class GlobalOrdinal,
2246 Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node2 > >
2247 MDVector< Scalar, Node >::
2248 getTpetraVectorView()
const 2251 TEUCHOS_TEST_FOR_EXCEPTION(
2253 MDMapNoncontiguousError,
2254 "This MDVector's MDMap is non-contiguous. This can happen when you take " 2255 "a slice of a parent MDVector.");
2258 Teuchos::RCP<
const Tpetra::Map< LocalOrdinal,
2260 Node2 > > tpetraMap =
2261 _mdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal, Node2 >(
true);
2264 return Teuchos::rcp(
new Tpetra::Vector< Scalar,
2268 _mdArrayView.arrayView()));
2273 template<
class Scalar,
class Node >
2274 template<
class LocalOrdinal,
2275 class GlobalOrdinal,
2277 Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node2 > >
2278 MDVector< Scalar, Node >::
2279 getTpetraVectorCopy()
const 2281 typedef typename MDArrayView< const Scalar >::iterator iterator;
2284 Teuchos::RCP<
const Tpetra::Map< LocalOrdinal,
2286 Node2 > > tpetraMap =
2287 _mdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal, Node2 >(
true);
2290 Teuchos::RCP< Tpetra::Vector< Scalar,
2294 Teuchos::rcp(
new Tpetra::Vector< Scalar,
2297 Node2 >(tpetraMap) );
2301 Teuchos::ArrayRCP< Scalar > tpetraVectorArray = result->getDataNonConst();
2303 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2304 tpetraVectorArray[ii++] = *it;
2312 template<
class Scalar,
class Node >
2313 template <
class LocalOrdinal,
2314 class GlobalOrdinal,
2316 Teuchos::RCP< Tpetra::MultiVector< Scalar,
2320 MDVector< Scalar, Node >::
2321 getTpetraMultiVectorView()
const 2324 int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
2325 int padding = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
2326 int commDim = getCommDim(vectorAxis);
2327 size_t numVectors = getGlobalDim(vectorAxis);
2330 Teuchos::RCP< const MDMap< Node > > newMdMap;
2331 if (padding == 0 && commDim == 1)
2332 newMdMap = Teuchos::rcp(
new MDMap< Node >(*_mdMap, vectorAxis, 0));
2338 TEUCHOS_TEST_FOR_EXCEPTION(
2339 ! newMdMap->isContiguous(),
2340 MDMapNoncontiguousError,
2341 "This MDVector's MDMap is non-contiguous. This can happen when you take " 2342 "a slice of a parent MDVector.");
2347 size_type stride = newMdMap->getLocalDim(0,
true);
2348 for (
int axis = 0; axis < newMdMap->numDims(); ++axis)
2349 stride *= newMdMap->getLocalDim(axis,
true);
2350 TEUCHOS_TEST_FOR_EXCEPTION(
2351 stride*numVectors > Teuchos::OrdinalTraits<GlobalOrdinal>::max(),
2353 "Buffer size " << stride*numVectors <<
" is too large for Tpetra " 2354 "GlobalOrdinal = " <<
typeid(GlobalOrdinal).name() );
2355 size_t lda = (size_t)stride;
2358 Teuchos::RCP<
const Tpetra::Map< LocalOrdinal,
2360 Node2> > tpetraMap =
2361 newMdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal, Node2 >(
true);
2364 return Teuchos::rcp(
new Tpetra::MultiVector< Scalar,
2368 _mdArrayView.arrayView(),
2375 template<
class Scalar,
class Node >
2376 template <
class LocalOrdinal,
2377 class GlobalOrdinal,
2379 Teuchos::RCP< Tpetra::MultiVector< Scalar,
2383 MDVector< Scalar, Node >::
2384 getTpetraMultiVectorCopy()
const 2386 typedef typename MDArrayView< Scalar >::iterator iterator;
2387 typedef typename MDArrayView< const Scalar >::iterator citerator;
2390 int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
2391 int padding = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
2392 int commDim = getCommDim(vectorAxis);
2393 int numVectors = getGlobalDim(vectorAxis);
2396 Teuchos::RCP< const MDMap< Node > > newMdMap;
2397 if (padding == 0 && commDim == 1)
2398 newMdMap = Teuchos::rcp(
new MDMap< Node >(*_mdMap, vectorAxis, 0));
2406 Teuchos::RCP<
const Tpetra::Map< LocalOrdinal,
2408 Node2 > > tpetraMap =
2409 newMdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal, Node2 >(
true);
2412 Teuchos::RCP< Tpetra::MultiVector< Scalar,
2416 Teuchos::rcp(
new Tpetra::MultiVector< Scalar,
2419 Node2 >(tpetraMap, numVectors));
2424 if (numVectors == 1)
2426 Teuchos::ArrayRCP< Scalar > tpetraVectorArray = result->getDataNonConst(0);
2427 for (citerator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2428 tpetraVectorArray[ii++] = (
double) *it;
2432 for (
int iv = 0; iv < numVectors; ++iv)
2435 Teuchos::ArrayRCP< Scalar > tpetraVectorArray =
2436 result->getDataNonConst(iv);
2438 for (iterator it = subArray.begin(); it != subArray.end(); ++it)
2439 tpetraVectorArray[ii++] = (
double) *it;
2451 template<
class Scalar,
2457 #ifdef DOMI_MDVECTOR_VERBOSE 2458 std::cout <<
"_mdArrayView.getRawPtr() = " << _mdArrayView.
getRawPtr()
2460 std::cout <<
"_mdArrayView = " << _mdArrayView << std::endl;
2463 return _mdArrayView;
2465 for (
int axis = 0; axis < numDims(); ++axis)
2467 int lo = getLowerPadSize(axis);
2468 int hi = getLocalDim(axis,
true) - getUpperPadSize(axis);
2476 template<
class Scalar,
2483 return _mdArrayView.getConst();
2485 for (
int axis = 0; axis < numDims(); ++axis)
2487 int lo = getLowerPadSize(axis);
2488 int hi = getLocalDim(axis,
true) - getUpperPadSize(axis);
2496 template<
class Scalar,
2504 Slice(getLowerPadSize(axis)));
2505 return newArrayView;
2510 template<
class Scalar,
2518 Slice(getLowerPadSize(axis)));
2519 return newArrayView;
2524 template<
class Scalar,
2530 dim_type n = getLocalDim(axis,
true);
2531 int pad = getUpperPadSize(axis);
2533 if (pad) slice =
Slice(n-pad,n);
2534 else slice =
Slice(n-1,n-1);
2538 return newArrayView;
2543 template<
class Scalar,
2549 dim_type n = getLocalDim(axis,
true);
2550 int pad = getUpperPadSize(axis);
2552 if (pad) slice =
Slice(n-pad,n);
2553 else slice =
Slice(n-1,n-1);
2557 return newArrayView;
2562 template<
class Scalar,
2570 TEUCHOS_TEST_FOR_EXCEPTION(
2571 ! _mdMap->isCompatible(*(a._mdMap)),
2573 "MDMap of calling MDVector and argument 'a' are incompatible");
2576 Scalar local_dot = 0;
2577 iterator a_it = aView.begin();
2578 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end();
2580 local_dot += *it * *a_it;
2581 Scalar global_dot = 0;
2582 Teuchos::reduceAll(*_teuchosComm,
2583 Teuchos::REDUCE_SUM,
2592 template<
class Scalar,
2594 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
2598 typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
2601 mag local_norm1 = 0;
2602 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2603 local_norm1 += std::abs(*it);
2604 mag global_norm1 = 0;
2605 Teuchos::reduceAll(*_teuchosComm,
2606 Teuchos::REDUCE_SUM,
2610 return global_norm1;
2615 template<
class Scalar,
2617 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
2621 typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
2622 mag norm2 = dot(*
this);
2623 return Teuchos::ScalarTraits<mag>::squareroot(norm2);
2628 template<
class Scalar,
2630 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
2634 typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
2637 mag local_normInf = 0;
2638 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2639 local_normInf = std::max(local_normInf, std::abs(*it));
2640 mag global_normInf = 0;
2641 Teuchos::reduceAll(*_teuchosComm,
2642 Teuchos::REDUCE_MAX,
2646 return global_normInf;
2651 template<
class Scalar,
2653 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
2657 typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
2660 TEUCHOS_TEST_FOR_EXCEPTION(
2661 ! _mdMap->isCompatible(*(weights._mdMap)),
2663 "MDMap of calling MDVector and argument 'weights' are incompatible");
2666 mag local_wNorm = 0;
2667 iterator w_it = wView.begin();
2668 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end();
2670 local_wNorm += *it * *it * *w_it;
2671 mag global_wNorm = 0;
2672 Teuchos::reduceAll(*_teuchosComm,
2673 Teuchos::REDUCE_SUM,
2677 Teuchos::Array< dim_type > dimensions(numDims());
2678 for (
int i = 0; i < numDims(); ++i)
2679 dimensions[i] = _mdMap->getGlobalDim(i);
2680 size_type n = computeSize(dimensions);
2681 if (n == 0)
return 0;
2682 return Teuchos::ScalarTraits<mag>::squareroot(global_wNorm / n);
2687 template<
class Scalar,
2693 typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
2697 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2700 Teuchos::reduceAll(*_teuchosComm,
2701 Teuchos::REDUCE_SUM,
2705 Teuchos::Array< dim_type > dimensions(numDims());
2706 for (
int i = 0; i < numDims(); ++i)
2707 dimensions[i] = _mdMap->getGlobalDim(i);
2708 size_type n = computeSize(dimensions);
2709 if (n == 0)
return 0;
2710 return global_sum / n;
2715 template<
class Scalar,
2721 using Teuchos::TypeNameTraits;
2723 Teuchos::Array< dim_type > dims(numDims());
2724 for (
int axis = 0; axis < numDims(); ++axis)
2725 dims[axis] = getGlobalDim(axis,
true);
2727 std::ostringstream oss;
2728 oss <<
"\"Domi::MDVector\": {" 2729 <<
"Template parameters: {" 2730 <<
"Scalar: " << TypeNameTraits<Scalar>::name()
2731 <<
", Node: " << TypeNameTraits< Node >::name()
2733 if (this->getObjectLabel() !=
"")
2734 oss <<
", Label: \"" << this->getObjectLabel () <<
"\", ";
2735 oss <<
"Global dimensions: " << dims <<
" }";
2741 template<
class Scalar,
2746 const Teuchos::EVerbosityLevel verbLevel)
const 2750 using Teuchos::Comm;
2752 using Teuchos::TypeNameTraits;
2753 using Teuchos::VERB_DEFAULT;
2754 using Teuchos::VERB_NONE;
2755 using Teuchos::VERB_LOW;
2756 using Teuchos::VERB_MEDIUM;
2757 using Teuchos::VERB_HIGH;
2758 using Teuchos::VERB_EXTREME;
2760 const Teuchos::EVerbosityLevel vl =
2761 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
2764 Teuchos::RCP< const Teuchos::Comm< int > > comm = mdMap.
getTeuchosComm();
2765 const int myImageID = comm->getRank();
2766 const int numImages = comm->getSize();
2767 Teuchos::OSTab tab0(out);
2769 if (vl != VERB_NONE)
2773 out <<
"\"Domi::MDVector\":" << endl;
2775 Teuchos::OSTab tab1(out);
2778 out <<
"Template parameters:";
2780 Teuchos::OSTab tab2(out);
2781 out <<
"Scalar: " << TypeNameTraits<Scalar>::name() << endl
2782 <<
"Node: " << TypeNameTraits< Node >::name() << endl;
2785 if (this->getObjectLabel() !=
"")
2787 out <<
"Label: \"" << getObjectLabel() <<
"\"" << endl;
2789 Teuchos::Array< dim_type > globalDims(numDims());
2790 for (
int axis = 0; axis < numDims(); ++axis)
2791 globalDims[axis] = getGlobalDim(axis,
true);
2792 out <<
"Global dimensions: " << globalDims << endl;
2794 for (
int imageCtr = 0; imageCtr < numImages; ++imageCtr)
2796 if (myImageID == imageCtr)
2801 out <<
"Process: " << myImageID << endl;
2802 Teuchos::OSTab tab2(out);
2803 Teuchos::Array< dim_type > localDims(numDims());
2804 for (
int axis = 0; axis < numDims(); ++axis)
2805 localDims[axis] = getLocalDim(axis,
true);
2806 out <<
"Local dimensions: " << localDims << endl;
2819 template<
class Scalar,
2824 bool includePadding)
2829 for (iterator it = newArray.
begin(); it != newArray.
end(); ++it)
2835 template<
class Scalar,
2843 Teuchos::ScalarTraits< Scalar >::seedrandom(time(NULL));
2844 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2845 *it = Teuchos::ScalarTraits< Scalar >::random();
2850 template<
class Scalar,
2856 for (
int axis = 0; axis < numDims(); ++axis)
2858 updateCommPad(axis);
2864 template<
class Scalar,
2870 startUpdateCommPad(axis);
2871 endUpdateCommPad(axis);
2876 template<
class Scalar,
2886 if (_sendMessages.empty()) initializeMessages();
2889 int rank = _teuchosComm->getRank();
2890 int numProc = _teuchosComm->getSize();
2896 Teuchos::RCP< const Teuchos::MpiComm< int > > mpiComm =
2897 Teuchos::rcp_dynamic_cast<
const Teuchos::MpiComm< int > >(_teuchosComm);
2898 const Teuchos::OpaqueWrapper< MPI_Comm > & communicator =
2899 *(mpiComm->getRawMpiComm());
2902 MPI_Request request;
2903 for (
int boundary = 0; boundary < 2; ++boundary)
2905 MessageInfo message = _sendMessages[axis][boundary];
2906 if (message.proc >= 0)
2908 tag = 2 * (rank * numProc + message.proc) + boundary;
2910 #ifdef DOMI_MDVECTOR_OUTPUT_UPDATECOMMPAD 2911 cout << rank <<
": post send for axis " << axis <<
", boundary " 2912 << boundary <<
", buffer = " << message.buffer <<
", proc = " 2913 << message.proc <<
", tag = " << tag << endl;
2916 if (MPI_Isend(message.buffer,
2918 *(message.datatype),
2923 throw std::runtime_error(
"Domi::MDVector: Error in MPI_Isend");
2924 _requests.push_back(request);
2929 for (
int boundary = 0; boundary < 2; ++boundary)
2931 MessageInfo message = _recvMessages[axis][boundary];
2932 if (message.proc >= 0)
2934 tag = 2 * (message.proc * numProc + rank) + (1-boundary);
2936 #ifdef DOMI_MDVECTOR_OUTPUT_UPDATECOMMPAD 2937 cout << rank <<
": post recv for axis " << axis <<
", boundary " 2938 << boundary <<
", buffer = " << message.buffer <<
", proc = " 2939 << message.proc <<
", tag = " << tag << endl;
2942 if (MPI_Irecv(message.buffer,
2944 *(message.datatype),
2949 throw std::runtime_error(
"Domi::MDVector: Error in MPI_Irecv");
2950 _requests.push_back(request);
2957 if (isPeriodic(axis))
2959 for (
int sendBndry = 0; sendBndry < 2; ++sendBndry)
2961 int recvBndry = 1 - sendBndry;
2971 for ( ; it_recv != recvView.
end(); ++it_recv, ++it_send)
2972 *it_recv = *it_send;
2980 template<
class Scalar,
2987 if (_requests.size() > 0)
2989 Teuchos::Array< MPI_Status > status(_requests.size());
2990 if (MPI_Waitall(_requests.size(),
2993 throw std::runtime_error(
"Domi::MDVector: Error in MPI_Waitall");
3001 template<
class Scalar,
3010 int newAxis = _nextAxis + 1;
3011 if (newAxis >= result.
numDims()) newAxis = 0;
3012 result._nextAxis = newAxis;
3018 template<
class Scalar,
3027 int newAxis = _nextAxis + 1;
3028 if (newAxis >= result.
numDims()) newAxis = 0;
3029 result._nextAxis = newAxis;
3035 template<
class Scalar,
3043 int ndims = numDims();
3044 Teuchos::Array<int> sizes(ndims);
3045 Teuchos::Array<int> subsizes(ndims);
3046 Teuchos::Array<int> starts(ndims);
3047 MessageInfo messageInfo;
3049 _sendMessages.resize(ndims);
3050 _recvMessages.resize(ndims);
3052 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3053 std::stringstream msg;
3054 int rank = _teuchosComm->getRank();
3058 int order = mpiOrder(getLayout());
3059 MPI_Datatype datatype = mpiType< Scalar >();
3063 for (
int msgAxis = 0; msgAxis < ndims; ++msgAxis)
3066 for (
int axis = 0; axis < ndims; ++axis)
3068 sizes[axis] = _mdArrayRcp.dimension(axis);
3069 subsizes[axis] = _mdArrayView.dimension(axis);
3077 int proc = getLowerNeighbor(msgAxis);
3079 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3080 msg << endl <<
"P" << rank <<
": axis " << msgAxis <<
", lower neighbor = " 3085 subsizes[msgAxis] = getLowerPadSize(msgAxis);
3087 if (subsizes[msgAxis] == 0) proc = -1;
3089 messageInfo.buffer = (
void*) getData().getRawPtr();
3090 messageInfo.proc = proc;
3091 messageInfo.axis = msgAxis;
3101 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3102 msg << endl <<
"P" << rank <<
": axis " << msgAxis
3103 <<
", Lower receive message:" << endl <<
" ndims = " << ndims
3104 << endl <<
" sizes = " << sizes << endl <<
" subsizes = " 3105 << subsizes << endl <<
" starts = " << starts << endl
3106 <<
" order = " << order << endl;
3108 Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(
new MPI_Datatype);
3109 MPI_Type_create_subarray(ndims,
3116 MPI_Type_commit(commPad.get());
3117 messageInfo.datatype = commPad;
3119 messageInfo.dataview = _mdArrayView;
3120 for (
int axis = 0; axis < numDims(); ++axis)
3122 Slice slice(starts[axis], starts[axis] + subsizes[axis]);
3130 _recvMessages[msgAxis][0] = messageInfo;
3136 starts[msgAxis] = getLowerPadSize(msgAxis);
3137 if (isReplicatedBoundary(msgAxis) && getCommIndex(msgAxis) == 0)
3138 starts[msgAxis] += 1;
3144 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3145 msg << endl <<
"P" << rank <<
": axis " << msgAxis
3146 <<
", Lower send message:" << endl <<
" ndims = " << ndims
3147 << endl <<
" sizes = " << sizes << endl <<
" subsizes = " 3148 << subsizes << endl <<
" starts = " << starts << endl
3149 <<
" order = " << order << endl;
3151 Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(
new MPI_Datatype);
3152 MPI_Type_create_subarray(ndims,
3159 MPI_Type_commit(commPad.get());
3160 messageInfo.datatype = commPad;
3162 messageInfo.dataview = _mdArrayView;
3163 for (
int axis = 0; axis < numDims(); ++axis)
3165 Slice slice(starts[axis], starts[axis] + subsizes[axis]);
3173 _sendMessages[msgAxis][0] = messageInfo;
3179 proc = getUpperNeighbor(msgAxis);
3181 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3182 msg << endl <<
"P" << rank <<
": axis " << msgAxis <<
", upper neighbor = " 3186 subsizes[msgAxis] = getUpperPadSize(msgAxis);
3187 starts[msgAxis] = _mdArrayView.dimension(msgAxis) -
3188 getUpperPadSize(msgAxis);
3189 if (subsizes[msgAxis] == 0) proc = -1;
3190 messageInfo.proc = proc;
3199 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3200 msg << endl <<
"P" << rank <<
": axis " << msgAxis
3201 <<
", Upper receive message:" << endl <<
" ndims = " << ndims
3202 << endl <<
" sizes = " << sizes << endl <<
" subsizes = " 3203 << subsizes << endl <<
" starts = " << starts << endl
3204 <<
" order = " << order << endl;
3206 Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(
new MPI_Datatype);
3207 MPI_Type_create_subarray(ndims,
3214 MPI_Type_commit(commPad.get());
3215 messageInfo.datatype = commPad;
3217 messageInfo.dataview = _mdArrayView;
3218 for (
int axis = 0; axis < numDims(); ++axis)
3220 Slice slice(starts[axis], starts[axis] + subsizes[axis]);
3227 _recvMessages[msgAxis][1] = messageInfo;
3233 starts[msgAxis] -= getUpperPadSize(msgAxis);
3234 if (isReplicatedBoundary(msgAxis) &&
3235 getCommIndex(msgAxis) == getCommDim(msgAxis)-1)
3236 starts[msgAxis] -= 1;
3242 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3243 msg << endl <<
"P" << rank <<
": axis " << msgAxis
3244 <<
", Upper send message:" << endl <<
" ndims = " << ndims
3245 << endl <<
" sizes = " << sizes << endl <<
" subsizes = " 3246 << subsizes << endl <<
" starts = " << starts << endl
3247 <<
" order = " << order << endl;
3249 Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(
new MPI_Datatype);
3250 MPI_Type_create_subarray(ndims,
3257 MPI_Type_commit(commPad.get());
3258 messageInfo.datatype = commPad;
3260 messageInfo.dataview = _mdArrayView;
3261 for (
int axis = 0; axis < numDims(); ++axis)
3263 Slice slice(starts[axis], starts[axis] + subsizes[axis]);
3271 _sendMessages[msgAxis][1] = messageInfo;
3274 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3275 for (
int proc = 0; proc < _teuchosComm->getSize(); ++proc)
3282 _teuchosComm->barrier();
3290 template<
class Scalar,
3295 bool includeBndryPad)
const 3302 int pid = _teuchosComm->getRank();
3305 datafile = fopen(filename.c_str(),
"w");
3308 _teuchosComm->barrier();
3311 const Scalar * buffer = getData(
true).getRawPtr();
3315 Teuchos::RCP< FileInfo > & fileInfo = computeFileInfo(includeBndryPad);
3324 Teuchos::RCP< const Teuchos::MpiComm< int > > mpiComm =
3325 Teuchos::rcp_dynamic_cast<
const Teuchos::MpiComm< int > >(_teuchosComm);
3326 const Teuchos::OpaqueWrapper< MPI_Comm > & communicator =
3327 *(mpiComm->getRawMpiComm());
3330 int access = MPI_MODE_WRONLY | MPI_MODE_CREATE;
3335 char * cstr =
new char[filename.size()+1];
3336 std::strcpy(cstr, filename.c_str());
3341 char datarep[7] =
"native";
3342 MPI_File_open(communicator(), cstr, access, MPI_INFO_NULL, &mpiFile);
3343 MPI_File_set_view(mpiFile, 0, mpiType< Scalar >(),
3344 *(fileInfo->filetype), datarep, MPI_INFO_NULL);
3345 MPI_File_write_all(mpiFile, (
void*)buffer, 1, *(fileInfo->datatype),
3347 MPI_File_close(&mpiFile);
3359 datafile = fopen(filename.c_str(),
"w");
3366 for (const_iterator it = mdArrayView.begin(); it != mdArrayView.end(); ++it)
3368 fwrite((
const void *) &(*it),
sizeof(Scalar), 1, datafile);
3380 template<
class Scalar,
3385 bool includeBndryPad)
3388 const Scalar * buffer = getDataNonConst(
true).getRawPtr();
3392 Teuchos::RCP< FileInfo > & fileInfo = computeFileInfo(includeBndryPad);
3401 Teuchos::RCP< const Teuchos::MpiComm< int > > mpiComm =
3402 Teuchos::rcp_dynamic_cast<
const Teuchos::MpiComm< int > >(_teuchosComm);
3403 const Teuchos::OpaqueWrapper< MPI_Comm > & communicator =
3404 *(mpiComm->getRawMpiComm());
3407 int access = MPI_MODE_RDONLY;
3412 char * cstr =
new char[filename.size()+1];
3413 std::strcpy(cstr, filename.c_str());
3418 char datarep[7] =
"native";
3419 MPI_File_open(communicator(), cstr, access, MPI_INFO_NULL, &mpiFile);
3420 MPI_File_set_view(mpiFile, 0, mpiType< Scalar >(),
3421 *(fileInfo->filetype), datarep, MPI_INFO_NULL);
3422 MPI_File_read_all(mpiFile, (
void*)buffer, 1, *(fileInfo->datatype),
3424 MPI_File_close(&mpiFile);
3433 int ndims = _mdMap->numDims();
3437 datafile = fopen(filename.c_str(),
"r");
3444 Teuchos::Array< Ordinal > index(3);
3445 for (
int axis = 0; axis < ndims; ++axis)
3446 index[axis] = fileInfo->dataStart[axis];
3450 for (iterator it = mdArrayView.begin(); it != mdArrayView.end(); ++it)
3452 fread(&(*it),
sizeof(Scalar), 1, datafile);
3464 template<
class Scalar,
3466 Teuchos::RCP< typename MDVector< Scalar, Node >::FileInfo > &
3472 Teuchos::RCP< MDVector< Scalar, Node >::FileInfo > & fileInfo =
3473 includeBndryPad ? _fileInfoWithBndry : _fileInfo;
3476 if (!fileInfo.is_null())
return fileInfo;
3479 int ndims = _mdMap->
numDims();
3481 fileInfo->fileShape.resize(ndims);
3482 fileInfo->bufferShape.resize(ndims);
3483 fileInfo->dataShape.resize(ndims);
3484 fileInfo->fileStart.resize(ndims);
3485 fileInfo->dataStart.resize(ndims);
3488 for (
int axis = 0; axis < ndims; ++axis)
3492 fileInfo->fileShape[axis] = getGlobalDim(axis,includeBndryPad);
3493 fileInfo->bufferShape[axis] = getLocalDim(axis,
true );
3494 fileInfo->dataShape[axis] = getLocalDim(axis,
false);
3495 fileInfo->fileStart[axis] = getGlobalRankBounds(axis,includeBndryPad).start();
3496 fileInfo->dataStart[axis] = getLocalBounds(axis).start();
3498 if (includeBndryPad)
3500 int commIndex = _mdMap->getCommIndex(axis);
3503 int pad = getLowerBndryPad(axis);
3504 fileInfo->dataShape[axis] += pad;
3505 fileInfo->dataStart[axis] -= pad;
3507 if (commIndex == _mdMap->getCommDim(axis)-1)
3509 fileInfo->dataShape[axis] += getUpperBndryPad(axis);
3514 #ifdef DOMI_MDVECTOR_DEBUG_IO 3515 cout << pid <<
": fileShape = " << fileInfo->fileShape() << endl;
3516 cout << pid <<
": bufferShape = " << fileInfo->bufferShape() << endl;
3517 cout << pid <<
": dataShape = " << fileInfo->dataShape() << endl;
3518 cout << pid <<
": fileStart = " << fileInfo->fileStart() << endl;
3519 cout << pid <<
": dataStart = " << fileInfo->dataStart() << endl;
3523 int mpi_order = getLayout() == C_ORDER ? MPI_ORDER_C : MPI_ORDER_FORTRAN;
3525 fileInfo->filetype = Teuchos::rcp(
new MPI_Datatype);
3526 MPI_Type_create_subarray(ndims,
3527 fileInfo->fileShape.getRawPtr(),
3528 fileInfo->dataShape.getRawPtr(),
3529 fileInfo->fileStart.getRawPtr(),
3531 mpiType< Scalar >(),
3532 fileInfo->filetype.get());
3533 MPI_Type_commit(fileInfo->filetype.get());
3536 fileInfo->datatype = Teuchos::rcp(
new MPI_Datatype);
3537 MPI_Type_create_subarray(ndims,
3538 fileInfo->bufferShape.getRawPtr(),
3539 fileInfo->dataShape.getRawPtr(),
3540 fileInfo->dataStart.getRawPtr(),
3542 mpiType< Scalar >(),
3543 fileInfo->datatype.get());
3544 MPI_Type_commit(fileInfo->datatype.get());
3545 #endif // DGM_PARALLEL void writeBinary(const std::string &filename, bool includeBndryPad=false) const
Write the MDVector to a binary file.
Definition: Domi_MDVector.hpp:3294
int getCommPadSize(int axis) const
Get the communication padding size along the given axis.
Definition: Domi_MDVector.hpp:1874
Scalar meanValue() const
Compute the mean (average) value of this MDVector.
Definition: Domi_MDVector.hpp:2691
MDVector< Scalar, Node > operator[](dim_type index) const
Sub-vector access operator.
Definition: Domi_MDVector.hpp:3005
virtual std::string description() const
A simple one-line description of this MDVector.
Definition: Domi_MDVector.hpp:2719
void resize(const Teuchos::ArrayView< dim_type > &dims)
Resize the MDArrayRCP based on the given dimensions.
Definition: Domi_MDArrayRCP.hpp:1684
int getUpperPadSize(int axis) const
Get the size of the upper padding along the given axis.
Definition: Domi_MDVector.hpp:1863
MDVector(const Teuchos::RCP< const MDMap< Node > > &mdMap, bool zeroOut=true)
Main constructor.
Definition: Domi_MDVector.hpp:1179
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to a FancyOStream.
Definition: Domi_MDVector.hpp:2745
int getLowerPadSize(int axis) const
Get the size of the lower padding along the given axis.
Definition: Domi_MDVector.hpp:1852
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
void assign(const T &value)
Assign a value to all elements of the MDArrayView
Definition: Domi_MDArrayView.hpp:1413
MDArrayView< Scalar > getLowerPadDataNonConst(int axis)
Get a non-const view of the lower padding data along the given axis as an MDArrayView.
Definition: Domi_MDVector.hpp:2500
Teuchos::ArrayView< const Slice > getLocalBounds() const
Get the local loop bounds along every axis.
Definition: Domi_MDVector.hpp:1808
void setLowerPad(int axis, const Scalar value)
Assign all elements of the lower pad a constant value.
Definition: Domi_MDVector.hpp:1918
Teuchos::ScalarTraits< Scalar >::magnitudeType norm1() const
Compute the 1-norm of this MDVector.
Definition: Domi_MDVector.hpp:2596
void setUpperPad(int axis, const Scalar value)
Assign all elements of the upper pad a constant value.
Definition: Domi_MDVector.hpp:1931
dim_type getLocalDim(int axis, bool withCommPad=false) const
Get the local dimension along the specified axis.
Definition: Domi_MDVector.hpp:1797
iterator begin()
Return the beginning iterator.
Definition: Domi_MDArrayView.hpp:960
MDArrayView< const Scalar > getLowerPadData(int axis) const
Get a const view of the lower padding data along the given axis as an MDArrayView.
Definition: Domi_MDVector.hpp:2514
bool onSubcommunicator() const
Query whether this processor is on the sub-communicator.
Definition: Domi_MDVector.hpp:1676
Iterator class suitable for multi-dimensional arrays.
Definition: Domi_MDIterator.hpp:100
iterator end()
Return the ending iterator.
Definition: Domi_MDArrayView.hpp:969
Scalar dot(const MDVector< Scalar, Node > &a) const
Compute the dot product of this MDVector and MDVector a.
Definition: Domi_MDVector.hpp:2566
void startUpdateCommPad(int axis)
Start an asyncronous update of the communication padding.
Definition: Domi_MDVector.hpp:2880
MDArrayView< const Scalar > getData(bool includePadding=true) const
Get a const view of the data as an MDArrayView.
Definition: Domi_MDVector.hpp:2480
bool isReplicatedBoundary(int axis) const
Return whether the given axis supports a replicated boundary.
Definition: Domi_MDVector.hpp:1944
int numDims() const
Get the number of dimensions.
Definition: Domi_MDVector.hpp:1698
Layout getLayout() const
Get the storage order.
Definition: Domi_MDVector.hpp:1955
dim_type getGlobalDim(int axis, bool withBndryPad=false) const
Get the global dimension along the specified axis.
Definition: Domi_MDVector.hpp:1764
int getCommDim(int axis) const
Get the communicator size along the given axis.
Definition: Domi_MDVector.hpp:1709
const T * getRawPtr() const
Return a const raw pointer to the beginning of the MDArrayView or NULL if unsized.
Definition: Domi_MDArrayView.hpp:1521
void readBinary(const std::string &filename, bool includeBndryPad=false)
Read the MDVector from a binary file.
Definition: Domi_MDVector.hpp:3384
int getCommIndex(int axis) const
Get the axis rank of this processor.
Definition: Domi_MDVector.hpp:1731
MDArrayView< Scalar > getUpperPadDataNonConst(int axis)
Get a non-const view of the upper padding data along the given axis as an MDArrayView.
Definition: Domi_MDVector.hpp:2528
const dim_type stop() const
Stop index.
Definition: Domi_Slice.hpp:438
bool hasPadding() const
Return true if there is any padding stored locally.
Definition: Domi_MDVector.hpp:1841
Teuchos::ScalarTraits< Scalar >::magnitudeType norm2() const
Compute the 2-norm of this MDVector.
Definition: Domi_MDVector.hpp:2619
iterator begin()
Return the beginning iterator.
Definition: Domi_MDArrayRCP.hpp:1065
A Slice contains a start, stop, and step index, describing a subset of an ordered container...
Definition: Domi_Slice.hpp:137
Teuchos::ScalarTraits< Scalar >::magnitudeType normInf() const
Compute the infinity-norm of this MDVector.
Definition: Domi_MDVector.hpp:2632
Teuchos::RCP< const Teuchos::Comm< int > > getTeuchosComm() const
Get the Teuchos communicator.
Definition: Domi_MDMap.hpp:2062
bool isPeriodic(int axis) const
Return the periodic flag for the given axis.
Definition: Domi_MDVector.hpp:1720
void endUpdateCommPad(int axis)
Complete an asyncronous update of the communication padding.
Definition: Domi_MDVector.hpp:2984
void clear()
Clear the MDArrayRCP
Definition: Domi_MDArrayRCP.hpp:1652
MDArrayView< const Scalar > getUpperPadData(int axis) const
Get a const view of the upper padding data along the given axis as an MDArrayView.
Definition: Domi_MDVector.hpp:2547
void randomize()
Set all values in the multivector to pseudorandom numbers.
Definition: Domi_MDVector.hpp:2839
MDMap Error exception type.
Definition: Domi_Exceptions.hpp:114
Definition: Domi_ConfigDefs.hpp:97
Type Error exception type.
Definition: Domi_Exceptions.hpp:126
MDVector< Scalar, Node > & operator=(const MDVector< Scalar, Node > &source)
Assignment operator.
Definition: Domi_MDVector.hpp:1636
virtual Slice bounds(dim_type size) const
Return a Slice with concrete start and stop values.
virtual ~MDVector()
Destructor.
Definition: Domi_MDVector.hpp:1655
Slice getLocalInteriorBounds(int axis) const
Get the local interior looping bounds along the specified axis.
Definition: Domi_MDVector.hpp:1830
void putScalar(const Scalar &value, bool includePadding=true)
Set all values in the multivector with the given value.
Definition: Domi_MDVector.hpp:2823
int numDims() const
Return the number of dimensions.
Definition: Domi_MDArrayRCP.hpp:999
int getBndryPadSize(int axis) const
Get the boundary padding size along the given axis.
Definition: Domi_MDVector.hpp:1907
Teuchos::ScalarTraits< Scalar >::magnitudeType normWeighted(const MDVector< Scalar, Node > &weights) const
Compute the weighted norm of this.
Definition: Domi_MDVector.hpp:2655
Slice getGlobalRankBounds(int axis, bool withBndryPad=false) const
Get the global loop bounds on this processor along the specified axis.
Definition: Domi_MDVector.hpp:1786
Multi-dimensional distributed vector.
Definition: Domi_MDVector.hpp:175
const_pointer getRawPtr() const
Return a const raw pointer to the beginning of the MDArrayRCP or NULL if unsized. ...
Definition: Domi_MDArrayRCP.hpp:1718
Slice getGlobalBounds(int axis, bool withBndryPadding=false) const
Get the bounds of the global problem.
Definition: Domi_MDVector.hpp:1775
Teuchos::RCP< const Teuchos::Comm< int > > getTeuchosComm() const
Get the Teuchos communicator.
Definition: Domi_MDVector.hpp:1687
void updateCommPad()
Sum values of a locally replicated multivector across all processes.
Definition: Domi_MDVector.hpp:2854
iterator end()
Return the ending iterator.
Definition: Domi_MDArrayRCP.hpp:1074
int getLowerNeighbor(int axis) const
Get the rank of the lower neighbor.
Definition: Domi_MDVector.hpp:1742
int getLowerBndryPad(int axis) const
Get the size of the lower boundary padding along the given axis.
Definition: Domi_MDVector.hpp:1885
MDArrayView< const T > getConst() const
Return an MDArrayView< const T > of an MDArrayView< T > object.
Definition: Domi_MDArrayView.hpp:1055
dim_type dimension(int axis) const
Return the dimension of the given axis.
Definition: Domi_MDArrayRCP.hpp:1017
const Teuchos::RCP< const MDMap< Node > > getMDMap() const
MDMap accessor method.
Definition: Domi_MDVector.hpp:1665
Multi-dimensional map.
Definition: Domi_MDMap.hpp:144
const dim_type start() const
Start index.
Definition: Domi_Slice.hpp:431
MDArrayView< Scalar > getDataNonConst(bool includePadding=true)
Get a non-const view of the data as an MDArrayView.
Definition: Domi_MDVector.hpp:2455
int getUpperBndryPad(int axis) const
Get the size of the upper boundary padding along the given axis.
Definition: Domi_MDVector.hpp:1896
Invalid argument exception type.
Definition: Domi_Exceptions.hpp:53
int getUpperNeighbor(int axis) const
Get the rank of the upper neighbor.
Definition: Domi_MDVector.hpp:1753
MDMap Error exception type.
Definition: Domi_Exceptions.hpp:102
bool isContiguous() const
True if there are no stride gaps on any processor.
Definition: Domi_MDVector.hpp:1966