portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-09-20 09:15:38 +01:00
Code Releases Activity

Finalising traits

Browse Source
This commit is contained in:
Michael Marshall 2019-02-14 19:05:35 +00:00
parent 59c8cc1588
commit bee24655cd
5 changed files with 267 additions and 205 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

249
Grid/serialisation/BaseIO.h
View File

@ -36,86 +36,193 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
#include <Grid/Eigen/unsupported/CXX11/Tensor> #include <Grid/Eigen/unsupported/CXX11/Tensor>
namespace Grid { namespace Grid {
// Abstract writer/reader classes //////////////////////////////////////////// namespace EigenIO {
// static polymorphism implemented using CRTP idiom template<typename T> struct is_complex : public std::false_type {};
class Serializable; template<typename T> struct is_complex<std::complex<T>>
: std::integral_constant<bool, std::is_arithmetic<T>::value> {};
// which types are supported scalar types for Eigen::Tensor // Eigen tensors can be composed of arithmetic scalar and complex types
template<typename T> struct is_eigen_tensor_scalar : std::integral_constant<bool, template<typename T> struct is_scalar : std::integral_constant<bool,
std::is_arithmetic<T>::value || Grid::is_complex<T>::value> {}; std::is_arithmetic<T>::value || is_complex<T>::value> {};
// Eigen tensors can also be composed of a limited number of containers
// NB: grid tensors (iScalar, iVector, iMatrix) have stricter limits on complex types
template <typename T> struct is_container : public std::false_type {};
template <typename T> struct is_container<iScalar<T>> : public std::true_type {};
template <typename T, int N> struct is_container<iVector<T, N>> : public std::true_type {};
template <typename T, int N> struct is_container<iMatrix<T, N>> : public std::true_type {};
template <typename T, std::size_t N> struct is_container<std::array<T, N>> : public std::true_type {};
// Is this an Eigen tensor
template<typename T> struct is_tensor : std::integral_constant<bool,
std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value> {};
// Is this an Eigen tensor of a supported scalar
template<typename T> struct is_tensor_of_scalar : std::integral_constant<bool,
is_tensor<T>::value && is_scalar<typename T::Scalar>::value> {};
// Is this an Eigen tensor of a supported container
template<typename T> struct is_tensor_of_container : std::integral_constant<bool,
is_tensor<T>::value && is_container<typename T::Scalar>::value> {};
// These traits describe the Eigen tensor scalar objects supported for IO
// Mostly intended for grid tensors, i.e. compositions of iScalar, iVector and iMatrix
// but also supports fixed size arrays
template <typename T, typename C = void> struct Traits {}; // C needed for specialisation
// This defines the bottom level - i.e. it's a description of the underlying scalar
template <typename T> struct Traits<T, typename std::enable_if<is_scalar<T>::value, void>::type> {
static constexpr unsigned int depth = 0; // How many levels of Grid Tensor there are (TensorLevel)
static constexpr unsigned int rank = 0; // The rank of the grid tensor (i.e. how many indices used)
static constexpr unsigned int rank_non_trivial = 0; // As per rank, but excludes those of dimension 1
static constexpr unsigned int count = 1; // total number of elements (i.e. product of dimensions)
using scalar_type = T; // Type of the underlying scalar
static constexpr std::size_t scalar_size = sizeof(T); // Size of the underlying scalar in bytes
static constexpr std::size_t size = scalar_size * count; // total size of elements in bytes
static constexpr std::size_t Dimension(unsigned int dim) { return 0; } // Dimension size
static constexpr std::size_t DimensionNT(unsigned int dim) { return 0; } // non-trivial dim size
// e.g. iScalar<iVector<Complex,1>>
// depth = 2
// rank = 1
// rank_non_trivial = 0
// count = 1
// e.g. iVector<iMatrix<Complex,3>,4>
// depth = 2
// rank = 3
// rank_non_trivial = 3
// count = 36
// e.g. iScalar<iVector<iMatrix<Complex,4>,3>>
// depth = 3
// rank = 3
// rank_non_trivial = 3
// count = 48
};
template <typename T> struct Traits<iScalar<T>> {
static constexpr unsigned int depth = 1 + Traits<T>::depth;
static constexpr unsigned int rank = 0 + Traits<T>::rank;
static constexpr unsigned int rank_non_trivial = 0 + Traits<T>::rank_non_trivial;
static constexpr unsigned int count = 1 * Traits<T>::count;
using scalar_type = typename Traits<T>::scalar_type;
static constexpr std::size_t scalar_size = Traits<T>::scalar_size;
static constexpr std::size_t size = scalar_size * count;
static constexpr std::size_t Dimension(unsigned int dim) {
return ( dim == 0 ) ? 1 : Traits<T>::Dimension(dim - 1); }
static constexpr std::size_t DimensionNT(unsigned int dim) {
return Traits<T>::DimensionNT(dim); }
};
template <typename T, int N> struct Traits<iVector<T, N>> {
static constexpr unsigned int depth = 1 + Traits<T>::depth;
static constexpr unsigned int rank = 1 + Traits<T>::rank;
static constexpr unsigned int rank_non_trivial = (N>1 ? 1 : 0) + Traits<T>::rank_non_trivial;
static constexpr unsigned int count = N * Traits<T>::count;
using scalar_type = typename Traits<T>::scalar_type;
static constexpr std::size_t scalar_size = Traits<T>::scalar_size;
static constexpr std::size_t size = scalar_size * count;
static constexpr std::size_t Dimension(unsigned int dim) {
return ( dim == 0 ) ? N : Traits<T>::Dimension(dim - 1); }
static constexpr std::size_t DimensionNT(unsigned int dim) {
if( N != 1 ) {
if( dim == 0 ) return N;
dim--;
}
return Traits<T>::DimensionNT(dim);
}
};
template <typename T, int N> struct Traits<iMatrix<T, N>> {
static constexpr unsigned int depth = 1 + Traits<T>::depth;
static constexpr unsigned int rank = 2 + Traits<T>::rank;
static constexpr unsigned int rank_non_trivial = (N>1 ? 2 : 0) + Traits<T>::rank_non_trivial;
static constexpr unsigned int count = N * N * Traits<T>::count;
using scalar_type = typename Traits<T>::scalar_type;
static constexpr std::size_t scalar_size = Traits<T>::scalar_size;
static constexpr std::size_t size = scalar_size * count;
static constexpr std::size_t Dimension(unsigned int dim) {
return ( dim == 0 || dim == 1 ) ? N : Traits<T>::Dimension(dim - 2); }
static constexpr std::size_t DimensionNT(unsigned int dim) {
if( N != 1 ) {
if( dim == 0 || dim == 1 ) return N;
dim -= 2;
}
return Traits<T>::DimensionNT(dim);
}
};
template <typename T, int N> struct Traits<std::array<T, N>> : Traits<iVector<T, N>> {};
// Tensors have the same traits as their top-level scalar
// Shouldn't be necessary ... but I make the mistake of getting traits of the tensor so often
// that I am tempted to define this.
// HOWEVER, Eigen tensors have a dynamic size flavour, but the scalars are (currently) all fixed size
//template <typename T> struct Traits<T, typename std::enable_if<is_tensor<T>::value, void>::type> : Traits<T> {};
}
// Helper to allow iteration through an Eigen::Tensor (using a lambda) // Helper to allow iteration through an Eigen::Tensor (using a lambda)
template <typename ETensor, typename Lambda> template <typename ETensor, typename Lambda>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && is_grid_tensor<typename ETensor::Scalar>::value, void>::type typename std::enable_if<EigenIO::is_tensor<ETensor>::value, void>::type
for_all( ETensor &ET, Lambda lambda ) for_all( ETensor &ET, Lambda lambda )
{ {
using Scalar = typename ETensor::Scalar::scalar_type; using Scalar = typename ETensor::Scalar; // This could be a Container - we'll check later
const std::size_t NumElements = ET.size(); const std::size_t NumScalars = ET.size();
assert( NumElements > 0 ); assert( NumScalars > 0 );
if( NumElements == 1 ) { // Assemble a vector containing all the non-trivial dimensions (i.e. dimensions > 1)
const auto MyRank{grid_tensor_att<typename ETensor::Scalar>::rank_non_trivial}; unsigned int ScalarElementCount{1};
std::vector<std::size_t> SubIndex(MyRank); std::vector<size_t> NonTrivialDims;
for( auto &idx : SubIndex ) idx = 0; NonTrivialDims.reserve(ET.NumDimensions + EigenIO::Traits<Scalar>::rank_non_trivial);
typename ETensor::Index n = 0;
for( Scalar &Source : * ET.data() ) {
lambda(Source, n++, &SubIndex[0] );
// Now increment SubIndex
for( auto i = MyRank - 1; i >= 0 && ++SubIndex[i] == 11/*ReducedDims[i]*/; i-- )
SubIndex[i] = 0;
}
}
else {
// We're only interested in non-trivial dimensions (i.e. dimensions > 1)
unsigned int TrivialDimCount{0};
std::vector<size_t> ReducedDims;
ReducedDims.reserve(ET.NumDimensions + grid_tensor_att<typename ETensor::Scalar>::rank_non_trivial); // Make sure we only do one malloc
for(auto i = 0; i < ET.NumDimensions; i++ ) { for(auto i = 0; i < ET.NumDimensions; i++ ) {
auto dim = ET.dimension(i); auto dim = ET.dimension(i);
if( dim <= 1 ) { if( dim <= 1 ) {
TrivialDimCount++;
assert( dim == 1 ); // Not expecting dimension to be <= 0 assert( dim == 1 ); // Not expecting dimension to be <= 0
} else { } else {
size_t s = static_cast<size_t>(dim); size_t s = static_cast<size_t>(dim);
assert( s == dim ); // check we didn't lose anything in the conversion assert( s == dim ); // check we didn't lose anything in the conversion
ReducedDims.push_back(s); NonTrivialDims.push_back(s);
ScalarElementCount *= s;
} }
} }
// NB: NumElements > 1 implies this is not a scalar, so some dims should be left // Check that the number of containers is correct
assert( ET.NumDimensions > TrivialDimCount ); assert( NumScalars == ScalarElementCount );
// Now add the extra dimensions, based on object zero // If the Scalar is actually a container, add the inner Scalar's non-trivial dimensions
typename TensorToVec<typename ETensor::Scalar>::type ttv = tensorToVec(* ET.data()); size_t InnerScalarCount{1};
Flatten<typename TensorToVec<typename ETensor::Scalar>::type> f(ttv); for(auto i = 0; i < EigenIO::Traits<Scalar>::rank_non_trivial; i++ ) {
const std::vector<size_t> & ExtraDims{f.getDim()}; auto dim = EigenIO::Traits<Scalar>::DimensionNT(i);
assert(ExtraDims.size() == grid_tensor_att<typename ETensor::Scalar>::rank_non_trivial); assert( dim > 1 );
size_t ExtraCount{1}; NonTrivialDims.push_back(dim);
for( auto i : ExtraDims ) { InnerScalarCount *= dim;
assert( i > 0 );
ExtraCount *= i;
ReducedDims.push_back(i);
} }
assert(grid_tensor_att<typename ETensor::Scalar>::count == ExtraCount); assert(EigenIO::Traits<Scalar>::count == InnerScalarCount);
assert(grid_tensor_att<typename ETensor::Scalar>::size == sizeof( typename ETensor::Scalar )); assert(EigenIO::Traits<Scalar>::size == sizeof( Scalar ));
const unsigned int ReducedDimsSize = static_cast<unsigned int>(ReducedDims.size()); const unsigned int NonTrivialDimsSize = static_cast<unsigned int>(NonTrivialDims.size());
assert( ReducedDimsSize == ReducedDims.size() ); assert( NonTrivialDimsSize == NonTrivialDims.size() );
const typename ETensor::Index TotalNumElements = NumElements * ExtraCount; //const typename ETensor::Index TotalNumElements = NumScalars * InnerScalarCount;
std::array<typename ETensor::Index, ETensor::NumIndices> MyIndex; using Index = typename ETensor::Index;
std::array<Index, ETensor::NumIndices> MyIndex;
for( auto &idx : MyIndex ) idx = 0; for( auto &idx : MyIndex ) idx = 0;
std::vector<std::size_t> SubIndex(ReducedDimsSize); std::vector<std::size_t> SubIndex(NonTrivialDimsSize);
for( auto &idx : SubIndex ) idx = 0; for( auto &idx : SubIndex ) idx = 0;
typename ETensor::Index n = 0; Index n = 0;
for( std::size_t j = 0; j < NumElements; j++ ) { for( std::size_t j = 0; j < NumScalars; j++ ) {
for( Scalar &Source : ET( MyIndex ) ) { // if constexpr is C++ 17 ... but otherwise need two specialisations (Container vs Scalar)
lambda(Source, n++, &SubIndex[0] ); if constexpr ( EigenIO::Traits<Scalar>::rank_non_trivial == 0 ) {
lambda(ET( MyIndex ), n++, &SubIndex[0] );
// Now increment SubIndex // Now increment SubIndex
for( long i = ReducedDimsSize - 1; i >= 0 && ++SubIndex[i] == ReducedDims[i]; i-- ) for( auto i = NonTrivialDimsSize - 1; i != -1 && ++SubIndex[i] == NonTrivialDims[i]; i-- )
SubIndex[i] = 0; SubIndex[i] = 0;
} }
else
{
for( typename Scalar::scalar_type &Source : ET( MyIndex ) ) {
lambda(Source, n++, &SubIndex[0] );
// Now increment SubIndex
for( auto i = NonTrivialDimsSize - 1; i != -1 && ++SubIndex[i] == NonTrivialDims[i]; i-- )
SubIndex[i] = 0;
}
}
// Now increment MyIndex // Now increment MyIndex
for( long i = ET.NumDimensions - 1; i >= 0 && ++MyIndex[i] == ET.dimension(i); i-- ) for( auto i = ET.NumDimensions - 1; i != -1 && ++MyIndex[i] == ET.dimension(i); i-- )
MyIndex[i] = 0; MyIndex[i] = 0;
} }
} }
}
// Abstract writer/reader classes ////////////////////////////////////////////
// static polymorphism implemented using CRTP idiom
class Serializable;
// Static abstract writer // Static abstract writer
template <typename T> template <typename T>
@ -140,10 +247,10 @@ namespace Grid {
template <typename U, int N> template <typename U, int N>
void write(const std::string &s, const iMatrix<U, N> &output); void write(const std::string &s, const iMatrix<U, N> &output);
template <typename ETensor> template <typename ETensor>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && is_eigen_tensor_scalar<typename ETensor::Scalar>::value, void>::type typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && EigenIO::is_scalar<typename ETensor::Scalar>::value, void>::type
write(const std::string &s, const ETensor &output); write(const std::string &s, const ETensor &output);
template <typename ETensor> template <typename ETensor>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && is_grid_tensor<typename ETensor::Scalar>::value, void>::type typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && EigenIO::is_container<typename ETensor::Scalar>::value, void>::type
write(const std::string &s, const ETensor &output); write(const std::string &s, const ETensor &output);
void scientificFormat(const bool set); void scientificFormat(const bool set);
@ -253,7 +360,7 @@ namespace Grid {
// Eigen::Tensors of arithmetic/complex base type // Eigen::Tensors of arithmetic/complex base type
template <typename T> template <typename T>
template <typename ETensor> template <typename ETensor>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && is_eigen_tensor_scalar<typename ETensor::Scalar>::value, void>::type typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && EigenIO::is_scalar<typename ETensor::Scalar>::value, void>::type
Writer<T>::write(const std::string &s, const ETensor &output) Writer<T>::write(const std::string &s, const ETensor &output)
{ {
const typename ETensor::Index NumElements{output.size()}; const typename ETensor::Index NumElements{output.size()};
@ -263,8 +370,8 @@ namespace Grid {
else { else {
// We're only interested in non-trivial dimensions (i.e. dimensions > 1) // We're only interested in non-trivial dimensions (i.e. dimensions > 1)
unsigned int TrivialDimCount{0}; unsigned int TrivialDimCount{0};
std::vector<size_t> ReducedDims; std::vector<size_t> NonTrivialDims;
ReducedDims.reserve(output.NumDimensions); // Make sure we only do one malloc NonTrivialDims.reserve(output.NumDimensions); // Make sure we only do one malloc
for(auto i = 0; i < output.NumDimensions; i++ ) { for(auto i = 0; i < output.NumDimensions; i++ ) {
auto dim = output.dimension(i); auto dim = output.dimension(i);
if( dim <= 1 ) { if( dim <= 1 ) {
@ -273,7 +380,7 @@ namespace Grid {
} else { } else {
size_t s = static_cast<size_t>(dim); size_t s = static_cast<size_t>(dim);
assert( s == dim ); // check we didn't lose anything in the conversion assert( s == dim ); // check we didn't lose anything in the conversion
ReducedDims.push_back(s); NonTrivialDims.push_back(s);
} }
} }
// NB: NumElements > 1 implies this is not a scalar, so some dims should be left // NB: NumElements > 1 implies this is not a scalar, so some dims should be left
@ -298,7 +405,7 @@ namespace Grid {
MyIndex[i] = 0; MyIndex[i] = 0;
} }
} }
upcast->template writeMultiDim<typename ETensor::Scalar>(s, ReducedDims, pWriteBuffer, NumElements); upcast->template writeMultiDim<typename ETensor::Scalar>(s, NonTrivialDims, pWriteBuffer, NumElements);
if( pCopyBuffer ) free( pCopyBuffer ); if( pCopyBuffer ) free( pCopyBuffer );
} }
} }
@ -306,7 +413,7 @@ namespace Grid {
// Eigen::Tensors of Grid tensors (iScalar, iVector, iMatrix) // Eigen::Tensors of Grid tensors (iScalar, iVector, iMatrix)
template <typename T> template <typename T>
template <typename ETensor/*, typename U, int N*/> template <typename ETensor/*, typename U, int N*/>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && is_grid_tensor<typename ETensor::Scalar>::value, void>::type typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && EigenIO::is_container<typename ETensor::Scalar>::value, void>::type
Writer<T>::write(const std::string &s, const ETensor &output) Writer<T>::write(const std::string &s, const ETensor &output)
{ {
const typename ETensor::Index NumElements{output.size()}; const typename ETensor::Index NumElements{output.size()};
@ -316,8 +423,8 @@ namespace Grid {
else { else {
// We're only interested in non-trivial dimensions (i.e. dimensions > 1) // We're only interested in non-trivial dimensions (i.e. dimensions > 1)
unsigned int TrivialDimCount{0}; unsigned int TrivialDimCount{0};
std::vector<size_t> ReducedDims; std::vector<size_t> NonTrivialDims;
ReducedDims.reserve(output.NumDimensions + grid_tensor_att<typename ETensor::Scalar>::rank_non_trivial); // Make sure we only do one malloc NonTrivialDims.reserve(output.NumDimensions + EigenIO::Traits<typename ETensor::Scalar>::rank_non_trivial); // Make sure we only do one malloc
for(auto i = 0; i < output.NumDimensions; i++ ) { for(auto i = 0; i < output.NumDimensions; i++ ) {
auto dim = output.dimension(i); auto dim = output.dimension(i);
if( dim <= 1 ) { if( dim <= 1 ) {
@ -326,7 +433,7 @@ namespace Grid {
} else { } else {
size_t s = static_cast<size_t>(dim); size_t s = static_cast<size_t>(dim);
assert( s == dim ); // check we didn't lose anything in the conversion assert( s == dim ); // check we didn't lose anything in the conversion
ReducedDims.push_back(s); NonTrivialDims.push_back(s);
} }
} }
// NB: NumElements > 1 implies this is not a scalar, so some dims should be left // NB: NumElements > 1 implies this is not a scalar, so some dims should be left
@ -335,15 +442,15 @@ namespace Grid {
typename TensorToVec<typename ETensor::Scalar>::type ttv = tensorToVec(* output.data()); typename TensorToVec<typename ETensor::Scalar>::type ttv = tensorToVec(* output.data());
Flatten<typename TensorToVec<typename ETensor::Scalar>::type> f(ttv); Flatten<typename TensorToVec<typename ETensor::Scalar>::type> f(ttv);
const std::vector<size_t> & ExtraDims{f.getDim()}; const std::vector<size_t> & ExtraDims{f.getDim()};
assert(ExtraDims.size() == grid_tensor_att<typename ETensor::Scalar>::rank_non_trivial); assert(ExtraDims.size() == EigenIO::Traits<typename ETensor::Scalar>::rank_non_trivial);
size_t ExtraCount{1}; size_t ExtraCount{1};
for( auto i : ExtraDims ) { for( auto i : ExtraDims ) {
assert( i > 0 ); assert( i > 0 );
ExtraCount *= i; ExtraCount *= i;
ReducedDims.push_back(i); NonTrivialDims.push_back(i);
} }
assert(grid_tensor_att<typename ETensor::Scalar>::count == ExtraCount); assert(EigenIO::Traits<typename ETensor::Scalar>::count == ExtraCount);
assert(grid_tensor_att<typename ETensor::Scalar>::size == sizeof( typename ETensor::Scalar )); assert(EigenIO::Traits<typename ETensor::Scalar>::size == sizeof( typename ETensor::Scalar ));
// If the Tensor isn't in Row-Major order, then we'll need to copy it's data // If the Tensor isn't in Row-Major order, then we'll need to copy it's data
const bool CopyData{ETensor::Layout != Eigen::StorageOptions::RowMajor}; const bool CopyData{ETensor::Layout != Eigen::StorageOptions::RowMajor};
using Scalar = typename ETensor::Scalar::scalar_type; using Scalar = typename ETensor::Scalar::scalar_type;
@ -368,7 +475,7 @@ namespace Grid {
MyIndex[i] = 0; MyIndex[i] = 0;
} }
} }
upcast->template writeMultiDim<Scalar>(s, ReducedDims, pWriteBuffer, TotalNumElements); upcast->template writeMultiDim<Scalar>(s, NonTrivialDims, pWriteBuffer, TotalNumElements);
if( pCopyBuffer ) free( pCopyBuffer ); if( pCopyBuffer ) free( pCopyBuffer );
} }
} }

72
Grid/tensors/Tensor_class.h
View File

@ -174,36 +174,36 @@ class iScalar {
}; };
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return &_internal; } strong_inline begin() const { return &_internal; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return _internal.begin(); } strong_inline begin() const { return _internal.begin(); }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return (&_internal) + grid_tensor_att<iScalar<T>>::count; } strong_inline end() const { return (&_internal) + 1; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return _internal.begin() + grid_tensor_att<iScalar<T>>::count; } strong_inline end() const { return _internal.begin() + sizeof(_internal)/sizeof(scalar_type); }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return &_internal; } strong_inline begin() { return &_internal; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return _internal.begin(); } strong_inline begin() { return _internal.begin(); }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return (&_internal) + grid_tensor_att<iScalar<T>>::count; } strong_inline end() { return (&_internal) + 1; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return _internal.begin() + grid_tensor_att<iScalar<T>>::count; } strong_inline end() { return _internal.begin() + sizeof(_internal)/sizeof(scalar_type); }
}; };
/////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////
// Allows to turn scalar<scalar<scalar<double>>>> back to double. // Allows to turn scalar<scalar<scalar<double>>>> back to double.
@ -335,36 +335,36 @@ class iVector {
// } // }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return _internal; } strong_inline begin() const { return _internal; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return _internal[0].begin(); } strong_inline begin() const { return _internal[0].begin(); }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return _internal + grid_tensor_att<iVector<T,N>>::count; } strong_inline end() const { return _internal + N; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return _internal[0].begin() + grid_tensor_att<iVector<T,N>>::count; } strong_inline end() const { return _internal[0].begin() + sizeof(_internal)/sizeof(scalar_type); }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return _internal; } strong_inline begin() { return _internal; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return _internal[0].begin(); } strong_inline begin() { return _internal[0].begin(); }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return _internal + grid_tensor_att<iVector<T,N>>::count; } strong_inline end() { return _internal + N; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return _internal[0].begin() + grid_tensor_att<iVector<T,N>>::count; } strong_inline end() { return _internal[0].begin() + sizeof(_internal)/sizeof(scalar_type); }
}; };
template <class vtype, int N> template <class vtype, int N>
@ -522,36 +522,36 @@ class iMatrix {
// } // }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return _internal[0]; } strong_inline begin() const { return _internal[0]; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return _internal[0][0].begin(); } strong_inline begin() const { return _internal[0][0].begin(); }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return _internal[0] + grid_tensor_att<iMatrix<T,N>>::count; } strong_inline end() const { return _internal[0] + N * N; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, const scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return _internal[0][0].begin() + grid_tensor_att<iMatrix<T,N>>::count; } strong_inline end() const { return _internal[0][0].begin() + sizeof(_internal)/sizeof(scalar_type); }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return _internal[0]; } strong_inline begin() { return _internal[0]; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return _internal[0][0].begin(); } strong_inline begin() { return _internal[0][0].begin(); }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<!is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return _internal[0] + grid_tensor_att<iMatrix<T,N>>::count; } strong_inline end() { return _internal[0] + N * N; }
template <typename T = vtype> template <typename T = vtype>
typename std::enable_if<is_grid_tensor<T>::value, scalar_type *>::type typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return _internal[0][0].begin() + grid_tensor_att<iMatrix<T,N>>::count; } strong_inline end() { return _internal[0][0].begin() + sizeof(_internal)/sizeof(scalar_type); }
}; };
template <class v> template <class v>

73
Grid/tensors/Tensor_traits.h
View File

@ -288,79 +288,6 @@ namespace Grid {
enum { value = sizeof(real_scalar_type)/sizeof(float) }; enum { value = sizeof(real_scalar_type)/sizeof(float) };
}; };
////////////////////////////////////////////////////////////////////////////
// Review with Peter - obvious (but partial) overlap with the above
// What of this is good and should be kept ... vs what functions should I really be using?
////////////////////////////////////////////////////////////////////////////
// Need some forward references or the below won't work
template <class vtype>
class iScalar;
template <class vtype, int N>
class iVector;
template <class vtype, int N>
class iMatrix;
// which types are grid tensors
template <typename T> struct is_grid_tensor : public std::false_type {};
template <typename T> struct is_grid_tensor<iScalar<T>> : public std::true_type {};
template <typename T, int N> struct is_grid_tensor<iVector<T, N>> : public std::true_type {};
template <typename T, int N> struct is_grid_tensor<iMatrix<T, N>> : public std::true_type {};
// Rank and dimension of grid tensors, i.e. compositions of iScalar, iVector and iMatrix
// This defines the bottom level - i.e. it's a description of the underlying scalar
template <typename T> struct grid_tensor_att {
static constexpr unsigned int depth = 0; // How many levels of Grid Tensor there are (TensorLevel)
static constexpr unsigned int rank = 0; // The rank of the grid tensor (i.e. how many indices used)
static constexpr unsigned int rank_non_trivial = 0; // As per rank, but excludes those of dimension 1
static constexpr unsigned int count = 1; // total number of elements (i.e. product of dimensions)
using scalar_type = T; // Type of the underlying scalar
static constexpr std::size_t scalar_size = sizeof(T); // Size of the underlying scalar in bytes
static constexpr std::size_t size = scalar_size * count; // total size of elements in bytes
// e.g. iScalar<iVector<Complex,1>>
// depth = 2
// rank = 1
// rank_non_trivial = 0
// count = 1
// e.g. iVector<iMatrix<Complex,3>,4>
// depth = 2
// rank = 3
// rank_non_trivial = 3
// count = 36
// e.g. iScalar<iVector<iMatrix<Complex,4>,3>>
// depth = 3
// rank = 3
// rank_non_trivial = 3
// count = 48
};
template <typename T> struct grid_tensor_att<iScalar<T>> {
static constexpr unsigned int depth = 1 + grid_tensor_att<T>::depth;
static constexpr unsigned int rank = 0 + grid_tensor_att<T>::rank;
static constexpr unsigned int rank_non_trivial = 0 + grid_tensor_att<T>::rank_non_trivial;
static constexpr unsigned int count = 1 * grid_tensor_att<T>::count;
using scalar_type = typename grid_tensor_att<T>::scalar_type;
static constexpr std::size_t scalar_size = grid_tensor_att<T>::scalar_size;
static constexpr std::size_t size = scalar_size * count;
};
template <typename T, int N> struct grid_tensor_att<iVector<T, N>> {
static constexpr unsigned int depth = 1 + grid_tensor_att<T>::depth;
static constexpr unsigned int rank = 1 + grid_tensor_att<T>::rank;
static constexpr unsigned int rank_non_trivial = (N>1 ? 1 : 0) + grid_tensor_att<T>::rank_non_trivial;
static constexpr unsigned int count = N * grid_tensor_att<T>::count;
using scalar_type = typename grid_tensor_att<T>::scalar_type;
static constexpr std::size_t scalar_size = grid_tensor_att<T>::scalar_size;
static constexpr std::size_t size = scalar_size * count;
};
template <typename T, int N> struct grid_tensor_att<iMatrix<T, N>> {
static constexpr unsigned int depth = 1 + grid_tensor_att<T>::depth;
static constexpr unsigned int rank = 2 + grid_tensor_att<T>::rank;
static constexpr unsigned int rank_non_trivial = (N>1 ? 2 : 0) + grid_tensor_att<T>::rank_non_trivial;
static constexpr unsigned int count = N * N * grid_tensor_att<T>::count;
using scalar_type = typename grid_tensor_att<T>::scalar_type;
static constexpr std::size_t scalar_size = grid_tensor_att<T>::scalar_size;
static constexpr std::size_t size = scalar_size * count;
};
} }
#endif #endif

28
tests/IO/Test_serialisation.cc
View File

@ -101,6 +101,20 @@ void ioTest(const std::string &filename, const O &object, const std::string &nam
} }
#ifdef DEBUG #ifdef DEBUG
template <typename T>
//typename std::enable_if<EigenIO::is_tensor<T>::value, void>::type
void dump_tensor(T & t)
{
using Traits = typename EigenIO::Traits<typename T::Scalar>;
for_all( t, [&](typename Traits::scalar_type &c, typename T::Index index, const std::size_t * pDims ){
std::cout << " ";
for( int i = 0 ; i < t.NumDimensions + Traits::rank_non_trivial; i++ )
std::cout << "[" << pDims[i] << "]";
std::cout << " = " << c << std::endl;
} );
std::cout << "========================================" << std::endl;
}
//typedef int TestScalar; //typedef int TestScalar;
typedef std::complex<double> TestScalar; typedef std::complex<double> TestScalar;
typedef Eigen::Tensor<iMatrix<TestScalar,1>, 6> TestTensorSingle; typedef Eigen::Tensor<iMatrix<TestScalar,1>, 6> TestTensorSingle;
@ -145,6 +159,8 @@ bool EigenIOTest(void) {
Val += Inc; Val += Inc;
} }
ioTest<Hdf5Writer, Hdf5Reader, TestTensor>("iotest_tensor.h5", t, "eigen_tensor_instance_name"); ioTest<Hdf5Writer, Hdf5Reader, TestTensor>("iotest_tensor.h5", t, "eigen_tensor_instance_name");
std::cout << "t:";
dump_tensor(t);
// Now serialise a fixed size tensor // Now serialise a fixed size tensor
using FixedTensor = Eigen::TensorFixedSize<TestScalar, Eigen::Sizes<8,4,3>>; using FixedTensor = Eigen::TensorFixedSize<TestScalar, Eigen::Sizes<8,4,3>>;
@ -157,6 +173,8 @@ bool EigenIOTest(void) {
Val += Inc; Val += Inc;
} }
ioTest<Hdf5Writer, Hdf5Reader, FixedTensor>("iotest_tensor_fixed.h5", tf, "eigen_tensor_fixed_name"); ioTest<Hdf5Writer, Hdf5Reader, FixedTensor>("iotest_tensor_fixed.h5", tf, "eigen_tensor_fixed_name");
std::cout << "tf:";
dump_tensor(tf);
ETSerClass o; ETSerClass o;
ioTest<Hdf5Writer, Hdf5Reader, ETSerClass>("iotest_object.h5", o, "ETSerClass_object_instance_name"); ioTest<Hdf5Writer, Hdf5Reader, ETSerClass>("iotest_object.h5", o, "ETSerClass_object_instance_name");
@ -174,14 +192,8 @@ bool EigenIOTest(void) {
Val += Inc; Val += Inc;
} }
ioTest<Hdf5Writer, Hdf5Reader, LSCTensor>("iotest_LSCTensor.h5", l, "LSCTensor_object_instance_name"); ioTest<Hdf5Writer, Hdf5Reader, LSCTensor>("iotest_LSCTensor.h5", l, "LSCTensor_object_instance_name");
std::cout << "t:"; std::cout << "l:";
for_all( l, [](std::complex<double> &c, LSCTensor::Index index, const std::size_t * pDims ){ dump_tensor(l);
std::cout << " ";
for( int i = 0 ; i < 5; i++ )
std::cout << "[" << pDims[i] << "]";
std::cout << " = " << c << std::endl;
} );
std::cout << std::endl;
// Tensor of spin colour // Tensor of spin colour
LCMTensor l2; LCMTensor l2;

30
tests/hadrons/Test_hadrons_distil.cc
View File

@ -467,13 +467,13 @@ bool DebugEigenTest()
template <typename T> template <typename T>
void DebugGridTensorTest_print( int i ) void DebugGridTensorTest_print( int i )
{ {
std::cout << i << " : " << is_grid_tensor<T>::value std::cout << i << " : " << EigenIO::is_tensor<T>::value
<< ", depth " << grid_tensor_att<T>::depth << ", depth " << EigenIO::Traits<T>::depth
<< ", rank " << grid_tensor_att<T>::rank << ", rank " << EigenIO::Traits<T>::rank
<< ", rank_non_trivial " << grid_tensor_att<T>::rank_non_trivial << ", rank_non_trivial " << EigenIO::Traits<T>::rank_non_trivial
<< ", count " << grid_tensor_att<T>::count << ", count " << EigenIO::Traits<T>::count
<< ", scalar_size " << grid_tensor_att<T>::scalar_size << ", scalar_size " << EigenIO::Traits<T>::scalar_size
<< ", size " << grid_tensor_att<T>::size << ", size " << EigenIO::Traits<T>::size
<< std::endl; << std::endl;
} }
@ -494,8 +494,24 @@ public:
inline value_type * end(void) { return m_p + N; } inline value_type * end(void) { return m_p + N; }
}; };
bool DebugFelixTensorTest( void )
{
unsigned int Nmom = 2;
unsigned int Nt = 2;
unsigned int N_1 = 2;
unsigned int N_2 = 2;
unsigned int N_3 = 2;
using BaryonTensorSet = Eigen::Tensor<Complex, 6, Eigen::RowMajor>;
BaryonTensorSet BField3(Nmom,4,Nt,N_1,N_2,N_3);
std::vector<Complex> Memory(Nmom * Nt * N_1 * N_2 * N_3 * 2);
using BaryonTensorMap = Eigen::TensorMap<BaryonTensorSet>;
BaryonTensorMap BField4 (&Memory[0], Nmom,4,Nt,N_1,N_2,N_3);
return true;
}
bool DebugGridTensorTest( void ) bool DebugGridTensorTest( void )
{ {
DebugFelixTensorTest();
typedef Complex t1; typedef Complex t1;
typedef iScalar<t1> t2; typedef iScalar<t1> t2;
typedef iVector<t1, Ns> t3; typedef iVector<t1, Ns> t3;