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 {};
175 template<
class Scalar,
176 class Node = DefaultNode::DefaultNodeType >
192 bool zeroOut =
true);
216 const dim_type leadingDim,
217 const dim_type trailingDim = 0,
218 bool zeroOut =
true);
253 Teuchos::DataAccess access = Teuchos::View);
268 MDVector(
const Teuchos::RCP<
const Teuchos::Comm< int > > teuchosComm,
269 Teuchos::ParameterList & plist);
283 MDVector(
const Teuchos::RCP< const MDComm > mdComm,
284 Teuchos::ParameterList & plist);
333 const Teuchos::ArrayView< Slice > & slices,
334 const Teuchos::ArrayView< int > & bndryPad =
335 Teuchos::ArrayView< int >());
355 inline const Teuchos::RCP< const MDMap< Node > >
372 inline Teuchos::RCP< const Teuchos::Comm< int > >
getTeuchosComm()
const;
460 inline dim_type
getGlobalDim(
int axis,
bool withBndryPad=
false)
const;
496 inline dim_type
getLocalDim(
int axis,
bool withCommPad=
false)
const;
670 Teuchos::RCP< Epetra_IntVector > getEpetraIntVectorView()
const;
687 Teuchos::RCP< Epetra_Vector > getEpetraVectorView()
const;
709 Teuchos::RCP< Epetra_MultiVector > getEpetraMultiVectorView()
const;
720 Teuchos::RCP< Epetra_IntVector > getEpetraIntVectorCopy()
const;
731 Teuchos::RCP< Epetra_Vector > getEpetraVectorCopy()
const;
749 Teuchos::RCP< Epetra_MultiVector > getEpetraMultiVectorCopy()
const;
765 template<
class LocalOrdinal,
766 class GlobalOrdinal = LocalOrdinal,
768 Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node2 > >
769 getTpetraVectorView()
const;
788 template <
class LocalOrdinal,
789 class GlobalOrdinal = LocalOrdinal,
791 Teuchos::RCP< Tpetra::MultiVector< Scalar,
795 getTpetraMultiVectorView()
const;
802 template<
class LocalOrdinal,
803 class GlobalOrdinal = LocalOrdinal,
805 Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node2 > >
806 getTpetraVectorCopy()
const;
820 template <
class LocalOrdinal,
821 class GlobalOrdinal = LocalOrdinal,
823 Teuchos::RCP< Tpetra::MultiVector< Scalar,
827 getTpetraMultiVectorCopy()
const;
892 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
norm1()
const;
896 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
norm2()
const;
900 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
normInf()
const;
906 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
930 describe(Teuchos::FancyOStream &out,
931 const Teuchos::EVerbosityLevel verbLevel =
932 Teuchos::Describable::verbLevel_default)
const;
947 bool includePadding =
true);
1047 bool includeBndryPad =
false)
const;
1056 void readBinary(
const std::string & filename,
1057 bool includeBndryPad =
false);
1066 Teuchos::RCP< const Teuchos::Comm< int > > _teuchosComm;
1070 Teuchos::RCP< const MDMap< Node > > _mdMap;
1091 Teuchos::Array< MPI_Request > _requests;
1104 Teuchos::RCP< MPI_Datatype > datatype;
1118 Teuchos::Array< Teuchos::Tuple< MessageInfo, 2 > > _sendMessages;
1123 Teuchos::Array< Teuchos::Tuple< MessageInfo, 2 > > _recvMessages;
1127 void initializeMessages();
1139 Teuchos::Array< dim_type > fileShape;
1140 Teuchos::Array< dim_type > bufferShape;
1141 Teuchos::Array< dim_type > dataShape;
1142 Teuchos::Array< dim_type > fileStart;
1143 Teuchos::Array< dim_type > dataStart;
1145 Teuchos::RCP< MPI_Datatype > filetype;
1146 Teuchos::RCP< MPI_Datatype > datatype;
1154 mutable Teuchos::RCP< FileInfo > _fileInfo;
1160 mutable Teuchos::RCP< FileInfo > _fileInfoWithBndry;
1168 Teuchos::RCP< FileInfo > & computeFileInfo(
bool includeBndryPad)
const;
1178 template<
class Scalar,
1183 _teuchosComm(mdMap->getTeuchosComm()),
1194 setObjectLabel(
"Domi::MDVector");
1197 int numDims = _mdMap->numDims();
1198 Teuchos::Array< dim_type > dims(
numDims);
1199 for (
int axis = 0; axis <
numDims; ++axis)
1200 dims[axis] = _mdMap->getLocalDim(axis,
true);
1203 _mdArrayRcp.
resize(dims);
1204 _mdArrayView = _mdArrayRcp();
1209 template<
class Scalar,
1213 const dim_type leadingDim,
1214 const dim_type trailingDim,
1216 _teuchosComm(mdMap->getTeuchosComm()),
1217 _mdMap(mdMap->getAugmentedMDMap(leadingDim, trailingDim)),
1227 setObjectLabel(
"Domi::MDVector");
1230 int numDims = _mdMap->numDims();
1231 Teuchos::Array< dim_type > dims(
numDims);
1232 for (
int axis = 0; axis <
numDims; ++axis)
1233 dims[axis] = _mdMap->getLocalDim(axis,
true);
1236 _mdArrayRcp.
resize(dims);
1237 _mdArrayView = _mdArrayRcp();
1242 template<
class Scalar,
1248 _mdArrayRcp(source),
1249 _mdArrayView(_mdArrayRcp()),
1257 setObjectLabel(
"Domi::MDVector");
1258 int numDims = _mdMap->numDims();
1259 TEUCHOS_TEST_FOR_EXCEPTION(
1262 "MDMap and source array do not have the same number of dimensions");
1264 for (
int axis = 0; axis <
numDims; ++axis)
1266 TEUCHOS_TEST_FOR_EXCEPTION(
1267 _mdMap->getLocalDim(axis) != _mdArrayRcp.
dimension(axis),
1269 "Axis " << axis <<
": MDMap dimension = " << _mdMap->getLocalDim(axis)
1270 <<
", MDArray dimension = " << _mdArrayRcp.
dimension(axis));
1276 template<
class Scalar,
1282 _mdArrayRcp(mdArrayRcp),
1283 _mdArrayView(_mdArrayRcp()),
1291 #ifdef DOMI_MDVECTOR_VERBOSE 1292 cout <<
"_mdArrayRcp = " << _mdArrayRcp << endl;
1293 cout <<
"_mdArrayView.getRawPtr() = " << _mdArrayView.
getRawPtr()
1294 <<
" (constructor)" << endl;
1295 cout <<
"_mdArrayView = " << _mdArrayView << endl;
1297 setObjectLabel(
"Domi::MDVector");
1298 int numDims = _mdMap->numDims();
1299 TEUCHOS_TEST_FOR_EXCEPTION(
1302 "MDMap and source array do not have the same number of dimensions");
1304 for (
int axis = 0; axis <
numDims; ++axis)
1306 TEUCHOS_TEST_FOR_EXCEPTION(
1307 _mdMap->getLocalDim(axis) != _mdArrayRcp.
dimension(axis),
1309 "Axis " << axis <<
": MDMap dimension = " << _mdMap->getLocalDim(axis)
1310 <<
", MDArray dimension = " << _mdArrayRcp.
dimension(axis));
1316 template<
class Scalar,
1320 Teuchos::DataAccess access) :
1321 _teuchosComm(source.getMDMap()->getTeuchosComm()),
1322 _mdMap(source.getMDMap()),
1323 _mdArrayRcp(source._mdArrayRcp),
1324 _mdArrayView(source._mdArrayView),
1332 setObjectLabel(
"Domi::MDVector");
1334 if (access == Teuchos::Copy)
1336 #ifdef DOMI_MDVECTOR_VERBOSE 1337 cout <<
"Inside MDVector copy constructor with copy access" << endl;
1340 int numDims = _mdMap->numDims();
1341 Teuchos::Array< dim_type > dims(
numDims);
1342 for (
int axis = 0; axis <
numDims; ++axis)
1343 dims[axis] = _mdMap->getLocalDim(axis,
true);
1347 _mdArrayView = _mdArrayRcp();
1352 const_iterator src = source.
getData().begin();
1353 for (iterator trg = _mdArrayView.
begin(); trg != _mdArrayView.
end(); ++trg)
1359 #ifdef DOMI_MDVECTOR_VERBOSE 1362 cout <<
"Inside MDVector copy constructor with view access" 1370 template<
class Scalar,
1373 MDVector(
const Teuchos::RCP<
const Teuchos::Comm< int > > teuchosComm,
1374 Teuchos::ParameterList & plist) :
1375 _teuchosComm(teuchosComm),
1386 setObjectLabel(
"Domi::MDVector");
1390 dim_type leadingDim = plist.get(
"leading dimension" , 0);
1391 dim_type trailingDim = plist.get(
"trailing dimension", 0);
1392 if (leadingDim + trailingDim > 0)
1398 _mdMap = Teuchos::rcp(myMdMap);
1401 int numDims = _mdMap->numDims();
1402 Teuchos::Array< dim_type > dims(
numDims);
1403 for (
int axis = 0; axis <
numDims; ++axis)
1404 dims[axis] = _mdMap->getLocalDim(axis,
true);
1407 _mdArrayRcp.
resize(dims);
1408 _mdArrayView = _mdArrayRcp();
1413 template<
class Scalar,
1417 Teuchos::ParameterList & plist) :
1418 _teuchosComm(mdComm->getTeuchosComm()),
1419 _mdMap(Teuchos::rcp(new
MDMap< Node >(mdComm, plist))),
1429 setObjectLabel(
"Domi::MDVector");
1433 dim_type leadingDim = plist.get(
"leading dimension" , 0);
1434 dim_type trailingDim = plist.get(
"trailing dimension", 0);
1435 if (leadingDim + trailingDim > 0)
1441 _mdMap = Teuchos::rcp(myMdMap);
1444 int numDims = _mdMap->numDims();
1445 Teuchos::Array< dim_type > dims(
numDims);
1446 for (
int axis = 0; axis <
numDims; ++axis)
1447 dims[axis] = _mdMap->getLocalDim(axis,
true);
1450 _mdArrayRcp.
resize(dims);
1451 _mdArrayView = _mdArrayRcp();
1456 template<
class Scalar,
1461 dim_type globalIndex) :
1462 _teuchosComm(parent._teuchosComm),
1464 _mdArrayRcp(parent._mdArrayRcp),
1465 _mdArrayView(parent._mdArrayView),
1473 setObjectLabel(
"Domi::MDVector");
1476 Teuchos::RCP< const MDMap< Node > > parentMdMap = parent.
getMDMap();
1484 if (_mdMap->onSubcommunicator())
1491 dim_type origin = parentMdMap->getGlobalRankBounds(axis,
false).start() -
1492 parentMdMap->getLowerPadSize(axis);
1497 dim_type localIndex = globalIndex - origin;
1501 _mdArrayView = newView;
1507 _mdArrayRcp.
clear();
1514 template<
class Scalar,
1519 const Slice & slice,
1523 _mdArrayRcp(parent._mdArrayRcp),
1524 _mdArrayView(parent._mdArrayView),
1532 #ifdef DOMI_MDVECTOR_VERBOSE 1533 cout <<
"slice axis " << axis << endl;
1534 cout <<
" _mdArrayRcp @ " << _mdArrayRcp.
getRawPtr() << endl;
1535 cout <<
" _mdArrayView @ " << _mdArrayView.
getRawPtr() << endl;
1538 setObjectLabel(
"Domi::MDVector");
1541 Teuchos::RCP< const MDMap< Node > > parentMdMap = parent.
getMDMap();
1548 _teuchosComm = _mdMap->getTeuchosComm();
1551 if (_mdMap->onSubcommunicator())
1554 Slice bounds = slice.
bounds(parentMdMap->getGlobalDim(axis,
true));
1561 dim_type origin = parentMdMap->getGlobalRankBounds(axis,
false).
start() -
1562 parentMdMap->getLowerPadSize(axis);
1569 dim_type start = std::max(0, bounds.
start() - origin - bndryPad);
1576 dim_type stop = std::min(bounds.
stop() - origin + bndryPad,
1577 parentMdMap->getLocalDim(axis,
true));
1581 _mdArrayView = newView;
1587 _mdArrayRcp.
clear();
1590 #ifdef DOMI_MDVECTOR_VERBOSE 1591 cout <<
" _mdArrayView @ " << _mdArrayView.
getRawPtr() << endl;
1592 cout <<
" offset = " << int(_mdArrayView.
getRawPtr() -
1599 template<
class Scalar,
1603 const Teuchos::ArrayView< Slice > & slices,
1604 const Teuchos::ArrayView< int > & bndryPad)
1606 setObjectLabel(
"Domi::MDVector");
1609 int numDims = parent.
numDims();
1612 TEUCHOS_TEST_FOR_EXCEPTION(
1613 (slices.size() != numDims),
1615 "number of slices = " << slices.size() <<
" != parent MDVector number of " 1616 "dimensions = " << numDims);
1620 for (
int axis = 0; axis < numDims; ++axis)
1622 int bndryPadding = (axis < bndryPad.size()) ? bndryPad[axis] : 0;
1627 tempMDVector1 = tempMDVector2;
1629 *
this = tempMDVector1;
1634 template<
class Scalar,
1640 _teuchosComm = source._teuchosComm;
1641 _mdMap = source._mdMap;
1642 _mdArrayRcp = source._mdArrayRcp;
1643 _mdArrayView = source._mdArrayView;
1644 _nextAxis = source._nextAxis;
1646 _requests = source._requests;
1648 _sendMessages = source._sendMessages;
1649 _recvMessages = source._recvMessages;
1655 template<
class Scalar,
1663 template<
class Scalar,
1665 const Teuchos::RCP< const MDMap< Node > >
1674 template<
class Scalar,
1680 return _mdMap->onSubcommunicator();
1685 template<
class Scalar,
1687 Teuchos::RCP< const Teuchos::Comm< int > >
1691 return _mdMap->getTeuchosComm();
1696 template<
class Scalar,
1702 return _mdMap->numDims();
1707 template<
class Scalar,
1713 return _mdMap->getCommDim(axis);
1718 template<
class Scalar,
1724 return _mdMap->isPeriodic(axis);
1729 template<
class Scalar,
1735 return _mdMap->getCommIndex(axis);
1740 template<
class Scalar,
1746 return _mdMap->getLowerNeighbor(axis);
1751 template<
class Scalar,
1757 return _mdMap->getUpperNeighbor(axis);
1762 template<
class Scalar,
1768 return _mdMap->getGlobalDim(axis, withBndryPad);
1773 template<
class Scalar,
1779 return _mdMap->getGlobalBounds(axis, withBndryPad);
1784 template<
class Scalar,
1790 return _mdMap->getGlobalRankBounds(axis, withBndryPad);
1795 template<
class Scalar,
1801 return _mdMap->getLocalDim(axis, withCommPad);
1806 template<
class Scalar,
1808 Teuchos::ArrayView< const Slice >
1812 return _mdMap->getLocalBounds();
1817 template<
class Scalar,
1823 return _mdMap->getLocalBounds(axis, withCommPad);
1828 template<
class Scalar,
1834 return _mdMap->getLocalInteriorBounds(axis);
1839 template<
class Scalar,
1845 return _mdMap->hasPadding();
1850 template<
class Scalar,
1856 return _mdMap->getLowerPadSize(axis);
1861 template<
class Scalar,
1867 return _mdMap->getUpperPadSize(axis);
1872 template<
class Scalar,
1878 return _mdMap->getCommPadSize(axis);
1883 template<
class Scalar,
1889 return _mdMap->getLowerBndryPad(axis);
1894 template<
class Scalar,
1900 return _mdMap->getUpperBndryPad(axis);
1905 template<
class Scalar,
1911 return _mdMap->getBndryPadSize(axis);
1916 template<
class Scalar,
1929 template<
class Scalar,
1942 template<
class Scalar,
1948 return _mdMap->isReplicatedBoundary(axis);
1953 template<
class Scalar,
1959 return _mdMap->getLayout();
1964 template<
class Scalar,
1970 return _mdMap->isContiguous();
1984 template<
class Scalar,
1986 Teuchos::RCP< Epetra_IntVector >
1991 TEUCHOS_TEST_FOR_EXCEPTION(
1992 typeid(Scalar) !=
typeid(
int),
1994 "MDVector is of scalar type '" <<
typeid(Scalar).name() <<
"', but " 1995 "Epetra_IntVector requires scalar type 'int'");
1998 TEUCHOS_TEST_FOR_EXCEPTION(
2001 "This MDVector's MDMap is non-contiguous. This can happen when you take " 2002 "a slice of a parent MDVector.");
2005 Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(
true);
2010 void * buffer = (
void*) _mdArrayView.getRawPtr();
2011 return Teuchos::rcp(
new Epetra_IntVector(View,
2025 template<
class Scalar,
2027 Teuchos::RCP< Epetra_Vector >
2028 MDVector< Scalar, Node >::
2029 getEpetraVectorView()
const 2032 TEUCHOS_TEST_FOR_EXCEPTION(
2033 typeid(Scalar) !=
typeid(
double),
2035 "MDVector is of scalar type '" <<
typeid(Scalar).name() <<
"', but " 2036 "Epetra_Vector requires scalar type 'double'");
2039 TEUCHOS_TEST_FOR_EXCEPTION(
2041 MDMapNoncontiguousError,
2042 "This MDVector's MDMap is non-contiguous. This can happen when you take " 2043 "a slice of a parent MDVector.");
2046 Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(
true);
2051 void * buffer = (
void*) _mdArrayView.getRawPtr();
2052 return Teuchos::rcp(
new Epetra_Vector(View,
2066 template<
class Scalar,
2068 Teuchos::RCP< Epetra_MultiVector >
2069 MDVector< Scalar, Node >::
2070 getEpetraMultiVectorView()
const 2073 TEUCHOS_TEST_FOR_EXCEPTION(
2074 typeid(Scalar) !=
typeid(
double),
2076 "MDVector is of scalar type '" <<
typeid(Scalar).name() <<
"', but " 2077 "Epetra_Vector requires scalar type 'double'");
2080 int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
2081 int padding = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
2082 int commDim = getCommDim(vectorAxis);
2083 int numVectors = getGlobalDim(vectorAxis);
2086 Teuchos::RCP< const MDMap< Node > > newMdMap;
2087 if (padding == 0 && commDim == 1)
2088 newMdMap = Teuchos::rcp(
new MDMap< Node >(*_mdMap, vectorAxis, 0));
2094 TEUCHOS_TEST_FOR_EXCEPTION(
2095 ! newMdMap->isContiguous(),
2096 MDMapNoncontiguousError,
2097 "This MDVector's MDMap is non-contiguous. This can happen when you take " 2098 "a slice of a parent MDVector.");
2103 size_type stride = newMdMap->getLocalDim(0,
true);
2104 for (
int axis = 0; axis < newMdMap->numDims(); ++axis)
2105 stride *= newMdMap->getLocalDim(axis,
true);
2106 TEUCHOS_TEST_FOR_EXCEPTION(
2107 stride*numVectors > Teuchos::OrdinalTraits<int>::max(),
2109 "Buffer size " << stride*numVectors <<
" is too large for Epetra int " 2111 int lda = (int)stride;
2114 Teuchos::RCP< const Epetra_Map > epetraMap = newMdMap->getEpetraMap(
true);
2119 void * buffer = (
void*) _mdArrayView.getRawPtr();
2120 return Teuchos::rcp(
new Epetra_MultiVector(View,
2129 template<
class Scalar,
2131 Teuchos::RCP< Epetra_IntVector >
2132 MDVector< Scalar, Node >::
2133 getEpetraIntVectorCopy()
const 2135 typedef typename MDArrayView< const Scalar >::iterator iterator;
2138 Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(
true);
2141 Teuchos::RCP< Epetra_IntVector > result =
2142 Teuchos::rcp(
new Epetra_IntVector(*epetraMap));
2147 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2148 result->operator[](ii++) = (int) *it;
2156 template<
class Scalar,
2158 Teuchos::RCP< Epetra_Vector >
2159 MDVector< Scalar, Node >::
2160 getEpetraVectorCopy()
const 2162 typedef typename MDArrayView< const Scalar >::iterator iterator;
2165 Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(
true);
2168 Teuchos::RCP< Epetra_Vector > result =
2169 Teuchos::rcp(
new Epetra_Vector(*epetraMap));
2174 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2175 result->operator[](ii++) = (double) *it;
2183 template<
class Scalar,
2185 Teuchos::RCP< Epetra_MultiVector >
2186 MDVector< Scalar, Node >::
2187 getEpetraMultiVectorCopy()
const 2189 typedef typename MDArrayView< Scalar >::iterator iterator;
2190 typedef typename MDArrayView< const Scalar >::iterator citerator;
2193 int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
2194 int padding = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
2195 int commDim = getCommDim(vectorAxis);
2196 int numVectors = getGlobalDim(vectorAxis);
2199 Teuchos::RCP< const MDMap< Node > > newMdMap;
2200 if (padding == 0 && commDim == 1)
2201 newMdMap = Teuchos::rcp(
new MDMap< Node >(*_mdMap, vectorAxis, 0));
2209 Teuchos::RCP< const Epetra_Map > epetraMap = newMdMap->getEpetraMap(
true);
2212 Teuchos::RCP< Epetra_MultiVector > result =
2213 Teuchos::rcp(
new Epetra_MultiVector(*epetraMap, numVectors));
2218 if (numVectors == 1)
2220 for (citerator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2221 result->operator[](0)[ii++] = (double) *it;
2225 for (
int iv = 0; iv < numVectors; ++iv)
2229 for (iterator it = subArray.begin(); it != subArray.end(); ++it)
2230 result->operator[](iv)[ii++] = (double) *it;
2244 template<
class Scalar,
class Node >
2245 template<
class LocalOrdinal,
2246 class GlobalOrdinal,
2248 Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node2 > >
2249 MDVector< Scalar, Node >::
2250 getTpetraVectorView()
const 2253 TEUCHOS_TEST_FOR_EXCEPTION(
2255 MDMapNoncontiguousError,
2256 "This MDVector's MDMap is non-contiguous. This can happen when you take " 2257 "a slice of a parent MDVector.");
2260 Teuchos::RCP<
const Tpetra::Map< LocalOrdinal,
2262 Node2 > > tpetraMap =
2263 _mdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal, Node2 >(
true);
2266 return Teuchos::rcp(
new Tpetra::Vector< Scalar,
2270 _mdArrayView.arrayView()));
2275 template<
class Scalar,
class Node >
2276 template<
class LocalOrdinal,
2277 class GlobalOrdinal,
2279 Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node2 > >
2280 MDVector< Scalar, Node >::
2281 getTpetraVectorCopy()
const 2283 typedef typename MDArrayView< const Scalar >::iterator iterator;
2286 Teuchos::RCP<
const Tpetra::Map< LocalOrdinal,
2288 Node2 > > tpetraMap =
2289 _mdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal, Node2 >(
true);
2292 Teuchos::RCP< Tpetra::Vector< Scalar,
2296 Teuchos::rcp(
new Tpetra::Vector< Scalar,
2299 Node2 >(tpetraMap) );
2303 Teuchos::ArrayRCP< Scalar > tpetraVectorArray = result->getDataNonConst();
2305 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2306 tpetraVectorArray[ii++] = *it;
2314 template<
class Scalar,
class Node >
2315 template <
class LocalOrdinal,
2316 class GlobalOrdinal,
2318 Teuchos::RCP< Tpetra::MultiVector< Scalar,
2322 MDVector< Scalar, Node >::
2323 getTpetraMultiVectorView()
const 2326 int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
2327 int padding = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
2328 int commDim = getCommDim(vectorAxis);
2329 size_t numVectors = getGlobalDim(vectorAxis);
2332 Teuchos::RCP< const MDMap< Node > > newMdMap;
2333 if (padding == 0 && commDim == 1)
2334 newMdMap = Teuchos::rcp(
new MDMap< Node >(*_mdMap, vectorAxis, 0));
2340 TEUCHOS_TEST_FOR_EXCEPTION(
2341 ! newMdMap->isContiguous(),
2342 MDMapNoncontiguousError,
2343 "This MDVector's MDMap is non-contiguous. This can happen when you take " 2344 "a slice of a parent MDVector.");
2349 size_type stride = newMdMap->getLocalDim(0,
true);
2350 for (
int axis = 0; axis < newMdMap->numDims(); ++axis)
2351 stride *= newMdMap->getLocalDim(axis,
true);
2352 TEUCHOS_TEST_FOR_EXCEPTION(
2353 stride*numVectors > Teuchos::OrdinalTraits<GlobalOrdinal>::max(),
2355 "Buffer size " << stride*numVectors <<
" is too large for Tpetra " 2356 "GlobalOrdinal = " <<
typeid(GlobalOrdinal).name() );
2357 size_t lda = (size_t)stride;
2360 Teuchos::RCP<
const Tpetra::Map< LocalOrdinal,
2362 Node2> > tpetraMap =
2363 newMdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal, Node2 >(
true);
2366 return Teuchos::rcp(
new Tpetra::MultiVector< Scalar,
2370 _mdArrayView.arrayView(),
2377 template<
class Scalar,
class Node >
2378 template <
class LocalOrdinal,
2379 class GlobalOrdinal,
2381 Teuchos::RCP< Tpetra::MultiVector< Scalar,
2385 MDVector< Scalar, Node >::
2386 getTpetraMultiVectorCopy()
const 2388 typedef typename MDArrayView< Scalar >::iterator iterator;
2389 typedef typename MDArrayView< const Scalar >::iterator citerator;
2392 int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
2393 int padding = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
2394 int commDim = getCommDim(vectorAxis);
2395 int numVectors = getGlobalDim(vectorAxis);
2398 Teuchos::RCP< const MDMap< Node > > newMdMap;
2399 if (padding == 0 && commDim == 1)
2400 newMdMap = Teuchos::rcp(
new MDMap< Node >(*_mdMap, vectorAxis, 0));
2408 Teuchos::RCP<
const Tpetra::Map< LocalOrdinal,
2410 Node2 > > tpetraMap =
2411 newMdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal, Node2 >(
true);
2414 Teuchos::RCP< Tpetra::MultiVector< Scalar,
2418 Teuchos::rcp(
new Tpetra::MultiVector< Scalar,
2421 Node2 >(tpetraMap, numVectors));
2426 if (numVectors == 1)
2428 Teuchos::ArrayRCP< Scalar > tpetraVectorArray = result->getDataNonConst(0);
2429 for (citerator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2430 tpetraVectorArray[ii++] = (
double) *it;
2434 for (
int iv = 0; iv < numVectors; ++iv)
2437 Teuchos::ArrayRCP< Scalar > tpetraVectorArray =
2438 result->getDataNonConst(iv);
2440 for (iterator it = subArray.begin(); it != subArray.end(); ++it)
2441 tpetraVectorArray[ii++] = (
double) *it;
2453 template<
class Scalar,
2459 #ifdef DOMI_MDVECTOR_VERBOSE 2460 cout <<
"_mdArrayView.getRawPtr() = " << _mdArrayView.
getRawPtr()
2462 cout <<
"_mdArrayView = " << _mdArrayView << endl;
2465 return _mdArrayView;
2467 for (
int axis = 0; axis < numDims(); ++axis)
2469 int lo = getLowerPadSize(axis);
2470 int hi = getLocalDim(axis,
true) - getUpperPadSize(axis);
2478 template<
class Scalar,
2485 return _mdArrayView.getConst();
2487 for (
int axis = 0; axis < numDims(); ++axis)
2489 int lo = getLowerPadSize(axis);
2490 int hi = getLocalDim(axis,
true) - getUpperPadSize(axis);
2498 template<
class Scalar,
2506 Slice(getLowerPadSize(axis)));
2507 return newArrayView;
2512 template<
class Scalar,
2520 Slice(getLowerPadSize(axis)));
2521 return newArrayView;
2526 template<
class Scalar,
2532 dim_type n = getLocalDim(axis,
true);
2533 int pad = getUpperPadSize(axis);
2535 if (pad) slice =
Slice(n-pad,n);
2536 else slice =
Slice(n-1,n-1);
2540 return newArrayView;
2545 template<
class Scalar,
2551 dim_type n = getLocalDim(axis,
true);
2552 int pad = getUpperPadSize(axis);
2554 if (pad) slice =
Slice(n-pad,n);
2555 else slice =
Slice(n-1,n-1);
2559 return newArrayView;
2564 template<
class Scalar,
2572 TEUCHOS_TEST_FOR_EXCEPTION(
2573 ! _mdMap->isCompatible(*(a._mdMap)),
2575 "MDMap of calling MDVector and argument 'a' are incompatible");
2578 Scalar local_dot = 0;
2579 iterator a_it = aView.begin();
2580 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end();
2582 local_dot += *it * *a_it;
2583 Scalar global_dot = 0;
2584 Teuchos::reduceAll(*_teuchosComm,
2585 Teuchos::REDUCE_SUM,
2594 template<
class Scalar,
2596 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
2600 typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
2603 mag local_norm1 = 0;
2604 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2605 local_norm1 += std::abs(*it);
2606 mag global_norm1 = 0;
2607 Teuchos::reduceAll(*_teuchosComm,
2608 Teuchos::REDUCE_SUM,
2612 return global_norm1;
2617 template<
class Scalar,
2619 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
2623 typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
2624 mag norm2 = dot(*
this);
2625 return Teuchos::ScalarTraits<mag>::squareroot(norm2);
2630 template<
class Scalar,
2632 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
2636 typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
2639 mag local_normInf = 0;
2640 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2641 local_normInf = std::max(local_normInf, std::abs(*it));
2642 mag global_normInf = 0;
2643 Teuchos::reduceAll(*_teuchosComm,
2644 Teuchos::REDUCE_MAX,
2648 return global_normInf;
2653 template<
class Scalar,
2655 typename Teuchos::ScalarTraits< Scalar >::magnitudeType
2659 typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
2662 TEUCHOS_TEST_FOR_EXCEPTION(
2663 ! _mdMap->isCompatible(*(weights._mdMap)),
2665 "MDMap of calling MDVector and argument 'weights' are incompatible");
2668 mag local_wNorm = 0;
2669 iterator w_it = wView.begin();
2670 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end();
2672 local_wNorm += *it * *it * *w_it;
2673 mag global_wNorm = 0;
2674 Teuchos::reduceAll(*_teuchosComm,
2675 Teuchos::REDUCE_SUM,
2679 Teuchos::Array< dim_type > dimensions(numDims());
2680 for (
int i = 0; i < numDims(); ++i)
2681 dimensions[i] = _mdMap->getGlobalDim(i);
2682 size_type n = computeSize(dimensions);
2683 if (n == 0)
return 0;
2684 return Teuchos::ScalarTraits<mag>::squareroot(global_wNorm / n);
2689 template<
class Scalar,
2695 typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
2699 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2702 Teuchos::reduceAll(*_teuchosComm,
2703 Teuchos::REDUCE_SUM,
2707 Teuchos::Array< dim_type > dimensions(numDims());
2708 for (
int i = 0; i < numDims(); ++i)
2709 dimensions[i] = _mdMap->getGlobalDim(i);
2710 size_type n = computeSize(dimensions);
2711 if (n == 0)
return 0;
2712 return global_sum / n;
2717 template<
class Scalar,
2723 using Teuchos::TypeNameTraits;
2725 Teuchos::Array< dim_type > dims(numDims());
2726 for (
int axis = 0; axis < numDims(); ++axis)
2727 dims[axis] = getGlobalDim(axis,
true);
2729 std::ostringstream oss;
2730 oss <<
"\"Domi::MDVector\": {" 2731 <<
"Template parameters: {" 2732 <<
"Scalar: " << TypeNameTraits<Scalar>::name()
2733 <<
", Node: " << TypeNameTraits< Node >::name()
2735 if (this->getObjectLabel() !=
"")
2736 oss <<
", Label: \"" << this->getObjectLabel () <<
"\", ";
2737 oss <<
"Global dimensions: " << dims <<
" }";
2743 template<
class Scalar,
2748 const Teuchos::EVerbosityLevel verbLevel)
const 2751 using Teuchos::Comm;
2753 using Teuchos::TypeNameTraits;
2754 using Teuchos::VERB_DEFAULT;
2755 using Teuchos::VERB_NONE;
2756 using Teuchos::VERB_LOW;
2757 using Teuchos::VERB_MEDIUM;
2758 using Teuchos::VERB_HIGH;
2759 using Teuchos::VERB_EXTREME;
2761 const Teuchos::EVerbosityLevel vl =
2762 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
2765 Teuchos::RCP< const Teuchos::Comm< int > > comm = mdMap.
getTeuchosComm();
2766 const int myImageID = comm->getRank();
2767 const int numImages = comm->getSize();
2768 Teuchos::OSTab tab0(out);
2770 if (vl != VERB_NONE)
2774 out <<
"\"Domi::MDVector\":" << endl;
2776 Teuchos::OSTab tab1(out);
2779 out <<
"Template parameters:";
2781 Teuchos::OSTab tab2(out);
2782 out <<
"Scalar: " << TypeNameTraits<Scalar>::name() << endl
2783 <<
"Node: " << TypeNameTraits< Node >::name() << endl;
2786 if (this->getObjectLabel() !=
"")
2788 out <<
"Label: \"" << getObjectLabel() <<
"\"" << endl;
2790 Teuchos::Array< dim_type > globalDims(numDims());
2791 for (
int axis = 0; axis < numDims(); ++axis)
2792 globalDims[axis] = getGlobalDim(axis,
true);
2793 out <<
"Global dimensions: " << globalDims << endl;
2795 for (
int imageCtr = 0; imageCtr < numImages; ++imageCtr)
2797 if (myImageID == imageCtr)
2802 out <<
"Process: " << myImageID << endl;
2803 Teuchos::OSTab tab2(out);
2804 Teuchos::Array< dim_type > localDims(numDims());
2805 for (
int axis = 0; axis < numDims(); ++axis)
2806 localDims[axis] = getLocalDim(axis,
true);
2807 out <<
"Local dimensions: " << localDims << endl;
2820 template<
class Scalar,
2825 bool includePadding)
2830 for (iterator it = newArray.
begin(); it != newArray.
end(); ++it)
2836 template<
class Scalar,
2844 Teuchos::ScalarTraits< Scalar >::seedrandom(time(NULL));
2845 for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
2846 *it = Teuchos::ScalarTraits< Scalar >::random();
2851 template<
class Scalar,
2857 for (
int axis = 0; axis < numDims(); ++axis)
2859 updateCommPad(axis);
2865 template<
class Scalar,
2871 startUpdateCommPad(axis);
2872 endUpdateCommPad(axis);
2877 template<
class Scalar,
2887 if (_sendMessages.empty()) initializeMessages();
2890 int rank = _teuchosComm->getRank();
2891 int numProc = _teuchosComm->getSize();
2897 Teuchos::RCP< const Teuchos::MpiComm< int > > mpiComm =
2898 Teuchos::rcp_dynamic_cast<
const Teuchos::MpiComm< int > >(_teuchosComm);
2899 const Teuchos::OpaqueWrapper< MPI_Comm > & communicator =
2900 *(mpiComm->getRawMpiComm());
2903 MPI_Request request;
2904 for (
int boundary = 0; boundary < 2; ++boundary)
2906 MessageInfo message = _sendMessages[axis][boundary];
2907 if (message.proc >= 0)
2909 tag = 2 * (rank * numProc + message.proc) + boundary;
2911 #ifdef DOMI_MDVECTOR_OUTPUT_UPDATECOMMPAD 2912 cout << rank <<
": post send for axis " << axis <<
", boundary " 2913 << boundary <<
", buffer = " << message.buffer <<
", proc = " 2914 << message.proc <<
", tag = " << tag << endl;
2917 if (MPI_Isend(message.buffer,
2919 *(message.datatype),
2924 throw std::runtime_error(
"Domi::MDVector: Error in MPI_Isend");
2925 _requests.push_back(request);
2930 for (
int boundary = 0; boundary < 2; ++boundary)
2932 MessageInfo message = _recvMessages[axis][boundary];
2933 if (message.proc >= 0)
2935 tag = 2 * (message.proc * numProc + rank) + (1-boundary);
2937 #ifdef DOMI_MDVECTOR_OUTPUT_UPDATECOMMPAD 2938 cout << rank <<
": post recv for axis " << axis <<
", boundary " 2939 << boundary <<
", buffer = " << message.buffer <<
", proc = " 2940 << message.proc <<
", tag = " << tag << endl;
2943 if (MPI_Irecv(message.buffer,
2945 *(message.datatype),
2950 throw std::runtime_error(
"Domi::MDVector: Error in MPI_Irecv");
2951 _requests.push_back(request);
2958 if (isPeriodic(axis))
2960 for (
int sendBndry = 0; sendBndry < 2; ++sendBndry)
2962 int recvBndry = 1 - sendBndry;
2972 for ( ; it_recv != recvView.
end(); ++it_recv, ++it_send)
2973 *it_recv = *it_send;
2981 template<
class Scalar,
2988 if (_requests.size() > 0)
2990 Teuchos::Array< MPI_Status > status(_requests.size());
2991 if (MPI_Waitall(_requests.size(),
2994 throw std::runtime_error(
"Domi::MDVector: Error in MPI_Waitall");
3002 template<
class Scalar,
3011 int newAxis = _nextAxis + 1;
3012 if (newAxis >= result.
numDims()) newAxis = 0;
3013 result._nextAxis = newAxis;
3019 template<
class Scalar,
3028 int newAxis = _nextAxis + 1;
3029 if (newAxis >= result.
numDims()) newAxis = 0;
3030 result._nextAxis = newAxis;
3036 template<
class Scalar,
3044 int ndims = numDims();
3045 Teuchos::Array<int> sizes(ndims);
3046 Teuchos::Array<int> subsizes(ndims);
3047 Teuchos::Array<int> starts(ndims);
3048 MessageInfo messageInfo;
3050 _sendMessages.resize(ndims);
3051 _recvMessages.resize(ndims);
3053 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3054 std::stringstream msg;
3055 int rank = _teuchosComm->getRank();
3059 int order = mpiOrder(getLayout());
3060 MPI_Datatype datatype = mpiType< Scalar >();
3064 for (
int msgAxis = 0; msgAxis < ndims; ++msgAxis)
3067 for (
int axis = 0; axis < ndims; ++axis)
3069 sizes[axis] = _mdArrayRcp.dimension(axis);
3070 subsizes[axis] = _mdArrayView.dimension(axis);
3078 int proc = getLowerNeighbor(msgAxis);
3080 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3081 msg << endl <<
"P" << rank <<
": axis " << msgAxis <<
", lower neighbor = " 3086 subsizes[msgAxis] = getLowerPadSize(msgAxis);
3088 if (subsizes[msgAxis] == 0) proc = -1;
3090 messageInfo.buffer = (
void*) getData().getRawPtr();
3091 messageInfo.proc = proc;
3092 messageInfo.axis = msgAxis;
3102 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3103 msg << endl <<
"P" << rank <<
": axis " << msgAxis
3104 <<
", Lower receive message:" << endl <<
" ndims = " << ndims
3105 << endl <<
" sizes = " << sizes << endl <<
" subsizes = " 3106 << subsizes << endl <<
" starts = " << starts << endl
3107 <<
" order = " << order << endl;
3109 Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(
new MPI_Datatype);
3110 MPI_Type_create_subarray(ndims,
3117 MPI_Type_commit(commPad.get());
3118 messageInfo.datatype = commPad;
3120 messageInfo.dataview = _mdArrayView;
3121 for (
int axis = 0; axis < numDims(); ++axis)
3123 Slice slice(starts[axis], starts[axis] + subsizes[axis]);
3131 _recvMessages[msgAxis][0] = messageInfo;
3137 starts[msgAxis] = getLowerPadSize(msgAxis);
3138 if (isReplicatedBoundary(msgAxis) && getCommIndex(msgAxis) == 0)
3139 starts[msgAxis] += 1;
3145 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3146 msg << endl <<
"P" << rank <<
": axis " << msgAxis
3147 <<
", Lower send message:" << endl <<
" ndims = " << ndims
3148 << endl <<
" sizes = " << sizes << endl <<
" subsizes = " 3149 << subsizes << endl <<
" starts = " << starts << endl
3150 <<
" order = " << order << endl;
3152 Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(
new MPI_Datatype);
3153 MPI_Type_create_subarray(ndims,
3160 MPI_Type_commit(commPad.get());
3161 messageInfo.datatype = commPad;
3163 messageInfo.dataview = _mdArrayView;
3164 for (
int axis = 0; axis < numDims(); ++axis)
3166 Slice slice(starts[axis], starts[axis] + subsizes[axis]);
3174 _sendMessages[msgAxis][0] = messageInfo;
3180 proc = getUpperNeighbor(msgAxis);
3182 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3183 msg << endl <<
"P" << rank <<
": axis " << msgAxis <<
", upper neighbor = " 3187 subsizes[msgAxis] = getUpperPadSize(msgAxis);
3188 starts[msgAxis] = _mdArrayView.dimension(msgAxis) -
3189 getUpperPadSize(msgAxis);
3190 if (subsizes[msgAxis] == 0) proc = -1;
3191 messageInfo.proc = proc;
3200 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3201 msg << endl <<
"P" << rank <<
": axis " << msgAxis
3202 <<
", Upper receive message:" << endl <<
" ndims = " << ndims
3203 << endl <<
" sizes = " << sizes << endl <<
" subsizes = " 3204 << subsizes << endl <<
" starts = " << starts << endl
3205 <<
" order = " << order << endl;
3207 Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(
new MPI_Datatype);
3208 MPI_Type_create_subarray(ndims,
3215 MPI_Type_commit(commPad.get());
3216 messageInfo.datatype = commPad;
3218 messageInfo.dataview = _mdArrayView;
3219 for (
int axis = 0; axis < numDims(); ++axis)
3221 Slice slice(starts[axis], starts[axis] + subsizes[axis]);
3228 _recvMessages[msgAxis][1] = messageInfo;
3234 starts[msgAxis] -= getUpperPadSize(msgAxis);
3235 if (isReplicatedBoundary(msgAxis) &&
3236 getCommIndex(msgAxis) == getCommDim(msgAxis)-1)
3237 starts[msgAxis] -= 1;
3243 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3244 msg << endl <<
"P" << rank <<
": axis " << msgAxis
3245 <<
", Upper send message:" << endl <<
" ndims = " << ndims
3246 << endl <<
" sizes = " << sizes << endl <<
" subsizes = " 3247 << subsizes << endl <<
" starts = " << starts << endl
3248 <<
" order = " << order << endl;
3250 Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(
new MPI_Datatype);
3251 MPI_Type_create_subarray(ndims,
3258 MPI_Type_commit(commPad.get());
3259 messageInfo.datatype = commPad;
3261 messageInfo.dataview = _mdArrayView;
3262 for (
int axis = 0; axis < numDims(); ++axis)
3264 Slice slice(starts[axis], starts[axis] + subsizes[axis]);
3272 _sendMessages[msgAxis][1] = messageInfo;
3275 #ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE 3276 for (
int proc = 0; proc < _teuchosComm->getSize(); ++proc)
3283 _teuchosComm->barrier();
3291 template<
class Scalar,
3296 bool includeBndryPad)
const 3303 int pid = _teuchosComm->getRank();
3306 datafile = fopen(filename.c_str(),
"w");
3309 _teuchosComm->barrier();
3312 const Scalar * buffer = getData(
true).getRawPtr();
3316 Teuchos::RCP< FileInfo > & fileInfo = computeFileInfo(includeBndryPad);
3325 Teuchos::RCP< const Teuchos::MpiComm< int > > mpiComm =
3326 Teuchos::rcp_dynamic_cast<
const Teuchos::MpiComm< int > >(_teuchosComm);
3327 const Teuchos::OpaqueWrapper< MPI_Comm > & communicator =
3328 *(mpiComm->getRawMpiComm());
3331 int access = MPI_MODE_WRONLY | MPI_MODE_CREATE;
3336 char * cstr =
new char[filename.size()+1];
3337 std::strcpy(cstr, filename.c_str());
3342 char datarep[7] =
"native";
3343 MPI_File_open(communicator(), cstr, access, MPI_INFO_NULL, &mpiFile);
3344 MPI_File_set_view(mpiFile, 0, mpiType< Scalar >(),
3345 *(fileInfo->filetype), datarep, MPI_INFO_NULL);
3346 MPI_File_write_all(mpiFile, (
void*)buffer, 1, *(fileInfo->datatype),
3348 MPI_File_close(&mpiFile);
3360 datafile = fopen(filename.c_str(),
"w");
3367 for (const_iterator it = mdArrayView.begin(); it != mdArrayView.end(); ++it)
3369 fwrite((
const void *) &(*it),
sizeof(Scalar), 1, datafile);
3381 template<
class Scalar,
3386 bool includeBndryPad)
3389 const Scalar * buffer = getDataNonConst(
true).getRawPtr();
3393 Teuchos::RCP< FileInfo > & fileInfo = computeFileInfo(includeBndryPad);
3402 Teuchos::RCP< const Teuchos::MpiComm< int > > mpiComm =
3403 Teuchos::rcp_dynamic_cast<
const Teuchos::MpiComm< int > >(_teuchosComm);
3404 const Teuchos::OpaqueWrapper< MPI_Comm > & communicator =
3405 *(mpiComm->getRawMpiComm());
3408 int access = MPI_MODE_RDONLY;
3413 char * cstr =
new char[filename.size()+1];
3414 std::strcpy(cstr, filename.c_str());
3419 char datarep[7] =
"native";
3420 MPI_File_open(communicator(), cstr, access, MPI_INFO_NULL, &mpiFile);
3421 MPI_File_set_view(mpiFile, 0, mpiType< Scalar >(),
3422 *(fileInfo->filetype), datarep, MPI_INFO_NULL);
3423 MPI_File_read_all(mpiFile, (
void*)buffer, 1, *(fileInfo->datatype),
3425 MPI_File_close(&mpiFile);
3434 int ndims = _mdMap->numDims();
3438 datafile = fopen(filename.c_str(),
"r");
3445 Teuchos::Array< Ordinal > index(3);
3446 for (
int axis = 0; axis < ndims; ++axis)
3447 index[axis] = fileInfo->dataStart[axis];
3451 for (iterator it = mdArrayView.begin(); it != mdArrayView.end(); ++it)
3453 fread(&(*it),
sizeof(Scalar), 1, datafile);
3465 template<
class Scalar,
3467 Teuchos::RCP< typename MDVector< Scalar, Node >::FileInfo > &
3473 Teuchos::RCP< MDVector< Scalar, Node >::FileInfo > & fileInfo =
3474 includeBndryPad ? _fileInfoWithBndry : _fileInfo;
3477 if (!fileInfo.is_null())
return fileInfo;
3480 int ndims = _mdMap->
numDims();
3482 fileInfo->fileShape.resize(ndims);
3483 fileInfo->bufferShape.resize(ndims);
3484 fileInfo->dataShape.resize(ndims);
3485 fileInfo->fileStart.resize(ndims);
3486 fileInfo->dataStart.resize(ndims);
3489 for (
int axis = 0; axis < ndims; ++axis)
3493 fileInfo->fileShape[axis] = getGlobalDim(axis,includeBndryPad);
3494 fileInfo->bufferShape[axis] = getLocalDim(axis,
true );
3495 fileInfo->dataShape[axis] = getLocalDim(axis,
false);
3496 fileInfo->fileStart[axis] = getGlobalRankBounds(axis,includeBndryPad).start();
3497 fileInfo->dataStart[axis] = getLocalBounds(axis).start();
3499 if (includeBndryPad)
3501 int commIndex = _mdMap->getCommIndex(axis);
3504 int pad = getLowerBndryPad(axis);
3505 fileInfo->dataShape[axis] += pad;
3506 fileInfo->dataStart[axis] -= pad;
3508 if (commIndex == _mdMap->getCommDim(axis)-1)
3510 fileInfo->dataShape[axis] += getUpperBndryPad(axis);
3515 #ifdef DOMI_MDVECTOR_DEBUG_IO 3516 cout << pid <<
": fileShape = " << fileInfo->fileShape() << endl;
3517 cout << pid <<
": bufferShape = " << fileInfo->bufferShape() << endl;
3518 cout << pid <<
": dataShape = " << fileInfo->dataShape() << endl;
3519 cout << pid <<
": fileStart = " << fileInfo->fileStart() << endl;
3520 cout << pid <<
": dataStart = " << fileInfo->dataStart() << endl;
3524 int mpi_order = getLayout() == C_ORDER ? MPI_ORDER_C : MPI_ORDER_FORTRAN;
3526 fileInfo->filetype = Teuchos::rcp(
new MPI_Datatype);
3527 MPI_Type_create_subarray(ndims,
3528 fileInfo->fileShape.getRawPtr(),
3529 fileInfo->dataShape.getRawPtr(),
3530 fileInfo->fileStart.getRawPtr(),
3532 mpiType< Scalar >(),
3533 fileInfo->filetype.get());
3534 MPI_Type_commit(fileInfo->filetype.get());
3537 fileInfo->datatype = Teuchos::rcp(
new MPI_Datatype);
3538 MPI_Type_create_subarray(ndims,
3539 fileInfo->bufferShape.getRawPtr(),
3540 fileInfo->dataShape.getRawPtr(),
3541 fileInfo->dataStart.getRawPtr(),
3543 mpiType< Scalar >(),
3544 fileInfo->datatype.get());
3545 MPI_Type_commit(fileInfo->datatype.get());
3546 #endif // DGM_PARALLEL void writeBinary(const std::string &filename, bool includeBndryPad=false) const
Write the MDVector to a binary file.
Definition: Domi_MDVector.hpp:3295
int getCommPadSize(int axis) const
Get the communication padding size along the given axis.
Definition: Domi_MDVector.hpp:1876
Scalar meanValue() const
Compute the mean (average) value of this MDVector.
Definition: Domi_MDVector.hpp:2693
MDVector< Scalar, Node > operator[](dim_type index) const
Sub-vector access operator.
Definition: Domi_MDVector.hpp:3006
virtual std::string description() const
A simple one-line description of this MDVector.
Definition: Domi_MDVector.hpp:2721
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:1865
MDVector(const Teuchos::RCP< const MDMap< Node > > &mdMap, bool zeroOut=true)
Main constructor.
Definition: Domi_MDVector.hpp:1181
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:2747
int getLowerPadSize(int axis) const
Get the size of the lower padding along the given axis.
Definition: Domi_MDVector.hpp:1854
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
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:2502
Teuchos::ArrayView< const Slice > getLocalBounds() const
Get the local loop bounds along every axis.
Definition: Domi_MDVector.hpp:1810
void setLowerPad(int axis, const Scalar value)
Assign all elements of the lower pad a constant value.
Definition: Domi_MDVector.hpp:1920
Teuchos::ScalarTraits< Scalar >::magnitudeType norm1() const
Compute the 1-norm of this MDVector.
Definition: Domi_MDVector.hpp:2598
void setUpperPad(int axis, const Scalar value)
Assign all elements of the upper pad a constant value.
Definition: Domi_MDVector.hpp:1933
dim_type getLocalDim(int axis, bool withCommPad=false) const
Get the local dimension along the specified axis.
Definition: Domi_MDVector.hpp:1799
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:2516
bool onSubcommunicator() const
Query whether this processor is on the sub-communicator.
Definition: Domi_MDVector.hpp:1678
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:2568
void startUpdateCommPad(int axis)
Start an asyncronous update of the communication padding.
Definition: Domi_MDVector.hpp:2881
MDArrayView< const Scalar > getData(bool includePadding=true) const
Get a const view of the data as an MDArrayView.
Definition: Domi_MDVector.hpp:2482
bool isReplicatedBoundary(int axis) const
Return whether the given axis supports a replicated boundary.
Definition: Domi_MDVector.hpp:1946
int numDims() const
Get the number of dimensions.
Definition: Domi_MDVector.hpp:1700
Layout getLayout() const
Get the storage order.
Definition: Domi_MDVector.hpp:1957
dim_type getGlobalDim(int axis, bool withBndryPad=false) const
Get the global dimension along the specified axis.
Definition: Domi_MDVector.hpp:1766
int getCommDim(int axis) const
Get the communicator size along the given axis.
Definition: Domi_MDVector.hpp:1711
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:3385
int getCommIndex(int axis) const
Get the axis rank of this processor.
Definition: Domi_MDVector.hpp:1733
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:2530
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:1843
Teuchos::ScalarTraits< Scalar >::magnitudeType norm2() const
Compute the 2-norm of this MDVector.
Definition: Domi_MDVector.hpp:2621
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:2634
Teuchos::RCP< const Teuchos::Comm< int > > getTeuchosComm() const
Get the Teuchos communicator.
Definition: Domi_MDMap.hpp:2068
bool isPeriodic(int axis) const
Return the periodic flag for the given axis.
Definition: Domi_MDVector.hpp:1722
void endUpdateCommPad(int axis)
Complete an asyncronous update of the communication padding.
Definition: Domi_MDVector.hpp:2985
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:2549
void randomize()
Set all values in the multivector to pseudorandom numbers.
Definition: Domi_MDVector.hpp:2840
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:1638
virtual Slice bounds(dim_type size) const
Return a Slice with concrete start and stop values.
virtual ~MDVector()
Destructor.
Definition: Domi_MDVector.hpp:1657
Slice getLocalInteriorBounds(int axis) const
Get the local interior looping bounds along the specified axis.
Definition: Domi_MDVector.hpp:1832
void putScalar(const Scalar &value, bool includePadding=true)
Set all values in the multivector with the given value.
Definition: Domi_MDVector.hpp:2824
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:1909
Teuchos::ScalarTraits< Scalar >::magnitudeType normWeighted(const MDVector< Scalar, Node > &weights) const
Compute the weighted norm of this.
Definition: Domi_MDVector.hpp:2657
Slice getGlobalRankBounds(int axis, bool withBndryPad=false) const
Get the global loop bounds on this processor along the specified axis.
Definition: Domi_MDVector.hpp:1788
Multi-dimensional distributed vector.
Definition: Domi_MDVector.hpp:177
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:1777
Teuchos::RCP< const Teuchos::Comm< int > > getTeuchosComm() const
Get the Teuchos communicator.
Definition: Domi_MDVector.hpp:1689
void updateCommPad()
Sum values of a locally replicated multivector across all processes.
Definition: Domi_MDVector.hpp:2855
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:1744
int getLowerBndryPad(int axis) const
Get the size of the lower boundary padding along the given axis.
Definition: Domi_MDVector.hpp:1887
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:1667
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:2457
int getUpperBndryPad(int axis) const
Get the size of the upper boundary padding along the given axis.
Definition: Domi_MDVector.hpp:1898
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:1755
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:1968