portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-10-31 03:54:33 +00:00
Code Releases Activity

Merge branch 'feature/distil' of github.com:mmphys/Grid into feature/distil

Browse Source
This commit is contained in:
ferben
2019-02-15 10:47:41 +00:00
parent 9f2ca98dfc e8bd8767c0
commit a111d814db
6 changed files with 584 additions and 367 deletions
Download Patch File Download Diff File Expand all files Collapse all files

397
Grid/serialisation/BaseIO.h
View File

@@ -36,75 +36,186 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
#include <Grid/Eigen/unsupported/CXX11/Tensor>
namespace Grid {
namespace EigenIO {
template<typename T> struct is_complex : public std::false_type {};
template<typename T> struct is_complex<std::complex<T>>
: std::integral_constant<bool, std::is_arithmetic<T>::value> {};
// Eigen tensors can be composed of arithmetic scalar and complex types
template<typename T> struct is_scalar : std::integral_constant<bool,
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) {
return ( N == 1 ) ? Traits<T>::DimensionNT(dim) : ( dim == 0 ) ? N : Traits<T>::DimensionNT(dim - 1);
}
};
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) {
return ( N == 1 ) ? Traits<T>::DimensionNT(dim) : ( dim == 0 || dim == 1 ) ? N : Traits<T>::DimensionNT(dim - 2);
}
};
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)
template <typename ETensor, typename Lambda>
typename std::enable_if<EigenIO::is_tensor<ETensor>::value, void>::type
for_all( ETensor &ET, Lambda lambda )
{
using Scalar = typename ETensor::Scalar; // This could be a Container - we'll check later
const std::size_t NumScalars = ET.size();
assert( NumScalars > 0 );
// Assemble a vector containing all the non-trivial dimensions (i.e. dimensions > 1)
unsigned int ScalarElementCount{1};
std::vector<size_t> NonTrivialDims;
NonTrivialDims.reserve(ET.NumDimensions + EigenIO::Traits<Scalar>::rank_non_trivial);
for(auto i = 0; i < ET.NumDimensions; i++ ) {
auto dim = ET.dimension(i);
if( dim <= 1 ) {
assert( dim == 1 ); // Not expecting dimension to be <= 0
} else {
size_t s = static_cast<size_t>(dim);
assert( s == dim ); // check we didn't lose anything in the conversion
NonTrivialDims.push_back(s);
ScalarElementCount *= s;
}
}
// Check that the number of containers is correct
assert( NumScalars == ScalarElementCount );
// If the Scalar is actually a container, add the inner Scalar's non-trivial dimensions
size_t InnerScalarCount{1};
for(auto i = 0; i < EigenIO::Traits<Scalar>::rank_non_trivial; i++ ) {
auto dim = EigenIO::Traits<Scalar>::DimensionNT(i);
assert( dim > 1 );
NonTrivialDims.push_back(dim);
InnerScalarCount *= dim;
}
assert(EigenIO::Traits<Scalar>::count == InnerScalarCount);
assert(EigenIO::Traits<Scalar>::size == sizeof( Scalar ));
const unsigned int NonTrivialDimsSize = static_cast<unsigned int>(NonTrivialDims.size());
assert( NonTrivialDimsSize == NonTrivialDims.size() );
//const typename ETensor::Index TotalNumElements = NumScalars * InnerScalarCount;
using Index = typename ETensor::Index;
std::array<Index, ETensor::NumIndices> MyIndex;
for( auto &idx : MyIndex ) idx = 0;
std::vector<std::size_t> SubIndex(NonTrivialDimsSize);
for( auto &idx : SubIndex ) idx = 0;
Index n = 0;
for( std::size_t j = 0; j < NumScalars; j++ ) {
// if constexpr is C++ 17 ... but otherwise need two specialisations (Container vs Scalar)
if constexpr ( EigenIO::Traits<Scalar>::rank_non_trivial == 0 ) {
lambda(ET( MyIndex ), n++, &SubIndex[0] );
// Now increment SubIndex
for( auto i = NonTrivialDimsSize - 1; i != -1 && ++SubIndex[i] == NonTrivialDims[i]; i-- )
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
for( auto i = ET.NumDimensions - 1; i != -1 && ++MyIndex[i] == ET.dimension(i); i-- )
MyIndex[i] = 0;
}
}
// Abstract writer/reader classes ////////////////////////////////////////////
// static polymorphism implemented using CRTP idiom
class Serializable;
// which types are supported scalar types for Eigen::Tensor
template<typename T> struct is_eigen_tensor_scalar : std::integral_constant<bool,
std::is_arithmetic<T>::value || is_complex<T>::value> {};
// which types are grid tensors
template <typename T> struct is_grid_tensor : public std::false_type {
static constexpr unsigned int rank = 0;
static constexpr unsigned int dim = 1;
};
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
template <typename T> struct grid_tensor_att {
static constexpr unsigned int depth = 0; // How many levels of Grid Tensor there are
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)
typedef T scalar_type; // Type of the underlying scalar
static constexpr size_t scalar_size = sizeof(T); // Size of the underlying scalar in bytes
static constexpr 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;
typedef typename grid_tensor_att<T>::scalar_type scalar_type;
static constexpr size_t scalar_size = grid_tensor_att<T>::scalar_size;
static constexpr 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;
typedef typename grid_tensor_att<T>::scalar_type scalar_type;
static constexpr size_t scalar_size = grid_tensor_att<T>::scalar_size;
static constexpr 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;
typedef typename grid_tensor_att<T>::scalar_type scalar_type;
static constexpr size_t scalar_size = grid_tensor_att<T>::scalar_size;
static constexpr size_t size = scalar_size * count;
};
// Static abstract writer
template <typename T>
class Writer
@@ -128,24 +239,11 @@ namespace Grid {
template <typename U, int N>
void write(const std::string &s, const iMatrix<U, N> &output);
template <typename ETensor>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && (std::is_arithmetic<typename ETensor::Scalar>::value || Grid::is_complex<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);
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
//&& !(std::is_arithmetic<typename ETensor::Scalar>::value || Grid::is_complex<typename ETensor::Scalar>::value)
/*&& ( std::is_base_of<typename ETensor::Scalar, iScalar<U> >::value
|| std::is_base_of<typename ETensor::Scalar, iVector<U, N>>::value
|| std::is_base_of<typename ETensor::Scalar, iMatrix<U, N>>::value )*/, void>::type
template <typename ETensor>
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);
/*template <typename ETensor, typename U, int N>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value
&& std::is_base_of<typename ETensor::Scalar, iVector<U, N>>::value, void>::type
write(const std::string &s, const ETensor &output);
template <typename ETensor, typename U, int N>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value
&& std::is_base_of<typename ETensor::Scalar, iMatrix<U, N>>::value, void>::type
write(const std::string &s, const ETensor &output);*/
void scientificFormat(const bool set);
bool isScientific(void);
@@ -254,18 +352,18 @@ namespace Grid {
// Eigen::Tensors of arithmetic/complex base type
template <typename T>
template <typename ETensor>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && (std::is_arithmetic<typename ETensor::Scalar>::value || Grid::is_complex<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)
{
std::cout << "Eigen::Tensors of arithmetic/complex base type" << std::endl;
const typename ETensor::Index NumElements{output.size()};
assert( NumElements > 0 );
if( NumElements == 1 )
upcast->writeDefault(s, * output.data());
else {
// We're not interested in trivial dimensions, i.e. dimensions = 1
// We're only interested in non-trivial dimensions (i.e. dimensions > 1)
unsigned int TrivialDimCount{0};
std::vector<size_t> ReducedDims;
std::vector<size_t> NonTrivialDims;
NonTrivialDims.reserve(output.NumDimensions); // Make sure we only do one malloc
for(auto i = 0; i < output.NumDimensions; i++ ) {
auto dim = output.dimension(i);
if( dim <= 1 ) {
@@ -274,46 +372,51 @@ namespace Grid {
} else {
size_t s = static_cast<size_t>(dim);
assert( s == dim ); // check we didn't lose anything in the conversion
ReducedDims.push_back(s);
NonTrivialDims.push_back(s);
}
}
assert( output.NumDimensions > TrivialDimCount > 0 ); // NB: NumElements > 1 implies this is not a scalar, so some dims should be left
// Create a single, flat vector to hold all the data
std::vector<typename ETensor::Scalar> flat(NumElements);
// Now copy all the data to my flat vector
// Regardless of the Eigen::Tensor storage order, the copy will be Row Major
std::array<typename ETensor::Index, ETensor::NumIndices> MyIndex;
for( int i = 0 ; i < output.NumDimensions ; i++ ) MyIndex[i] = 0;
for( typename ETensor::Index n = 0; n < NumElements; n++ ) {
flat[n] = output( MyIndex );
// Now increment the index
for( int i = output.NumDimensions - 1; i >= 0 && ++MyIndex[i] == output.dimension(i); i-- )
MyIndex[i] = 0;
// NB: NumElements > 1 implies this is not a scalar, so some dims should be left
assert( output.NumDimensions > TrivialDimCount );
// 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};
using Scalar = typename ETensor::Scalar;
const Scalar * pWriteBuffer;
Scalar * pCopyBuffer = nullptr;
if( !CopyData )
pWriteBuffer = output.data();
else {
// Regardless of the Eigen::Tensor storage order, the copy will be Row Major
pCopyBuffer = static_cast<Scalar *>(malloc(sizeof(Scalar) * NumElements));
pWriteBuffer = pCopyBuffer;
std::array<typename ETensor::Index, ETensor::NumIndices> MyIndex;
for( auto &idx : MyIndex ) idx = 0;
for( typename ETensor::Index n = 0; n < NumElements; n++ ) {
pCopyBuffer[n] = output( MyIndex );
// Now increment the index
for( int i = output.NumDimensions - 1; i >= 0 && ++MyIndex[i] == output.dimension(i); i-- )
MyIndex[i] = 0;
}
}
upcast->template writeMultiDim<typename ETensor::Scalar>(s, ReducedDims, flat);
upcast->template writeMultiDim<typename ETensor::Scalar>(s, NonTrivialDims, pWriteBuffer, NumElements);
if( pCopyBuffer ) free( pCopyBuffer );
}
}
// Eigen::Tensors of iScalar<U>
// Eigen::Tensors of Grid tensors (iScalar, iVector, iMatrix)
template <typename T>
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
//&& !(std::is_arithmetic<typename ETensor::Scalar>::value || Grid::is_complex<typename ETensor::Scalar>::value)
/*&& ( std::is_base_of<typename ETensor::Scalar, iScalar<U> >::value
|| std::is_base_of<typename ETensor::Scalar, iVector<U, N>>::value
|| std::is_base_of<typename ETensor::Scalar, iMatrix<U, N>>::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)
{
std::cout << "Eigen::Tensors of iScalar<U>" << std::endl;
const typename ETensor::Index NumElements{output.size()};
assert( NumElements > 0 );
if( NumElements == 1 )
upcast->writeDefault(s, tensorToVec(* output.data()));
else {
// We're not interested in trivial dimensions, i.e. dimensions = 1
// We're only interested in non-trivial dimensions (i.e. dimensions > 1)
unsigned int TrivialDimCount{0};
std::vector<size_t> ReducedDims;
std::vector<size_t> NonTrivialDims;
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++ ) {
auto dim = output.dimension(i);
if( dim <= 1 ) {
@@ -322,63 +425,53 @@ namespace Grid {
} else {
size_t s = static_cast<size_t>(dim);
assert( s == dim ); // check we didn't lose anything in the conversion
ReducedDims.push_back(s);
NonTrivialDims.push_back(s);
}
}
assert( output.NumDimensions > TrivialDimCount > 0 ); // 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
assert( output.NumDimensions > TrivialDimCount );
// Now add the extra dimensions, based on object zero
typename TensorToVec<typename ETensor::Scalar>::type ttv = tensorToVec(* output.data());
Flatten<typename TensorToVec<typename ETensor::Scalar>::type> f(ttv);
const std::vector<size_t> & ExtraDims{f.getDim()};
assert(ExtraDims.size() == EigenIO::Traits<typename ETensor::Scalar>::rank_non_trivial);
size_t ExtraCount{1};
for( auto i : ExtraDims ) {
assert( i > 0 );
ExtraCount *= i;
ReducedDims.push_back(i);
NonTrivialDims.push_back(i);
}
typedef typename ETensor::Scalar::scalar_type Scalar;
assert( sizeof( typename ETensor::Scalar ) == ExtraCount * sizeof( Scalar ) );
// Create a single, flat vector to hold all the data
assert(EigenIO::Traits<typename ETensor::Scalar>::count == ExtraCount);
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
const bool CopyData{ETensor::Layout != Eigen::StorageOptions::RowMajor};
using Scalar = typename ETensor::Scalar::scalar_type;
const Scalar * pWriteBuffer;
Scalar * pCopyBuffer = nullptr;
const typename ETensor::Index TotalNumElements = NumElements * ExtraCount;
std::vector<Scalar> flat(TotalNumElements);
// Now copy all the data to my flat vector
// Regardless of the Eigen::Tensor storage order, the copy will be Row Major
std::array<typename ETensor::Index, ETensor::NumIndices> MyIndex;
for( int i = 0 ; i < output.NumDimensions ; i++ ) MyIndex[i] = 0;
for( typename ETensor::Index n = 0; n < TotalNumElements; ) {
const Scalar * p = reinterpret_cast<const Scalar *>( &output( MyIndex ));
for( auto j = 0; j < ExtraCount ; j++ )
flat[n++] = * p++;
// Now increment the index
for( int i = output.NumDimensions - 1; i >= 0 && ++MyIndex[i] == output.dimension(i); i-- )
MyIndex[i] = 0;
if( !CopyData )
pWriteBuffer = output.data()->begin();
else {
// Regardless of the Eigen::Tensor storage order, the copy will be Row Major
pCopyBuffer = static_cast<Scalar *>(malloc(TotalNumElements * sizeof(Scalar)));
pWriteBuffer = pCopyBuffer;
Scalar * pCopy = pCopyBuffer;
std::array<typename ETensor::Index, ETensor::NumIndices> MyIndex;
for( auto &idx : MyIndex ) idx = 0;
for( typename ETensor::Index n = 0; n < NumElements; n++ ) {
// Copy the grid tensor
for( const Scalar &Source : output( MyIndex ) )
* pCopy ++ = Source;
// Now increment the index
for( int i = output.NumDimensions - 1; i >= 0 && ++MyIndex[i] == output.dimension(i); i-- )
MyIndex[i] = 0;
}
}
upcast->template writeMultiDim<Scalar>(s, ReducedDims, flat);
upcast->template writeMultiDim<Scalar>(s, NonTrivialDims, pWriteBuffer, TotalNumElements);
if( pCopyBuffer ) free( pCopyBuffer );
}
}
// Eigen::Tensors of iVector<U, N>
/*template <typename T>
template <typename ETensor, typename U, int N>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value
&& std::is_base_of<typename ETensor::Scalar, iVector<U, N>>::value, void>::type
Writer<T>::write(const std::string &s, const ETensor &output)
{
//upcast->writeDefault(s, tensorToVec(output));
std::cout << "I really should add code to write Eigen::Tensor (iVector) ..." << std::endl;
}
// Eigen::Tensors of iMatrix<U, N>
template <typename T>
template <typename ETensor, typename U, int N>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value
&& std::is_base_of<typename ETensor::Scalar, iMatrix<U, N>>::value, void>::type
Writer<T>::write(const std::string &s, const ETensor &output)
{
//upcast->writeDefault(s, tensorToVec(output));
std::cout << "I really should add code to write Eigen::Tensor (iMatrix) ..." << std::endl;
}*/
template <typename T>
void Writer<T>::scientificFormat(const bool set)
{
@@ -516,7 +609,7 @@ namespace Grid {
template <typename T>
static inline typename std::enable_if<std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value, bool>::type
CompareMember(const T &lhs, const T &rhs) {
Eigen::Tensor<bool, 0> bResult = (lhs == rhs).all();
Eigen::Tensor<bool, 0, T::Options> bResult = (lhs == rhs).all();
return bResult(0);
}
@@ -526,7 +619,7 @@ namespace Grid {
const auto NumElements{lhs.size()};
bool bResult = ( NumElements == rhs.size() );
for( auto i = 0 ; i < NumElements && bResult ; i++ ) {
Eigen::Tensor<bool, 0> b = (lhs[i] == rhs[i]).all();
Eigen::Tensor<bool, 0, T::Options> b = (lhs[i] == rhs[i]).all();
bResult = b(0);
}
return bResult;

27
Grid/serialisation/Hdf5IO.h
View File

@@ -39,7 +39,7 @@ namespace Grid
typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
writeDefault(const std::string &s, const std::vector<U> &x);
template <typename U>
void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const std::vector<U> & DataRowMajor);
void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
H5NS::Group & getGroup(void);
private:
template <typename U>
@@ -104,31 +104,32 @@ namespace Grid
void Hdf5Writer::writeDefault(const std::string &s, const std::string &x);
template <typename U>
void Hdf5Writer::writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const std::vector<U> & DataRowMajor)
void Hdf5Writer::writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements)
{
// Hdf5 needs the dimensions as hsize_t
std::vector<hsize_t> dim;
for (auto &d: Dimensions)
dim.push_back(d);
int rank = static_cast<int>(Dimensions.size());
std::vector<hsize_t> dim(rank);
for(int i = 0; i < rank; i++)
dim[i] = Dimensions[i];
// write to file
H5NS::DataSpace dataSpace(dim.size(), dim.data());
if (DataRowMajor.size() > dataSetThres_)
H5NS::DataSpace dataSpace(rank, dim.data());
size_t DataSize = NumElements * sizeof(U);
if (DataSize > dataSetThres_)
{
H5NS::DataSet dataSet;
H5NS::DSetCreatPropList plist;
plist.setChunk(dim.size(), dim.data());
plist.setChunk(rank, dim.data());
plist.setFletcher32();
dataSet = group_.createDataSet(s, Hdf5Type<U>::type(), dataSpace, plist);
dataSet.write(&DataRowMajor[0], Hdf5Type<U>::type());
dataSet.write(pDataRowMajor, Hdf5Type<U>::type());
}
else
{
H5NS::Attribute attribute;
attribute = group_.createAttribute(s, Hdf5Type<U>::type(), dataSpace);
attribute.write(Hdf5Type<U>::type(), &DataRowMajor[0]);
attribute.write(Hdf5Type<U>::type(), pDataRowMajor);
}
}
@@ -145,7 +146,7 @@ namespace Grid
const auto &flatx = flat.getFlatVector();
for (auto &d: flat.getDim())
dim.push_back(d);
writeMultiDim<Element>(s, dim, flatx);
writeMultiDim<Element>(s, dim, &flatx[0], flatx.size());
}
template <typename U>

94
Grid/tensors/Tensor_class.h
View File

@@ -173,7 +173,37 @@ class iScalar {
return stream;
};
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return &_internal; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return _internal.begin(); }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return (&_internal) + 1; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return _internal.begin() + sizeof(_internal)/sizeof(scalar_type); }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return &_internal; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return _internal.begin(); }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return (&_internal) + 1; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return _internal.begin() + sizeof(_internal)/sizeof(scalar_type); }
};
///////////////////////////////////////////////////////////
// Allows to turn scalar<scalar<scalar<double>>>> back to double.
@@ -303,6 +333,38 @@ class iVector {
// strong_inline vtype && operator ()(int i) {
// return _internal[i];
// }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return _internal; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return _internal[0].begin(); }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return _internal + N; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return _internal[0].begin() + sizeof(_internal)/sizeof(scalar_type); }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return _internal; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return _internal[0].begin(); }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return _internal + N; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return _internal[0].begin() + sizeof(_internal)/sizeof(scalar_type); }
};
template <class vtype, int N>
@@ -458,6 +520,38 @@ class iMatrix {
// strong_inline vtype && operator ()(int i,int j) {
// return _internal[i][j];
// }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return _internal[0]; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline begin() const { return _internal[0][0].begin(); }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return _internal[0] + N * N; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
strong_inline end() const { return _internal[0][0].begin() + sizeof(_internal)/sizeof(scalar_type); }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return _internal[0]; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline begin() { return _internal[0][0].begin(); }
template <typename T = vtype>
typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return _internal[0] + N * N; }
template <typename T = vtype>
typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
strong_inline end() { return _internal[0][0].begin() + sizeof(_internal)/sizeof(scalar_type); }
};
template <class v>

89
Hadrons/Modules/MDistil/BContraction.hpp
View File

@@ -18,6 +18,20 @@ BEGIN_HADRONS_NAMESPACE
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MDistil)
// general baryon tensor set based on Eigen tensors and Grid-allocated memory
// Dimensions:
// 0 - ext - external field (momentum, EM field, ...)
// 1 - str - spin-color structure
// 2 - t - timeslice
// 3 - i - left distillation mode index
// 4 - j - middle distillation mode index
// 5 - k - left distillation mode index
// template <typename T>
// using BaryonTensorSet = Eigen::TensorMap<Eigen::Tensor<T, 6, Eigen::RowMajor>>;
class BContractionPar: Serializable
{
public:
@@ -141,6 +155,7 @@ void TBContraction<FImpl>::execute(void)
SpinVector tmp222;
SpinVector tmp111;
assert(parity == 1 || parity == -1);
std::vector<std::vector<int>> epsilon = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
@@ -159,61 +174,11 @@ void TBContraction<FImpl>::execute(void)
Gamma::Algebra::GammaYGamma5, // i gamma_4 C gamma_5 = i gamma_2 gamma_5
};
std::vector<Complex> factor23{{0.,-1.},{0.,1.},{0.,1.}};
if((0)){
for (int i1=0 ; i1 < N_1 ; i1++){
for (int i2=0 ; i2 < N_2 ; i2++){
for (int i3=0 ; i3 < N_3 ; i3++){
for (int imom=0 ; imom < Nmom ; imom++){
for (int t=0 ; t < Nt ; t++){
Bindex = i1 + N_1*(i2 + N_2*(i3 + N_3*(imom+Nmom*t)));
ExtractSliceLocal(tmp1,one[i1],0,t,3);
ExtractSliceLocal(tmp2,two[i2],0,t,3);
ExtractSliceLocal(tmp3,three[i3],0,t,3);
parallel_for (unsigned int sU = 0; sU < grid3d->oSites(); ++sU)
{
for (int ie=0 ; ie < 6 ; ie++){
// Why does peekColour not work????
for (int is=0 ; is < 4 ; is++){
tmp11s()(is)() = tmp11[sU]()(is)(epsilon[ie][0]);
tmp22s()(is)() = tmp22[sU]()(is)(epsilon[ie][1]);
tmp33s()(is)() = tmp33[sU]()(is)(epsilon[ie][2]);
}
tmp333 = Gamma(gamma23_[0])*tmp33s;
// this should be outerProduct??? Does not work.
for (int isl=0 ; isl < 4 ; isl++){
for (int isr=0 ; isr < 4 ; isr++){
diquark()(isl,isr)() = factor23[0]*tmp22s()(isl)()*tmp333()(isr)();
}
}
// Is there a way to compute gamma * SpinMatrix (left component)???
for (int isr=0 ; isr < 4 ; isr++){
for (int isl=0 ; isl < 4 ; isl++){
tmp222()(isl)() = diquark()(isl,isr)();
}
tmp111 = Gamma(gamma12_[0]) * tmp222;
for (int isl=0 ; isl < 4 ; isl++){
g_diquark()(isl,isr)() = tmp111()(isl)();
}
}
for (int is=0 ; is < 4 ; is++){
BField[Bindex]+=(double)epsilon_sgn[ie]*tmp11s()(is)()*g_diquark()(is,is)();
}
}
}
}
}
}
}
}
//BaryonTensorSet<Complex> BField3(storage,Nmom,4,Nt,N_1,N_2,N_3);
using BaryonTensorSet = Eigen::Tensor<Complex, 6, Eigen::RowMajor>;
BaryonTensorSet BField3(Nmom,4,Nt,N_1,N_2,N_3);
for (int t=0 ; t < Nt ; t++){
Bindex = 0 + N_1*(0 + N_2*(0 + N_3*(0+Nmom*t)));
std::cout << "BaryonField(t=" << t << ") = " << BField[Bindex] << std::endl;
}
} // end if 0
// ONLY THIS IS CORRECT?
std::vector<SpinVector> BField2(Nmom*Nt*N_1*N_2*N_3);
std::vector<SpinVector> BField2(Nmom*Nt*N_1*N_2*N_3);
Complex diquark2;
for (int i1=0 ; i1 < N_1 ; i1++){
for (int i2=0 ; i2 < N_2 ; i2++){
@@ -238,7 +203,10 @@ void TBContraction<FImpl>::execute(void)
tmp222 = g4*tmp111;
tmp111 = 0.5*(double)parity*(tmp111 + tmp222); // P_\pm * ...
diquark2 = factor23[0]*innerProduct(tmp22s,tmp333);
BField2[Bindex]+=(double)epsilon_sgn[ie]*tmp111*diquark2;
//BField2[Bindex]+=(double)epsilon_sgn[ie]*tmp111*diquark2;
for (int is=0 ; is < 4 ; is++){
BField3(imom,is,t,i1,i2,i3)+=(double)epsilon_sgn[ie]*tmp111()(is)()*diquark2;
}
}
}
}
@@ -248,10 +216,17 @@ void TBContraction<FImpl>::execute(void)
}
for (int is=0 ; is < 4 ; is++){
for (int t=0 ; t < Nt ; t++){
Bindex = 0 + N_1*(0 + N_2*(0 + N_3*(0+Nmom*t)));
std::cout << "BaryonField(is=" << is << ",t=" << t << ") = " << BField2[Bindex]()(is)() << std::endl;
// Bindex = 0 + N_1*(0 + N_2*(0 + N_3*(0+Nmom*t)));
// std::cout << "BaryonField(is=" << is << ",t=" << t << ") = " << BField2[Bindex]()(is)() << std::endl;
std::cout << "BaryonField(is=" << is << ",t=" << t << ") = " << BField3(0,is,t,0,0,0) << std::endl;
}
}
//Product ijk * ijk
// for ijk * jik: (4,5),(5,4),(6,6) z.b.
Eigen::array<Eigen::IndexPair<int>, 3> product_dims = { Eigen::IndexPair<int>(4, 4),Eigen::IndexPair<int>(5, 5) ,Eigen::IndexPair<int>(6, 6) };
// Whycan't I choose the dimension to be 3??? Want to sum over them, not save each element!
Eigen::Tensor<Complex,6,Eigen::RowMajor> C2 = BField3.contract(BField3,product_dims);
}
END_MODULE_NAMESPACE

139
tests/IO/Test_serialisation.cc
View File

@@ -91,15 +91,131 @@ void ioTest(const std::string &filename, const O &object, const std::string &nam
R reader(filename);
O buf;
bool good;
#ifndef DEBUG
read(reader, "testobject", buf);
good = (object == buf);
std::cout << name << " IO test: " << (good ? "success" : "failure");
std::cout << std::endl;
bool good = (object == buf);
std::cout << name << " IO test: " << std::endl;
if (!good) exit(EXIT_FAILURE);
#endif
}
#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 std::complex<double> TestScalar;
typedef Eigen::Tensor<iMatrix<TestScalar,1>, 6> TestTensorSingle;
typedef Eigen::Tensor<TestScalar, 3, Eigen::StorageOptions::RowMajor> TestTensor;
typedef Eigen::TensorFixedSize<TestScalar, Eigen::Sizes<9,4,2>, Eigen::StorageOptions::RowMajor> TestTensorFixed;
typedef std::vector<TestTensorFixed> aTestTensorFixed;
typedef Eigen::TensorFixedSize<SpinColourVector, Eigen::Sizes<11,3,2>> LSCTensor;
typedef Eigen::TensorFixedSize<LorentzColourMatrix, Eigen::Sizes<5,7,2>> LCMTensor;
// From Test_serialisation.cc
class ETSerClass: Serializable {
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(ETSerClass
, SpinColourVector, scv
, SpinColourMatrix, scm
, TestTensor, Critter
, TestTensorFixed, FixedCritter
, aTestTensorFixed, aFixedCritter
, LSCTensor, MyLSCTensor
, LCMTensor, MyLCMTensor
);
ETSerClass() : Critter(7,3,2), aFixedCritter(3) {}
};
bool EigenIOTest(void) {
constexpr TestScalar Inc{1,-1};
TestTensorSingle ts(1,1,1,1,1,1);
ts(0,0,0,0,0,0) = Inc * 3.1415927;
ioTest<Hdf5Writer, Hdf5Reader, TestTensorSingle>("iotest_single.h5", ts, "Singlet");
SpinColourVector scv, scv2;
scv2 = scv;
ioTest<Hdf5Writer, Hdf5Reader, SpinColourVector>("iotest_vector.h5", scv, "SpinColourVector");
SpinColourMatrix scm;
ioTest<Hdf5Writer, Hdf5Reader, SpinColourMatrix>("iotest_matrix.h5", scm, "SpinColourMatrix");
TestTensor t(6,3,2);
TestScalar Val{Inc};
for( int i = 0 ; i < t.dimension(0) ; i++)
for( int j = 0 ; j < t.dimension(1) ; j++)
for( int k = 0 ; k < t.dimension(2) ; k++) {
t(i,j,k) = Val;
Val += Inc;
}
ioTest<Hdf5Writer, Hdf5Reader, TestTensor>("iotest_tensor.h5", t, "eigen_tensor_instance_name");
std::cout << "t:";
dump_tensor(t);
// Now serialise a fixed size tensor
using FixedTensor = Eigen::TensorFixedSize<TestScalar, Eigen::Sizes<8,4,3>>;
FixedTensor tf;
Val = Inc;
for( int i = 0 ; i < tf.dimension(0) ; i++)
for( int j = 0 ; j < tf.dimension(1) ; j++)
for( int k = 0 ; k < tf.dimension(2) ; k++) {
tf(i,j,k) = Val;
Val += Inc;
}
ioTest<Hdf5Writer, Hdf5Reader, FixedTensor>("iotest_tensor_fixed.h5", tf, "eigen_tensor_fixed_name");
std::cout << "tf:";
dump_tensor(tf);
ETSerClass o;
ioTest<Hdf5Writer, Hdf5Reader, ETSerClass>("iotest_object.h5", o, "ETSerClass_object_instance_name");
// Tensor of spin colour
LSCTensor l;
Val = 0;
for( int i = 0 ; i < l.dimension(0) ; i++)
for( int j = 0 ; j < l.dimension(1) ; j++)
for( int k = 0 ; k < l.dimension(2) ; k++)
for( int s = 0 ; s < Ns ; s++ )
for( int c = 0 ; c < Nc ; c++ )
{
l(i,j,k)()(s)(c) = Val;
Val += Inc;
}
ioTest<Hdf5Writer, Hdf5Reader, LSCTensor>("iotest_LSCTensor.h5", l, "LSCTensor_object_instance_name");
std::cout << "l:";
dump_tensor(l);
// Tensor of spin colour
LCMTensor l2;
Val = 0;
for( int i = 0 ; i < l2.dimension(0) ; i++)
for( int j = 0 ; j < l2.dimension(1) ; j++)
for( int k = 0 ; k < l2.dimension(2) ; k++)
for( int l = 0 ; l < Ns ; l++ )
for( int c = 0 ; c < Nc ; c++ )
for( int c2 = 0 ; c2 < Nc ; c2++ )
{
l2(i,j,k)(l)()(c,c2) = Val;
Val += Inc;
}
ioTest<Hdf5Writer, Hdf5Reader, LCMTensor>("iotest_LCMTensor.h5", l2, "LCMTensor_object_instance_name");
std::cout << "Wow!" << std::endl;
return true;
}
#endif
template <typename T>
void tensorConvTestFn(GridSerialRNG &rng, const std::string label)
{
@@ -121,12 +237,13 @@ void tensorConvTestFn(GridSerialRNG &rng, const std::string label)
int main(int argc,char **argv)
{
Grid_init(&argc,&argv);
std::cout << std::boolalpha << "==== basic IO" << std::endl; // display true / false for boolean
#ifndef DEBUG
GridSerialRNG rng;
rng.SeedFixedIntegers(std::vector<int>({42,10,81,9}));
std::cout << "==== basic IO" << std::endl;
XmlWriter WR("bother.xml");
// test basic type writing
@@ -160,8 +277,8 @@ int main(int argc,char **argv)
std::cout << "-- serialisable class writing to std::cout:" << std::endl;
std::cout << obj << std::endl;
std::cout << "-- serialisable class comparison:" << std::endl;
std::cout << "vec[0] == obj: " << ((vec[0] == obj) ? "true" : "false") << std::endl;
std::cout << "vec[1] == obj: " << ((vec[1] == obj) ? "true" : "false") << std::endl;
std::cout << "vec[0] == obj: " << (vec[0] == obj) << std::endl;
std::cout << "vec[1] == obj: " << (vec[1] == obj) << std::endl;
std::cout << "-- pair writing to std::cout:" << std::endl;
std::cout << pair << std::endl;
@@ -227,4 +344,8 @@ int main(int argc,char **argv)
tensorConvTest(rng, ColourVector);
tensorConvTest(rng, SpinMatrix);
tensorConvTest(rng, SpinVector);
#elif HAVE_HDF5
if(! EigenIOTest() ) exit(EXIT_FAILURE);
#endif
Grid_finalize();
}

205
tests/hadrons/Test_hadrons_distil.cc
View File

@@ -466,153 +466,62 @@ bool DebugEigenTest()
return true;
}
template <typename W, typename R, typename O>
bool ioTest(const std::string &filename, const O &object, const std::string &name)
{
// writer needs to be destroyed so that writing physically happens
{
W writer(filename);
write(writer, "testobject", object);
}
/*R reader(filename);
O buf;
bool good;
read(reader, "testobject", buf);
good = (object == buf);
std::cout << name << " IO test: " << (good ? "success" : "failure");
std::cout << std::endl;
return good;*/
return true;
}
//typedef int TestScalar;
typedef std::complex<double> TestScalar;
typedef Eigen::Tensor<TestScalar, 3> TestTensor;
typedef Eigen::TensorFixedSize<TestScalar, Eigen::Sizes<9,4,2>> TestTensorFixed;
typedef std::vector<TestTensorFixed> aTestTensorFixed;
typedef Eigen::TensorFixedSize<SpinColourVector, Eigen::Sizes<11,3,2>> LSCTensor;
typedef Eigen::TensorFixedSize<LorentzColourMatrix, Eigen::Sizes<5,7,2>> LCMTensor;
// From Test_serialisation.cc
class myclass: Serializable {
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(myclass
, SpinColourVector, scv
, SpinColourMatrix, scm
, TestTensor, Critter
, TestTensorFixed, FixedCritter
, aTestTensorFixed, aFixedCritter
, LSCTensor, MyLSCTensor
, LCMTensor, MyLCMTensor
);
myclass() : Critter(7,3,2), aFixedCritter(3) {}
};
bool DebugIOTest(void) {
SpinColourVector scv, scv2;
scv2 = scv;
ioTest<Hdf5Writer, Hdf5Reader, SpinColourVector>("iotest_vector.h5", scv, "SpinColourVector");
SpinColourMatrix scm;
ioTest<Hdf5Writer, Hdf5Reader, SpinColourMatrix>("iotest_matrix.h5", scm, "SpinColourMatrix");
constexpr TestScalar Inc{1,-1};
TestTensor t(3,6,2);
TestScalar Val{Inc};
for( int i = 0 ; i < t.dimension(0) ; i++)
for( int j = 0 ; j < t.dimension(1) ; j++)
for( int k = 0 ; k < t.dimension(2) ; k++) {
t(i,j,k) = Val;
Val += Inc;
}
ioTest<Hdf5Writer, Hdf5Reader, TestTensor>("iotest_tensor.h5", t, "eigen_tensor_instance_name");
// Now serialise a fixed size tensor
using FixedTensor = Eigen::TensorFixedSize<TestScalar, Eigen::Sizes<8,4,3>>;
FixedTensor tf;
Val = Inc;
for( int i = 0 ; i < tf.dimension(0) ; i++)
for( int j = 0 ; j < tf.dimension(1) ; j++)
for( int k = 0 ; k < tf.dimension(2) ; k++) {
tf(i,j,k) = Val;
Val += Inc;
}
ioTest<Hdf5Writer, Hdf5Reader, FixedTensor>("iotest_tensor_fixed.h5", tf, "eigen_tensor_fixed_name");
myclass o;
ioTest<Hdf5Writer, Hdf5Reader, myclass>("iotest_object.h5", o, "myclass_object_instance_name");
// Tensor of spin colour
LSCTensor l;
Val = 0;
for( int i = 0 ; i < l.dimension(0) ; i++)
for( int j = 0 ; j < l.dimension(1) ; j++)
for( int k = 0 ; k < l.dimension(2) ; k++)
for( int s = 0 ; s < Ns ; s++ )
for( int c = 0 ; c < Nc ; c++ )
{
l(i,j,k)()(s)(c) = Val;
Val += Inc;
}
ioTest<Hdf5Writer, Hdf5Reader, LSCTensor>("iotest_LSCTensor.h5", l, "LSCTensor_object_instance_name");
// Tensor of spin colour
LCMTensor l2;
Val = 0;
for( int i = 0 ; i < l2.dimension(0) ; i++)
for( int j = 0 ; j < l2.dimension(1) ; j++)
for( int k = 0 ; k < l2.dimension(2) ; k++)
for( int l = 0 ; l < Ns ; l++ )
for( int c = 0 ; c < Nc ; c++ )
for( int c2 = 0 ; c2 < Nc ; c2++ )
{
l2(i,j,k)(l)()(c,c2) = Val;
Val += Inc;
}
ioTest<Hdf5Writer, Hdf5Reader, LCMTensor>("iotest_LCMTensor.h5", l2, "LCMTensor_object_instance_name");
std::cout << "Wow!" << std::endl;
return true;
}
//template <typename T> ReallyIsGridTensor struct {
//false_type;
//}
/*template<typename T>
struct GridSize : public std::integral_constant<std::size_t, 0> {};
template<typename T>
struct GridRank<iScalar<T>> : public std::integral_constant<std::size_t, rank<T>::value + 1> {};
template<class T, std::size_t N>
struct rank<T[N]> : public std::integral_constant<std::size_t, rank<T>::value + 1> {};
*/
template <typename T>
void DebugGridTensorTest_print( int i )
{
std::cout << i << " : " << is_grid_tensor<T>::value
<< ", depth " << grid_tensor_att<T>::depth
<< ", rank " << grid_tensor_att<T>::rank
<< ", rank_non_trivial " << grid_tensor_att<T>::rank_non_trivial
<< ", count " << grid_tensor_att<T>::count
<< ", scalar_size " << grid_tensor_att<T>::scalar_size
<< ", size " << grid_tensor_att<T>::size
std::cout << i << " : " << EigenIO::is_tensor<T>::value
<< ", depth " << EigenIO::Traits<T>::depth
<< ", rank " << EigenIO::Traits<T>::rank
<< ", rank_non_trivial " << EigenIO::Traits<T>::rank_non_trivial
<< ", count " << EigenIO::Traits<T>::count
<< ", scalar_size " << EigenIO::Traits<T>::scalar_size
<< ", size " << EigenIO::Traits<T>::size
<< std::endl;
}
// begin() and end() are the minimum necessary to support range-for loops
// should really turn this into an iterator ...
template<typename T, int N>
class TestObject {
public:
using value_type = T;
private:
value_type * m_p;
public:
TestObject() {
m_p = reinterpret_cast<value_type *>(std::malloc(N * sizeof(value_type)));
}
~TestObject() { std::free(m_p); }
inline value_type * begin(void) { return m_p; }
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 )
{
DebugFelixTensorTest();
typedef Complex t1;
typedef iScalar<t1> t2;
typedef iVector<t1, Ns> t3;
typedef iMatrix<t1, Nc> t4;
typedef iVector<iMatrix<t1,1>,4> t5;
typedef iScalar<t5> t6;
typedef iMatrix<iVector<iScalar<iMatrix<t6, 1>>,2>,7> t7;
typedef iMatrix<t6, 3> t7;
typedef iMatrix<iVector<iScalar<t7>,4>,2> t8;
int i = 1;
DebugGridTensorTest_print<t1>( i++ );
DebugGridTensorTest_print<t2>( i++ );
@@ -621,6 +530,27 @@ bool DebugGridTensorTest( void )
DebugGridTensorTest_print<t5>( i++ );
DebugGridTensorTest_print<t6>( i++ );
DebugGridTensorTest_print<t7>( i++ );
DebugGridTensorTest_print<t8>( i++ );
//using TOC7 = TestObject<std::complex<double>, 7>;
using TOC7 = t7;
TOC7 toc7;
constexpr std::complex<double> Inc{1,-1};
std::complex<double> Start{Inc};
for( auto &x : toc7 ) {
x = Start;
Start += Inc;
}
i = 0;
for( auto x : toc7 ) std::cout << "toc7[" << i++ << "] = " << x << std::endl;
t2 o2;
auto a2 = TensorRemove(o2);
//t3 o3;
//t4 o4;
//auto a3 = TensorRemove(o3);
//auto a4 = TensorRemove(o4);
return true;
}
#endif
@@ -629,11 +559,14 @@ int main(int argc, char *argv[])
{
#ifdef DEBUG
// Debug only - test of Eigen::Tensor
std::cout << "sizeof(std::streamsize) = " << sizeof(std::streamsize) << std::endl;
std::cout << "sizeof(Eigen::Index) = " << sizeof(Eigen::Index) << std::endl;
std::cout << "sizeof(int) = " << sizeof(int)
<< ", sizeof(long) = " << sizeof(long)
<< ", sizeof(size_t) = " << sizeof(size_t)
<< ", sizeof(std::size_t) = " << sizeof(std::size_t)
<< ", sizeof(std::streamsize) = " << sizeof(std::streamsize)
<< ", sizeof(Eigen::Index) = " << sizeof(Eigen::Index) << std::endl;
//if( DebugEigenTest() ) return 0;
if(DebugGridTensorTest()) return 0;
//if(DebugIOTest()) return 0;
#endif
// Decode command-line parameters. 1st one is which test to run