diff --git a/.gitignore b/.gitignore index 45a3ea53..5338acb9 100644 --- a/.gitignore +++ b/.gitignore @@ -114,3 +114,4 @@ gh-pages/ ##################### Grid/qcd/spin/gamma-gen/*.h Grid/qcd/spin/gamma-gen/*.cc +Grid/util/Version.h diff --git a/Grid/Grid.h b/Grid/Grid.h index 9dcc207b..70859c81 100644 --- a/Grid/Grid.h +++ b/Grid/Grid.h @@ -42,6 +42,7 @@ Author: paboyle #include #include #include +#include #include #include #include diff --git a/Grid/parallelIO/BinaryIO.h b/Grid/parallelIO/BinaryIO.h index 1895dc3e..144ff29f 100644 --- a/Grid/parallelIO/BinaryIO.h +++ b/Grid/parallelIO/BinaryIO.h @@ -619,6 +619,7 @@ PARALLEL_CRITICAL { std::cout << GridLogMessage << "writeLatticeObject: read test checksum failure, re-writing (" << attemptsLeft << " attempt(s) remaining)" << std::endl; offset = offsetCopy; + parallel_for(uint64_t x=0;x_fdimensions[nu+shift]))); } - in_buf = exp((Real)(2.0*M_PI)*ci*ph*(-1.0))*in; + in_buf = exp(Scalar(2.0*M_PI)*ci*ph*(-1.0))*in; if(fiveD){//FFT only on temporal and spatial dimensions std::vector mask(Nd+1,1); mask[0] = 0; @@ -78,7 +78,7 @@ namespace Grid { } //phase for boundary condition - out = out * exp((Real)(2.0*M_PI)*ci*ph); + out = out * exp(Scalar(2.0*M_PI)*ci*ph); }; virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector twist) { diff --git a/Grid/qcd/action/fermion/FermionOperator.h b/Grid/qcd/action/fermion/FermionOperator.h index 55942dea..436808e6 100644 --- a/Grid/qcd/action/fermion/FermionOperator.h +++ b/Grid/qcd/action/fermion/FermionOperator.h @@ -106,21 +106,21 @@ namespace Grid { ComplexField coor(in._grid); ComplexField ph(in._grid); ph = zero; FermionField in_buf(in._grid); in_buf = zero; - Complex ci(0.0,1.0); + Scalar ci(0.0,1.0); assert(twist.size() == Nd);//check that twist is Nd for(unsigned int nu = 0; nu < Nd; nu++) { LatticeCoordinate(coor, nu); ph = ph + twist[nu]*coor*((1./(in._grid->_fdimensions[nu]))); } - in_buf = exp((Real)(2.0*M_PI)*ci*ph*(-1.0))*in; + in_buf = exp(Scalar(-2.0*M_PI)*ci*ph)*in; theFFT.FFT_all_dim(in_k,in_buf,FFT::forward); this->MomentumSpacePropagator(prop_k,in_k,mass,twist); theFFT.FFT_all_dim(out,prop_k,FFT::backward); //phase for boundary condition - out = out * exp((Real)(2.0*M_PI)*ci*ph); + out = out * exp(Scalar(2.0*M_PI)*ci*ph); }; virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h index 95680900..52e1e100 100644 --- a/Grid/qcd/action/gauge/GaugeImplTypes.h +++ b/Grid/qcd/action/gauge/GaugeImplTypes.h @@ -46,6 +46,7 @@ namespace QCD { #define INHERIT_GIMPL_TYPES(GImpl) \ typedef typename GImpl::Simd Simd; \ + typedef typename GImpl::Scalar Scalar; \ typedef typename GImpl::LinkField GaugeLinkField; \ typedef typename GImpl::Field GaugeField; \ typedef typename GImpl::ComplexField ComplexField;\ @@ -63,7 +64,8 @@ namespace QCD { template class GaugeImplTypes { public: typedef S Simd; - + typedef typename Simd::scalar_type scalar_type; + typedef scalar_type Scalar; template using iImplScalar = iScalar > >; template using iImplGaugeLink = iScalar > >; template using iImplGaugeField = iVector >, Nd>; diff --git a/Grid/qcd/action/gauge/Photon.h b/Grid/qcd/action/gauge/Photon.h index f059fcf3..9afafe6c 100644 --- a/Grid/qcd/action/gauge/Photon.h +++ b/Grid/qcd/action/gauge/Photon.h @@ -38,6 +38,7 @@ namespace QCD{ { public: typedef S Simd; + typedef typename Simd::scalar_type Scalar; template using iImplGaugeLink = iScalar>>; diff --git a/Grid/qcd/utils/CovariantSmearing.h b/Grid/qcd/utils/CovariantSmearing.h new file mode 100644 index 00000000..eae5ed71 --- /dev/null +++ b/Grid/qcd/utils/CovariantSmearing.h @@ -0,0 +1,87 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/scalar/CovariantLaplacian.h + +Copyright (C) 2016 + +Author: Azusa Yamaguchi + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +#pragma once + +namespace Grid { +namespace QCD { + +template class CovariantSmearing : public Gimpl +{ +public: + INHERIT_GIMPL_TYPES(Gimpl); + + typedef typename Gimpl::GaugeLinkField GaugeMat; + typedef typename Gimpl::GaugeField GaugeLorentz; + + template + static void GaussianSmear(const std::vector& U, + T& chi, + const Real& width, int Iterations, int orthog) + { + GridBase *grid = chi._grid; + T psi(grid); + + //////////////////////////////////////////////////////////////////////////////////// + // Follow Chroma conventions for width to keep compatibility with previous data + // Free field iterates + // chi = (1 - w^2/4N p^2)^N chi + // + // ~ (e^(-w^2/4N p^2)^N chi + // ~ (e^(-w^2/4 p^2) chi + // ~ (e^(-w'^2/2 p^2) chi [ w' = w/sqrt(2) ] + // + // Which in coordinate space is proportional to + // + // e^(-x^2/w^2) = e^(-x^2/2w'^2) + // + // The 4 is a bit unconventional from Gaussian width perspective, but... it's Chroma convention. + // 2nd derivative approx d^2/dx^2 = x+mu + x-mu - 2x + // + // d^2/dx^2 = - p^2 + // + // chi = ( 1 + w^2/4N d^2/dx^2 )^N chi + // + //////////////////////////////////////////////////////////////////////////////////// + Real coeff = (width*width) / Real(4*Iterations); + + int dims = Nd; + if( orthog < Nd ) dims=Nd-1; + + for(int n = 0; n < Iterations; ++n) { + psi = (-2.0*dims)*chi; + for(int mu=0;mu #include #include #include +#include namespace Grid { + namespace EigenIO { + // EigenIO works for scalars that are not just Grid supported scalars + template struct is_complex : public std::false_type {}; + // Support all complex types (not just Grid complex types) - even if the definitions overlap (!) + template struct is_complex< T , typename + std::enable_if< ::Grid::is_complex< T >::value>::type> : public std::true_type {}; + template struct is_complex, typename + std::enable_if>::value>::type> : public std::true_type {}; + + // Helpers to support I/O for Eigen tensors of arithmetic scalars, complex types, or Grid tensors + template struct is_scalar : public std::false_type {}; + template struct is_scalar::value || is_complex::value>::type> : public std::true_type {}; + + // Is this an Eigen tensor + template struct is_tensor : std::integral_constant, T>::value> {}; + + // Is this an Eigen tensor of a supported scalar + template struct is_tensor_of_scalar : public std::false_type {}; + template struct is_tensor_of_scalar::value && is_scalar::value>::type> : public std::true_type {}; + + // Is this an Eigen tensor of a supported container + template struct is_tensor_of_container : public std::false_type {}; + template struct is_tensor_of_container::value && isGridTensor::value>::type> : public std::true_type {}; + + // These traits describe the scalars inside Eigen tensors + // I wish I could define these in reference to the scalar type (so there would be fewer traits defined) + // but I'm unable to find a syntax to make this work + template struct Traits {}; + // Traits are the default for scalars, or come from GridTypeMapper for GridTensors + template struct Traits::value>::type> + : public GridTypeMapper_Base { + using scalar_type = typename T::Scalar; // ultimate base scalar + static constexpr bool is_complex = ::Grid::EigenIO::is_complex::value; + }; + // Traits are the default for scalars, or come from GridTypeMapper for GridTensors + template struct Traits::value>::type> { + using BaseTraits = GridTypeMapper; + using scalar_type = typename BaseTraits::scalar_type; // ultimate base scalar + static constexpr bool is_complex = ::Grid::EigenIO::is_complex::value; + static constexpr int TensorLevel = BaseTraits::TensorLevel; + static constexpr int Rank = BaseTraits::Rank; + static constexpr std::size_t count = BaseTraits::count; + static constexpr int Dimension(int dim) { return BaseTraits::Dimension(dim); } + }; + + // Is this a fixed-size Eigen tensor + template struct is_tensor_fixed : public std::false_type {}; + template + struct is_tensor_fixed> + : public std::true_type {}; + template class MapPointer_> + struct is_tensor_fixed, MapOptions_, MapPointer_>> + : public std::true_type {}; + + // Is this a variable-size Eigen tensor + template struct is_tensor_variable : public std::false_type {}; + template struct is_tensor_variable::value + && !is_tensor_fixed::value>::type> : public std::true_type {}; + } + // Abstract writer/reader classes //////////////////////////////////////////// // static polymorphism implemented using CRTP idiom class Serializable; - + // Static abstract writer template class Writer @@ -49,10 +113,10 @@ namespace Grid { void push(const std::string &s); void pop(void); template - typename std::enable_if::value, void>::type + typename std::enable_if::value>::type write(const std::string& s, const U &output); template - typename std::enable_if::value, void>::type + typename std::enable_if::value && !EigenIO::is_tensor::value>::type write(const std::string& s, const U &output); template void write(const std::string &s, const iScalar &output); @@ -60,6 +124,42 @@ namespace Grid { void write(const std::string &s, const iVector &output); template void write(const std::string &s, const iMatrix &output); + template + typename std::enable_if::value>::type + write(const std::string &s, const ETensor &output); + + // Helper functions for Scalar vs Container specialisations + template + inline typename std::enable_if::value, + const typename ETensor::Scalar *>::type + getFirstScalar(const ETensor &output) + { + return output.data(); + } + + template + inline typename std::enable_if::value, + const typename EigenIO::Traits::scalar_type *>::type + getFirstScalar(const ETensor &output) + { + return output.data()->begin(); + } + + template + inline typename std::enable_if::value, void>::type + copyScalars(S * &pCopy, const S &Source) + { + * pCopy ++ = Source; + } + + template + inline typename std::enable_if::value, void>::type + copyScalars(typename GridTypeMapper::scalar_type * &pCopy, const S &Source) + { + for( const typename GridTypeMapper::scalar_type &item : Source ) + * pCopy ++ = item; + } + void scientificFormat(const bool set); bool isScientific(void); void setPrecision(const unsigned int prec); @@ -83,7 +183,8 @@ namespace Grid { typename std::enable_if::value, void>::type read(const std::string& s, U &output); template - typename std::enable_if::value, void>::type + typename std::enable_if::value + && !EigenIO::is_tensor::value, void>::type read(const std::string& s, U &output); template void read(const std::string &s, iScalar &output); @@ -91,6 +192,32 @@ namespace Grid { void read(const std::string &s, iVector &output); template void read(const std::string &s, iMatrix &output); + template + typename std::enable_if::value, void>::type + read(const std::string &s, ETensor &output); + template + typename std::enable_if::value, void>::type + Reshape(ETensor &t, const std::array &dims ); + template + typename std::enable_if::value, void>::type + Reshape(ETensor &t, const std::array &dims ); + + // Helper functions for Scalar vs Container specialisations + template + inline typename std::enable_if::value, void>::type + copyScalars(S &Dest, const S * &pSource) + { + Dest = * pSource ++; + } + + template + inline typename std::enable_if::value, void>::type + copyScalars(S &Dest, const typename GridTypeMapper::scalar_type * &pSource) + { + for( typename GridTypeMapper::scalar_type &item : Dest ) + item = * pSource ++; + } + protected: template void fromString(U &output, const std::string &s); @@ -135,12 +262,14 @@ namespace Grid { template template - typename std::enable_if::value, void>::type + typename std::enable_if::value + && !EigenIO::is_tensor::value, void>::type Writer::write(const std::string &s, const U &output) { upcast->writeDefault(s, output); } + template template void Writer::write(const std::string &s, const iScalar &output) @@ -161,6 +290,57 @@ namespace Grid { { upcast->writeDefault(s, tensorToVec(output)); } + + // Eigen::Tensors of Grid tensors (iScalar, iVector, iMatrix) + template + template + typename std::enable_if::value, void>::type + Writer::write(const std::string &s, const ETensor &output) + { + using Index = typename ETensor::Index; + using Container = typename ETensor::Scalar; // NB: could be same as scalar + using Traits = EigenIO::Traits; + using Scalar = typename Traits::scalar_type; // type of the underlying scalar + constexpr unsigned int TensorRank{ETensor::NumIndices}; + constexpr unsigned int ContainerRank{Traits::Rank}; // Only non-zero for containers + constexpr unsigned int TotalRank{TensorRank + ContainerRank}; + const Index NumElements{output.size()}; + assert( NumElements > 0 ); + + // Get the dimensionality of the tensor + std::vector TotalDims(TotalRank); + for(auto i = 0; i < TensorRank; i++ ) { + auto dim = output.dimension(i); + TotalDims[i] = static_cast(dim); + assert( TotalDims[i] == dim ); // check we didn't lose anything in the conversion + } + for(auto i = 0; i < ContainerRank; i++ ) + TotalDims[TensorRank + i] = Traits::Dimension(i); + + // If the Tensor isn't in Row-Major order, then we'll need to copy it's data + const bool CopyData{NumElements > 1 && ETensor::Layout != Eigen::StorageOptions::RowMajor}; + const Scalar * pWriteBuffer; + std::vector CopyBuffer; + const Index TotalNumElements = NumElements * Traits::count; + if( !CopyData ) { + pWriteBuffer = getFirstScalar( output ); + } else { + // Regardless of the Eigen::Tensor storage order, the copy will be Row Major + CopyBuffer.resize( TotalNumElements ); + Scalar * pCopy = &CopyBuffer[0]; + pWriteBuffer = pCopy; + std::array MyIndex; + for( auto &idx : MyIndex ) idx = 0; + for( auto n = 0; n < NumElements; n++ ) { + const Container & c = output( MyIndex ); + copyScalars( pCopy, c ); + // Now increment the index + for( int i = output.NumDimensions - 1; i >= 0 && ++MyIndex[i] == output.dimension(i); i-- ) + MyIndex[i] = 0; + } + } + upcast->template writeMultiDim(s, TotalDims, pWriteBuffer, TotalNumElements); + } template void Writer::scientificFormat(const bool set) @@ -215,7 +395,8 @@ namespace Grid { template template - typename std::enable_if::value, void>::type + typename std::enable_if::value + && !EigenIO::is_tensor::value, void>::type Reader::read(const std::string &s, U &output) { upcast->readDefault(s, output); @@ -251,6 +432,79 @@ namespace Grid { vecToTensor(output, v); } + template + template + typename std::enable_if::value, void>::type + Reader::read(const std::string &s, ETensor &output) + { + using Index = typename ETensor::Index; + using Container = typename ETensor::Scalar; // NB: could be same as scalar + using Traits = EigenIO::Traits; + using Scalar = typename Traits::scalar_type; // type of the underlying scalar + constexpr unsigned int TensorRank{ETensor::NumIndices}; + constexpr unsigned int ContainerRank{Traits::Rank}; // Only non-zero for containers + constexpr unsigned int TotalRank{TensorRank + ContainerRank}; + using ETDims = std::array; // Dimensions of the tensor + + // read the (flat) data and dimensionality + std::vector dimData; + std::vector buf; + upcast->readMultiDim( s, buf, dimData ); + assert(dimData.size() == TotalRank && "EigenIO: Tensor rank mismatch" ); + // Make sure that the number of elements read matches dimensions read + std::size_t NumContainers = 1; + for( auto i = 0 ; i < TensorRank ; i++ ) + NumContainers *= dimData[i]; + // If our scalar object is a Container, make sure it's dimensions match what we read back + std::size_t ElementsPerContainer = 1; + for( auto i = 0 ; i < ContainerRank ; i++ ) { + assert( dimData[TensorRank+i] == Traits::Dimension(i) && "Tensor Container dimensions don't match data" ); + ElementsPerContainer *= dimData[TensorRank+i]; + } + assert( NumContainers * ElementsPerContainer == buf.size() && "EigenIO: Number of elements != product of dimensions" ); + // Now see whether the tensor is the right shape, or can be made to be + const auto & dims = output.dimensions(); + bool bShapeOK = (output.data() != nullptr); + for( auto i = 0; bShapeOK && i < TensorRank ; i++ ) + if( dims[i] != dimData[i] ) + bShapeOK = false; + // Make the tensor the same size as the data read + ETDims MyIndex; + if( !bShapeOK ) { + for( auto i = 0 ; i < TensorRank ; i++ ) + MyIndex[i] = dimData[i]; + Reshape(output, MyIndex); + } + // Copy the data into the tensor + for( auto &d : MyIndex ) d = 0; + const Scalar * pSource = &buf[0]; + for( std::size_t n = 0 ; n < NumContainers ; n++ ) { + Container & c = output( MyIndex ); + copyScalars( c, pSource ); + // Now increment the index + for( int i = TensorRank - 1; i != -1 && ++MyIndex[i] == dims[i]; i-- ) + MyIndex[i] = 0; + } + assert( pSource == &buf[NumContainers * ElementsPerContainer] ); + } + + template + template + typename std::enable_if::value, void>::type + Reader::Reshape(ETensor &t, const std::array &dims ) + { + assert( 0 && "EigenIO: Fixed tensor dimensions can't be changed" ); + } + + template + template + typename std::enable_if::value, void>::type + Reader::Reshape(ETensor &t, const std::array &dims ) + { + //t.reshape( dims ); + t.resize( dims ); + } + template template void Reader::fromString(U &output, const std::string &s) @@ -289,8 +543,70 @@ namespace Grid { { return os; } + + template + static inline typename std::enable_if::value || !EigenIO::is_tensor::value, bool>::type + CompareMember(const T1 &lhs, const T2 &rhs) { + return lhs == rhs; + } + + template + static inline typename std::enable_if::value && EigenIO::is_tensor::value, bool>::type + CompareMember(const T1 &lhs, const T2 &rhs) { + // First check whether dimensions match (Eigen tensor library will assert if they don't match) + bool bReturnValue = (T1::NumIndices == T2::NumIndices); + for( auto i = 0 ; bReturnValue && i < T1::NumIndices ; i++ ) + bReturnValue = ( lhs.dimension(i) == rhs.dimension(i) ); + if( bReturnValue ) { + Eigen::Tensor bResult = (lhs == rhs).all(); + bReturnValue = bResult(0); + } + return bReturnValue; + } + + template + static inline typename std::enable_if::value, bool>::type + CompareMember(const std::vector &lhs, const std::vector &rhs) { + const auto NumElements = lhs.size(); + bool bResult = ( NumElements == rhs.size() ); + for( auto i = 0 ; i < NumElements && bResult ; i++ ) + bResult = CompareMember(lhs[i], rhs[i]); + return bResult; + } + + template + static inline typename std::enable_if::value, void>::type + WriteMember(std::ostream &os, const T &object) { + os << object; + } + + template + static inline typename std::enable_if::value, void>::type + WriteMember(std::ostream &os, const T &object) { + using Index = typename T::Index; + const Index NumElements{object.size()}; + assert( NumElements > 0 ); + Index count = 1; + os << "T<"; + for( int i = 0; i < T::NumIndices; i++ ) { + Index dim = object.dimension(i); + count *= dim; + if( i ) + os << ","; + os << dim; + } + assert( count == NumElements && "Number of elements doesn't match tensor dimensions" ); + os << ">{"; + const typename T::Scalar * p = object.data(); + for( Index i = 0; i < count; i++ ) { + if( i ) + os << ","; + os << *p++; + } + os << "}"; + } }; - + // Generic writer interface ////////////////////////////////////////////////// template inline void push(Writer &w, const std::string &s) { diff --git a/Grid/serialisation/BinaryIO.h b/Grid/serialisation/BinaryIO.h index 757753c7..bf011454 100644 --- a/Grid/serialisation/BinaryIO.h +++ b/Grid/serialisation/BinaryIO.h @@ -51,6 +51,8 @@ namespace Grid { template void writeDefault(const std::string &s, const std::vector &x); void writeDefault(const std::string &s, const char *x); + template + void writeMultiDim(const std::string &s, const std::vector & Dimensions, const U * pDataRowMajor, size_t NumElements); private: std::ofstream file_; }; @@ -66,6 +68,8 @@ namespace Grid { void readDefault(const std::string &s, U &output); template void readDefault(const std::string &s, std::vector &output); + template + void readMultiDim(const std::string &s, std::vector &buf, std::vector &dim); private: std::ifstream file_; }; @@ -92,6 +96,27 @@ namespace Grid { } } + template + void BinaryWriter::writeMultiDim(const std::string &s, const std::vector & Dimensions, const U * pDataRowMajor, size_t NumElements) + { + uint64_t rank = static_cast( Dimensions.size() ); + uint64_t tmp = 1; + for( auto i = 0 ; i < rank ; i++ ) + tmp *= Dimensions[i]; + assert( tmp == NumElements && "Dimensions don't match size of data being written" ); + // Total number of elements + write("", tmp); + // Number of dimensions + write("", rank); + // Followed by each dimension + for( auto i = 0 ; i < rank ; i++ ) { + tmp = Dimensions[i]; + write("", tmp); + } + for( auto i = 0; i < NumElements; ++i) + write("", pDataRowMajor[i]); + } + // Reader template implementation //////////////////////////////////////////// template void BinaryReader::readDefault(const std::string &s, U &output) @@ -114,6 +139,30 @@ namespace Grid { read("", output[i]); } } + + template + void BinaryReader::readMultiDim(const std::string &s, std::vector &buf, std::vector &dim) + { + // Number of elements + uint64_t NumElements; + read("", NumElements); + // Number of dimensions + uint64_t rank; + read("", rank); + // Followed by each dimension + uint64_t count = 1; + dim.resize(rank); + uint64_t tmp; + for( auto i = 0 ; i < rank ; i++ ) { + read("", tmp); + dim[i] = tmp; + count *= tmp; + } + assert( count == NumElements && "Dimensions don't match size of data being read" ); + buf.resize(count); + for( auto i = 0; i < count; ++i) + read("", buf[i]); + } } #endif diff --git a/Grid/serialisation/Hdf5IO.h b/Grid/serialisation/Hdf5IO.h index 59804240..19537599 100644 --- a/Grid/serialisation/Hdf5IO.h +++ b/Grid/serialisation/Hdf5IO.h @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -38,6 +39,8 @@ namespace Grid template typename std::enable_if>::is_number, void>::type writeDefault(const std::string &s, const std::vector &x); + template + void writeMultiDim(const std::string &s, const std::vector & Dimensions, const U * pDataRowMajor, size_t NumElements); H5NS::Group & getGroup(void); private: template @@ -48,7 +51,7 @@ namespace Grid std::vector path_; H5NS::H5File file_; H5NS::Group group_; - unsigned int dataSetThres_{HDF5_DEF_DATASET_THRES}; + const unsigned int dataSetThres_{HDF5_DEF_DATASET_THRES}; }; class Hdf5Reader: public Reader @@ -66,6 +69,8 @@ namespace Grid template typename std::enable_if>::is_number, void>::type readDefault(const std::string &s, std::vector &x); + template + void readMultiDim(const std::string &s, std::vector &buf, std::vector &dim); H5NS::Group & getGroup(void); private: template @@ -101,6 +106,75 @@ namespace Grid template <> void Hdf5Writer::writeDefault(const std::string &s, const std::string &x); + template + void Hdf5Writer::writeMultiDim(const std::string &s, const std::vector & Dimensions, const U * pDataRowMajor, size_t NumElements) + { + // Hdf5 needs the dimensions as hsize_t + const int rank = static_cast(Dimensions.size()); + std::vector dim(rank); + for(int i = 0; i < rank; i++) + dim[i] = Dimensions[i]; + // write the entire dataset to file + H5NS::DataSpace dataSpace(rank, dim.data()); + + if (NumElements > dataSetThres_) + { + // Make sure 1) each dimension; and 2) chunk size is < 4GB + const hsize_t MaxElements = ( sizeof( U ) == 1 ) ? 0xffffffff : 0x100000000 / sizeof( U ); + hsize_t ElementsPerChunk = 1; + bool bTooBig = false; + for( int i = rank - 1 ; i != -1 ; i-- ) { + auto &d = dim[i]; + if( bTooBig ) + d = 1; // Chunk size is already as big as can be - remaining dimensions = 1 + else { + // If individual dimension too big, reduce by prime factors if possible + while( d > MaxElements && ( d & 1 ) == 0 ) + d >>= 1; + const char ErrorMsg[] = " dimension > 4GB and not divisible by 2^n. " + "Hdf5IO chunk size will be inefficient. NB Serialisation is not intended for large datasets - please consider alternatives."; + if( d > MaxElements ) { + std::cout << GridLogWarning << "Individual" << ErrorMsg << std::endl; + hsize_t quotient = d / MaxElements; + if( d % MaxElements ) + quotient++; + d /= quotient; + } + // Now make sure overall size is not too big + hsize_t OverflowCheck = ElementsPerChunk; + ElementsPerChunk *= d; + assert( OverflowCheck == ElementsPerChunk / d && "Product of dimensions overflowed hsize_t" ); + // If product of dimensions too big, reduce by prime factors + while( ElementsPerChunk > MaxElements && ( ElementsPerChunk & 1 ) == 0 ) { + bTooBig = true; + d >>= 1; + ElementsPerChunk >>= 1; + } + if( ElementsPerChunk > MaxElements ) { + std::cout << GridLogWarning << "Product of" << ErrorMsg << std::endl; + hsize_t quotient = ElementsPerChunk / MaxElements; + if( ElementsPerChunk % MaxElements ) + quotient++; + d /= quotient; + ElementsPerChunk /= quotient; + } + } + } + H5NS::DataSet dataSet; + H5NS::DSetCreatPropList plist; + plist.setChunk(rank, dim.data()); + plist.setFletcher32(); + dataSet = group_.createDataSet(s, Hdf5Type::type(), dataSpace, plist); + dataSet.write(pDataRowMajor, Hdf5Type::type()); + } + else + { + H5NS::Attribute attribute; + attribute = group_.createAttribute(s, Hdf5Type::type(), dataSpace); + attribute.write(Hdf5Type::type(), pDataRowMajor); + } + } + template typename std::enable_if>::is_number, void>::type Hdf5Writer::writeDefault(const std::string &s, const std::vector &x) @@ -110,34 +184,11 @@ namespace Grid // flatten the vector and getting dimensions Flatten> flat(x); - std::vector dim; + std::vector dim; const auto &flatx = flat.getFlatVector(); - for (auto &d: flat.getDim()) - { dim.push_back(d); - } - - // write to file - H5NS::DataSpace dataSpace(dim.size(), dim.data()); - - if (flatx.size() > dataSetThres_) - { - H5NS::DataSet dataSet; - H5NS::DSetCreatPropList plist; - - plist.setChunk(dim.size(), dim.data()); - plist.setFletcher32(); - dataSet = group_.createDataSet(s, Hdf5Type::type(), dataSpace, plist); - dataSet.write(flatx.data(), Hdf5Type::type()); - } - else - { - H5NS::Attribute attribute; - - attribute = group_.createAttribute(s, Hdf5Type::type(), dataSpace); - attribute.write(Hdf5Type::type(), flatx.data()); - } + writeMultiDim(s, dim, &flatx[0], flatx.size()); } template @@ -173,10 +224,9 @@ namespace Grid template <> void Hdf5Reader::readDefault(const std::string &s, std::string &x); - + template - typename std::enable_if>::is_number, void>::type - Hdf5Reader::readDefault(const std::string &s, std::vector &x) + void Hdf5Reader::readMultiDim(const std::string &s, std::vector &buf, std::vector &dim) { // alias to element type typedef typename element>::type Element; @@ -184,7 +234,6 @@ namespace Grid // read the dimensions H5NS::DataSpace dataSpace; std::vector hdim; - std::vector dim; hsize_t size = 1; if (group_.attrExists(s)) @@ -204,8 +253,8 @@ namespace Grid } // read the flat vector - std::vector buf(size); - + buf.resize(size); + if (size > dataSetThres_) { H5NS::DataSet dataSet; @@ -220,7 +269,19 @@ namespace Grid attribute = group_.openAttribute(s); attribute.read(Hdf5Type::type(), buf.data()); } - + } + + template + typename std::enable_if>::is_number, void>::type + Hdf5Reader::readDefault(const std::string &s, std::vector &x) + { + // alias to element type + typedef typename element>::type Element; + + std::vector dim; + std::vector buf; + readMultiDim( s, buf, dim ); + // reconstruct the multidimensional vector Reconstruct> r(buf, dim); diff --git a/Grid/serialisation/MacroMagic.h b/Grid/serialisation/MacroMagic.h index 95cd1c3c..ce305673 100644 --- a/Grid/serialisation/MacroMagic.h +++ b/Grid/serialisation/MacroMagic.h @@ -109,8 +109,8 @@ THE SOFTWARE. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #define GRID_MACRO_MEMBER(A,B) A B; -#define GRID_MACRO_COMP_MEMBER(A,B) result = (result and (lhs. B == rhs. B)); -#define GRID_MACRO_OS_WRITE_MEMBER(A,B) os<< #A <<" " #B << " = " << obj. B << " ; " < void writeDefault(const std::string &s, const std::vector &x); + template + void writeMultiDim(const std::string &s, const std::vector & Dimensions, const U * pDataRowMajor, size_t NumElements); private: void indent(void); private: @@ -69,6 +71,8 @@ namespace Grid void readDefault(const std::string &s, U &output); template void readDefault(const std::string &s, std::vector &output); + template + void readMultiDim(const std::string &s, std::vector &buf, std::vector &dim); private: void checkIndent(void); private: @@ -95,7 +99,18 @@ namespace Grid write(s, x[i]); } } - + + template + void TextWriter::writeMultiDim(const std::string &s, const std::vector & Dimensions, const U * pDataRowMajor, size_t NumElements) + { + uint64_t Rank = Dimensions.size(); + write(s, Rank); + for( uint64_t d : Dimensions ) + write(s, d); + while( NumElements-- ) + write(s, *pDataRowMajor++); + } + // Reader template implementation //////////////////////////////////////////// template void TextReader::readDefault(const std::string &s, U &output) @@ -121,6 +136,23 @@ namespace Grid read("", output[i]); } } + + template + void TextReader::readMultiDim(const std::string &s, std::vector &buf, std::vector &dim) + { + const char sz[] = ""; + uint64_t Rank; + read(sz, Rank); + dim.resize( Rank ); + size_t NumElements = 1; + for( auto &d : dim ) { + read(sz, d); + NumElements *= d; + } + buf.resize( NumElements ); + for( auto &x : buf ) + read(s, x); + } } #endif diff --git a/Grid/serialisation/VectorUtils.h b/Grid/serialisation/VectorUtils.h index b6b95c10..a5a73992 100644 --- a/Grid/serialisation/VectorUtils.h +++ b/Grid/serialisation/VectorUtils.h @@ -1,3 +1,32 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./Grid/serialisation/VectorUtils.h + + Copyright (C) 2015 + + Author: Antonin Portelli + Author: Peter Boyle + Author: paboyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ +/* END LEGAL */ #ifndef GRID_SERIALISATION_VECTORUTILS_H #define GRID_SERIALISATION_VECTORUTILS_H @@ -53,6 +82,17 @@ namespace Grid { return os; } + // std::vector> nested to specified Rank ////////////////////////////////// + template + struct NestedStdVector { + typedef typename std::vector::type> type; + }; + + template + struct NestedStdVector { + typedef T type; + }; + // Grid scalar tensors to nested std::vectors ////////////////////////////////// template struct TensorToVec @@ -436,4 +476,4 @@ std::string vecToStr(const std::vector &v) return sstr.str(); } -#endif \ No newline at end of file +#endif diff --git a/Grid/serialisation/XmlIO.h b/Grid/serialisation/XmlIO.h index a26a33c5..40d5f2bb 100644 --- a/Grid/serialisation/XmlIO.h +++ b/Grid/serialisation/XmlIO.h @@ -57,6 +57,8 @@ namespace Grid void writeDefault(const std::string &s, const U &x); template void writeDefault(const std::string &s, const std::vector &x); + template + void writeMultiDim(const std::string &s, const std::vector & Dimensions, const U * pDataRowMajor, size_t NumElements); std::string docString(void); std::string string(void); private: @@ -79,6 +81,8 @@ namespace Grid void readDefault(const std::string &s, U &output); template void readDefault(const std::string &s, std::vector &output); + template + void readMultiDim(const std::string &s, std::vector &buf, std::vector &dim); void readCurrentSubtree(std::string &s); private: void checkParse(const pugi::xml_parse_result &result, const std::string name); @@ -122,13 +126,45 @@ namespace Grid void XmlWriter::writeDefault(const std::string &s, const std::vector &x) { push(s); - for (auto &x_i: x) + for( auto &u : x ) { - write("elem", x_i); + write("elem", u); } pop(); } - + + template + void XmlWriter::writeMultiDim(const std::string &s, const std::vector & Dimensions, const U * pDataRowMajor, size_t NumElements) + { + push(s); + size_t count = 1; + const int Rank = static_cast( Dimensions.size() ); + write("rank", Rank ); + std::vector MyIndex( Rank ); + for( auto d : Dimensions ) { + write("dim", d); + count *= d; + } + assert( count == NumElements && "XmlIO : element count doesn't match dimensions" ); + static const char sName[] = "tensor"; + for( int i = 0 ; i < Rank ; i++ ) { + MyIndex[i] = 0; + push(sName); + } + while (NumElements--) { + write("elem", *pDataRowMajor++); + int i; + for( i = Rank - 1 ; i != -1 && ++MyIndex[i] == Dimensions[i] ; i-- ) + MyIndex[i] = 0; + int Rollover = Rank - 1 - i; + for( i = 0 ; i < Rollover ; i++ ) + pop(); + for( i = 0 ; NumElements && i < Rollover ; i++ ) + push(sName); + } + pop(); + } + // Reader template implementation //////////////////////////////////////////// template void XmlReader::readDefault(const std::string &s, U &output) @@ -145,25 +181,66 @@ namespace Grid template void XmlReader::readDefault(const std::string &s, std::vector &output) { - std::string buf; - unsigned int i = 0; - if (!push(s)) { std::cout << GridLogWarning << "XML: cannot open node '" << s << "'"; std::cout << std::endl; - - return; + } else { + for(unsigned int i = 0; node_.child("elem"); ) + { + output.resize(i + 1); + read("elem", output[i++]); + node_.child("elem").set_name("elem-done"); + } + pop(); + } + } + + template + void XmlReader::readMultiDim(const std::string &s, std::vector &buf, std::vector &dim) + { + if (!push(s)) + { + std::cout << GridLogWarning << "XML: cannot open node '" << s << "'"; + std::cout << std::endl; + } else { + static const char sName[] = "tensor"; + static const char sNameDone[] = "tensor-done"; + int Rank; + read("rank", Rank); + dim.resize( Rank ); + size_t NumElements = 1; + for( auto &d : dim ) + { + read("dim", d); + node_.child("dim").set_name("dim-done"); + NumElements *= d; + } + buf.resize( NumElements ); + std::vector MyIndex( Rank ); + for( int i = 0 ; i < Rank ; i++ ) { + MyIndex[i] = 0; + push(sName); + } + + for( auto &x : buf ) + { + NumElements--; + read("elem", x); + node_.child("elem").set_name("elem-done"); + int i; + for( i = Rank - 1 ; i != -1 && ++MyIndex[i] == dim[i] ; i-- ) + MyIndex[i] = 0; + int Rollover = Rank - 1 - i; + for( i = 0 ; i < Rollover ; i++ ) { + node_.set_name(sNameDone); + pop(); + } + for( i = 0 ; NumElements && i < Rollover ; i++ ) + push(sName); + } + pop(); } - while (node_.child("elem")) - { - output.resize(i + 1); - read("elem", output[i]); - node_.child("elem").set_name("elem-done"); - i++; - } - pop(); } - } #endif diff --git a/Grid/simd/Grid_vector_types.h b/Grid/simd/Grid_vector_types.h index c67e74cb..707af211 100644 --- a/Grid/simd/Grid_vector_types.h +++ b/Grid/simd/Grid_vector_types.h @@ -10,6 +10,7 @@ Author: Azusa Yamaguchi Author: Guido Cossu Author: Peter Boyle Author: neo +Author: Michael Marshall This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -89,17 +90,25 @@ template using NotEnableIf = Invoke struct is_complex : public std::false_type {}; -template <> struct is_complex > : public std::true_type {}; -template <> struct is_complex > : public std::true_type {}; +template <> struct is_complex : public std::true_type {}; +template <> struct is_complex : public std::true_type {}; -template using IfReal = Invoke::value, int> >; +template struct is_real : public std::false_type {}; +template struct is_real::value, + void>::type> : public std::true_type {}; + +template struct is_integer : public std::false_type {}; +template struct is_integer::value, + void>::type> : public std::true_type {}; + +template using IfReal = Invoke::value, int> >; template using IfComplex = Invoke::value, int> >; -template using IfInteger = Invoke::value, int> >; +template using IfInteger = Invoke::value, int> >; template using IfSame = Invoke::value, int> >; -template using IfNotReal = Invoke::value, int> >; +template using IfNotReal = Invoke::value, int> >; template using IfNotComplex = Invoke::value, int> >; -template using IfNotInteger = Invoke::value, int> >; +template using IfNotInteger = Invoke::value, int> >; template using IfNotSame = Invoke::value, int> >; //////////////////////////////////////////////////////// @@ -857,8 +866,10 @@ template struct is_simd : public std::false_type {}; template <> struct is_simd : public std::true_type {}; template <> struct is_simd : public std::true_type {}; +template <> struct is_simd : public std::true_type {}; template <> struct is_simd : public std::true_type {}; template <> struct is_simd : public std::true_type {}; +template <> struct is_simd : public std::true_type {}; template <> struct is_simd : public std::true_type {}; template using IfSimd = Invoke::value, int> >; diff --git a/Grid/tensors/Tensor_class.h b/Grid/tensors/Tensor_class.h index c7f868db..d59640df 100644 --- a/Grid/tensors/Tensor_class.h +++ b/Grid/tensors/Tensor_class.h @@ -5,6 +5,7 @@ Copyright (C) 2015 Author: Azusa Yamaguchi Author: Peter Boyle +Author: Michael Marshall This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -42,27 +43,26 @@ namespace Grid { // class GridTensorBase {}; +// Too late to remove these traits from Grid Tensors, so inherit from GridTypeMapper +#define GridVector_CopyTraits \ + using element = vtype; \ + using scalar_type = typename Traits::scalar_type; \ + using vector_type = typename Traits::vector_type; \ + using vector_typeD = typename Traits::vector_typeD; \ + using tensor_reduced = typename Traits::tensor_reduced; \ + using scalar_object = typename Traits::scalar_object; \ + using Complexified = typename Traits::Complexified; \ + using Realified = typename Traits::Realified; \ + using DoublePrecision = typename Traits::DoublePrecision; \ + static constexpr int TensorLevel = Traits::TensorLevel + template class iScalar { public: vtype _internal; - typedef vtype element; - typedef typename GridTypeMapper::scalar_type scalar_type; - typedef typename GridTypeMapper::vector_type vector_type; - typedef typename GridTypeMapper::vector_typeD vector_typeD; - typedef typename GridTypeMapper::tensor_reduced tensor_reduced_v; - typedef typename GridTypeMapper::scalar_object recurse_scalar_object; - typedef iScalar tensor_reduced; - typedef iScalar scalar_object; - // substitutes a real or complex version with same tensor structure - typedef iScalar::Complexified> Complexified; - typedef iScalar::Realified> Realified; - - // get double precision version - typedef iScalar::DoublePrecision> DoublePrecision; - - enum { TensorLevel = GridTypeMapper::TensorLevel + 1 }; + using Traits = GridTypeMapper >; + GridVector_CopyTraits; // Scalar no action // template using tensor_reduce_level = typename @@ -173,7 +173,10 @@ class iScalar { return stream; }; - + strong_inline const scalar_type * begin() const { return reinterpret_cast(&_internal); } + strong_inline scalar_type * begin() { return reinterpret_cast< scalar_type *>(&_internal); } + strong_inline const scalar_type * end() const { return begin() + Traits::count; } + strong_inline scalar_type * end() { return begin() + Traits::count; } }; /////////////////////////////////////////////////////////// // Allows to turn scalar>>> back to double. @@ -194,22 +197,9 @@ class iVector { public: vtype _internal[N]; - typedef vtype element; - typedef typename GridTypeMapper::scalar_type scalar_type; - typedef typename GridTypeMapper::vector_type vector_type; - typedef typename GridTypeMapper::vector_typeD vector_typeD; - typedef typename GridTypeMapper::tensor_reduced tensor_reduced_v; - typedef typename GridTypeMapper::scalar_object recurse_scalar_object; - typedef iScalar tensor_reduced; - typedef iVector scalar_object; + using Traits = GridTypeMapper >; + GridVector_CopyTraits; - // substitutes a real or complex version with same tensor structure - typedef iVector::Complexified, N> Complexified; - typedef iVector::Realified, N> Realified; - - // get double precision version - typedef iVector::DoublePrecision, N> DoublePrecision; - template ::value, T>::type * = nullptr> strong_inline auto operator=(T arg) -> iVector { @@ -218,7 +208,6 @@ class iVector { return *this; } - enum { TensorLevel = GridTypeMapper::TensorLevel + 1 }; iVector(const Zero &z) { *this = zero; }; iVector() = default; /* @@ -303,6 +292,11 @@ class iVector { // strong_inline vtype && operator ()(int i) { // return _internal[i]; // } + + strong_inline const scalar_type * begin() const { return reinterpret_cast(_internal); } + strong_inline scalar_type * begin() { return reinterpret_cast< scalar_type *>(_internal); } + strong_inline const scalar_type * end() const { return begin() + Traits::count; } + strong_inline scalar_type * end() { return begin() + Traits::count; } }; template @@ -310,25 +304,8 @@ class iMatrix { public: vtype _internal[N][N]; - typedef vtype element; - typedef typename GridTypeMapper::scalar_type scalar_type; - typedef typename GridTypeMapper::vector_type vector_type; - typedef typename GridTypeMapper::vector_typeD vector_typeD; - typedef typename GridTypeMapper::tensor_reduced tensor_reduced_v; - typedef typename GridTypeMapper::scalar_object recurse_scalar_object; - - // substitutes a real or complex version with same tensor structure - typedef iMatrix::Complexified, N> Complexified; - typedef iMatrix::Realified, N> Realified; - - // get double precision version - typedef iMatrix::DoublePrecision, N> DoublePrecision; - - // Tensor removal - typedef iScalar tensor_reduced; - typedef iMatrix scalar_object; - - enum { TensorLevel = GridTypeMapper::TensorLevel + 1 }; + using Traits = GridTypeMapper >; + GridVector_CopyTraits; iMatrix(const Zero &z) { *this = zero; }; iMatrix() = default; @@ -458,6 +435,11 @@ class iMatrix { // strong_inline vtype && operator ()(int i,int j) { // return _internal[i][j]; // } + + strong_inline const scalar_type * begin() const { return reinterpret_cast(_internal[0]); } + strong_inline scalar_type * begin() { return reinterpret_cast< scalar_type *>(_internal[0]); } + strong_inline const scalar_type * end() const { return begin() + Traits::count; } + strong_inline scalar_type * end() { return begin() + Traits::count; } }; template @@ -480,6 +462,3 @@ void vprefetch(const iMatrix &vv) { } } #endif - - - diff --git a/Grid/tensors/Tensor_traits.h b/Grid/tensors/Tensor_traits.h index c1ef397a..45a1d807 100644 --- a/Grid/tensors/Tensor_traits.h +++ b/Grid/tensors/Tensor_traits.h @@ -5,6 +5,7 @@ Author: Azusa Yamaguchi Author: Peter Boyle Author: Christopher Kelly +Author: Michael Marshall This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or @@ -26,6 +27,17 @@ Author: Christopher Kelly namespace Grid { + // Forward declarations + template class iScalar; + template class iVector; + template class iMatrix; + + // These are the Grid tensors + template struct isGridTensor : public std::false_type { static constexpr bool notvalue = true; }; + template struct isGridTensor> : public std::true_type { static constexpr bool notvalue = false; }; + template struct isGridTensor> : public std::true_type { static constexpr bool notvalue = false; }; + template struct isGridTensor> : public std::true_type { static constexpr bool notvalue = false; }; + ////////////////////////////////////////////////////////////////////////////////// // Want to recurse: GridTypeMapper >::scalar_type == ComplexD. // Use of a helper class like this allows us to template specialise and "dress" @@ -40,25 +52,26 @@ namespace Grid { // to study C++11's type_traits.h file. (std::enable_if >) // ////////////////////////////////////////////////////////////////////////////////// - - template class GridTypeMapper { - public: - typedef typename T::scalar_type scalar_type; - typedef typename T::vector_type vector_type; - typedef typename T::vector_typeD vector_typeD; - typedef typename T::tensor_reduced tensor_reduced; - typedef typename T::scalar_object scalar_object; - typedef typename T::Complexified Complexified; - typedef typename T::Realified Realified; - typedef typename T::DoublePrecision DoublePrecision; - enum { TensorLevel = T::TensorLevel }; + + // This saves repeating common properties for supported Grid Scalar types + // TensorLevel How many nested grid tensors + // Rank Rank of the grid tensor + // count Total number of elements, i.e. product of dimensions + // Dimension(dim) Size of dimension dim + struct GridTypeMapper_Base { + static constexpr int TensorLevel = 0; + static constexpr int Rank = 0; + static constexpr std::size_t count = 1; + static constexpr int Dimension(int dim) { return 0; } }; ////////////////////////////////////////////////////////////////////////////////// // Recursion stops with these template specialisations ////////////////////////////////////////////////////////////////////////////////// - template<> class GridTypeMapper { - public: + + template struct GridTypeMapper {}; + + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealF scalar_type; typedef RealF vector_type; typedef RealD vector_typeD; @@ -67,10 +80,8 @@ namespace Grid { typedef ComplexF Complexified; typedef RealF Realified; typedef RealD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealD scalar_type; typedef RealD vector_type; typedef RealD vector_typeD; @@ -79,10 +90,8 @@ namespace Grid { typedef ComplexD Complexified; typedef RealD Realified; typedef RealD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexF scalar_type; typedef ComplexF vector_type; typedef ComplexD vector_typeD; @@ -91,10 +100,8 @@ namespace Grid { typedef ComplexF Complexified; typedef RealF Realified; typedef ComplexD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexD scalar_type; typedef ComplexD vector_type; typedef ComplexD vector_typeD; @@ -103,10 +110,8 @@ namespace Grid { typedef ComplexD Complexified; typedef RealD Realified; typedef ComplexD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef Integer scalar_type; typedef Integer vector_type; typedef Integer vector_typeD; @@ -115,11 +120,9 @@ namespace Grid { typedef void Complexified; typedef void Realified; typedef void DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealF scalar_type; typedef vRealF vector_type; typedef vRealD vector_typeD; @@ -128,10 +131,8 @@ namespace Grid { typedef vComplexF Complexified; typedef vRealF Realified; typedef vRealD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealD scalar_type; typedef vRealD vector_type; typedef vRealD vector_typeD; @@ -140,10 +141,20 @@ namespace Grid { typedef vComplexD Complexified; typedef vRealD Realified; typedef vRealD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { + // Fixme this is incomplete until Grid supports fp16 or bfp16 arithmetic types + typedef RealF scalar_type; + typedef vRealH vector_type; + typedef vRealD vector_typeD; + typedef vRealH tensor_reduced; + typedef RealF scalar_object; + typedef vComplexH Complexified; + typedef vRealH Realified; + typedef vRealD DoublePrecision; + }; + template<> struct GridTypeMapper : public GridTypeMapper_Base { + // Fixme this is incomplete until Grid supports fp16 or bfp16 arithmetic types typedef ComplexF scalar_type; typedef vComplexH vector_type; typedef vComplexD vector_typeD; @@ -152,10 +163,8 @@ namespace Grid { typedef vComplexH Complexified; typedef vRealH Realified; typedef vComplexD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexF scalar_type; typedef vComplexF vector_type; typedef vComplexD vector_typeD; @@ -164,10 +173,8 @@ namespace Grid { typedef vComplexF Complexified; typedef vRealF Realified; typedef vComplexD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexD scalar_type; typedef vComplexD vector_type; typedef vComplexD vector_typeD; @@ -176,10 +183,8 @@ namespace Grid { typedef vComplexD Complexified; typedef vRealD Realified; typedef vComplexD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef Integer scalar_type; typedef vInteger vector_type; typedef vInteger vector_typeD; @@ -188,57 +193,52 @@ namespace Grid { typedef void Complexified; typedef void Realified; typedef void DoublePrecision; - enum { TensorLevel = 0 }; }; - // First some of my own traits - template struct isGridTensor { - static const bool value = true; - static const bool notvalue = false; +#define GridTypeMapper_RepeatedTypes \ + using BaseTraits = GridTypeMapper; \ + using scalar_type = typename BaseTraits::scalar_type; \ + using vector_type = typename BaseTraits::vector_type; \ + using vector_typeD = typename BaseTraits::vector_typeD; \ + static constexpr int TensorLevel = BaseTraits::TensorLevel + 1 + + template struct GridTypeMapper> { + GridTypeMapper_RepeatedTypes; + using tensor_reduced = iScalar; + using scalar_object = iScalar; + using Complexified = iScalar; + using Realified = iScalar; + using DoublePrecision = iScalar; + static constexpr int Rank = BaseTraits::Rank + 1; + static constexpr std::size_t count = BaseTraits::count; + static constexpr int Dimension(int dim) { + return ( dim == 0 ) ? 1 : BaseTraits::Dimension(dim - 1); } }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; + + template struct GridTypeMapper> { + GridTypeMapper_RepeatedTypes; + using tensor_reduced = iScalar; + using scalar_object = iVector; + using Complexified = iVector; + using Realified = iVector; + using DoublePrecision = iVector; + static constexpr int Rank = BaseTraits::Rank + 1; + static constexpr std::size_t count = BaseTraits::count * N; + static constexpr int Dimension(int dim) { + return ( dim == 0 ) ? N : BaseTraits::Dimension(dim - 1); } }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; + + template struct GridTypeMapper> { + GridTypeMapper_RepeatedTypes; + using tensor_reduced = iScalar; + using scalar_object = iMatrix; + using Complexified = iMatrix; + using Realified = iMatrix; + using DoublePrecision = iMatrix; + static constexpr int Rank = BaseTraits::Rank + 2; + static constexpr std::size_t count = BaseTraits::count * N * N; + static constexpr int Dimension(int dim) { + return ( dim == 0 || dim == 1 ) ? N : BaseTraits::Dimension(dim - 2); } }; // Match the index @@ -263,20 +263,13 @@ namespace Grid { typedef T type; }; - //Query if a tensor or Lattice is SIMD vector or scalar - template - class isSIMDvectorized{ - template - static typename std::enable_if< !std::is_same< typename GridTypeMapper::type>::scalar_type, - typename GridTypeMapper::type>::vector_type>::value, char>::type test(void *); + //Query whether a tensor or Lattice is SIMD vector or scalar + template struct isSIMDvectorized : public std::false_type {}; + template struct isSIMDvectorized::type>::scalar_type, + typename GridTypeMapper::type>::vector_type>::value, void>::type> + : public std::true_type {}; - template - static double test(...); - - public: - enum {value = sizeof(test(0)) == sizeof(char) }; - }; - //Get the precision of a Lattice, tensor or scalar type in units of sizeof(float) template class getPrecision{ diff --git a/Hadrons/Application.cc b/Hadrons/Application.cc index be2e87a1..e6fb85fa 100644 --- a/Hadrons/Application.cc +++ b/Hadrons/Application.cc @@ -48,28 +48,32 @@ Application::Application(void) { initLogger(); auto dim = GridDefaultLatt(), mpi = GridDefaultMpi(), loc(dim); - locVol_ = 1; - for (unsigned int d = 0; d < dim.size(); ++d) + + if (dim.size()) { - loc[d] /= mpi[d]; - locVol_ *= loc[d]; + locVol_ = 1; + for (unsigned int d = 0; d < dim.size(); ++d) + { + loc[d] /= mpi[d]; + locVol_ *= loc[d]; + } + LOG(Message) << "====== HADRONS APPLICATION INITIALISATION ======" << std::endl; + LOG(Message) << "** Dimensions" << std::endl; + LOG(Message) << "Global lattice: " << dim << std::endl; + LOG(Message) << "MPI partition : " << mpi << std::endl; + LOG(Message) << "Local lattice : " << loc << std::endl; + LOG(Message) << std::endl; + LOG(Message) << "** Default parameters (and associated C macros)" << std::endl; + LOG(Message) << "ASCII output precision : " << MACOUT(DEFAULT_ASCII_PREC) << std::endl; + LOG(Message) << "Fermion implementation : " << MACOUTS(FIMPLBASE) << std::endl; + LOG(Message) << "z-Fermion implementation: " << MACOUTS(ZFIMPLBASE) << std::endl; + LOG(Message) << "Scalar implementation : " << MACOUTS(SIMPLBASE) << std::endl; + LOG(Message) << "Gauge implementation : " << MACOUTS(GIMPLBASE) << std::endl; + LOG(Message) << "Eigenvector base size : " + << MACOUT(HADRONS_DEFAULT_LANCZOS_NBASIS) << std::endl; + LOG(Message) << "Schur decomposition : " << MACOUTS(HADRONS_DEFAULT_SCHUR) << std::endl; + LOG(Message) << std::endl; } - LOG(Message) << "====== HADRONS APPLICATION INITIALISATION ======" << std::endl; - LOG(Message) << "** Dimensions" << std::endl; - LOG(Message) << "Global lattice: " << dim << std::endl; - LOG(Message) << "MPI partition : " << mpi << std::endl; - LOG(Message) << "Local lattice : " << loc << std::endl; - LOG(Message) << std::endl; - LOG(Message) << "** Default parameters (and associated C macros)" << std::endl; - LOG(Message) << "ASCII output precision : " << MACOUT(DEFAULT_ASCII_PREC) << std::endl; - LOG(Message) << "Fermion implementation : " << MACOUTS(FIMPLBASE) << std::endl; - LOG(Message) << "z-Fermion implementation: " << MACOUTS(ZFIMPLBASE) << std::endl; - LOG(Message) << "Scalar implementation : " << MACOUTS(SIMPLBASE) << std::endl; - LOG(Message) << "Gauge implementation : " << MACOUTS(GIMPLBASE) << std::endl; - LOG(Message) << "Eigenvector base size : " - << MACOUT(HADRONS_DEFAULT_LANCZOS_NBASIS) << std::endl; - LOG(Message) << "Schur decomposition : " << MACOUTS(HADRONS_DEFAULT_SCHUR) << std::endl; - LOG(Message) << std::endl; } Application::Application(const Application::GlobalPar &par) diff --git a/Hadrons/Environment.cc b/Hadrons/Environment.cc index e7a7ac55..540b3620 100644 --- a/Hadrons/Environment.cc +++ b/Hadrons/Environment.cc @@ -43,15 +43,13 @@ HADRONS_ERROR_REF(ObjectDefinition, "no object with address " + std::to_string(a // constructor ///////////////////////////////////////////////////////////////// Environment::Environment(void) { - dim_ = GridDefaultLatt(); - nd_ = dim_.size(); - createGrid(1); + dim_ = GridDefaultLatt(); + nd_ = dim_.size(); vol_ = 1.; for (auto d: dim_) { vol_ *= d; } - rng4d_.reset(new GridParallelRNG(getGrid())); } // grids /////////////////////////////////////////////////////////////////////// @@ -76,8 +74,13 @@ double Environment::getVolume(void) const } // random number generator ///////////////////////////////////////////////////// -GridParallelRNG * Environment::get4dRng(void) const +GridParallelRNG * Environment::get4dRng(void) { + if (rng4d_ == nullptr) + { + rng4d_.reset(new GridParallelRNG(getGrid())); + } + return rng4d_.get(); } diff --git a/Hadrons/Environment.hpp b/Hadrons/Environment.hpp index b77345c8..24f6b31e 100644 --- a/Hadrons/Environment.hpp +++ b/Hadrons/Environment.hpp @@ -113,7 +113,7 @@ public: unsigned int getNd(void) const; double getVolume(void) const; // random number generator - GridParallelRNG * get4dRng(void) const; + GridParallelRNG * get4dRng(void); // general memory management void addObject(const std::string name, const int moduleAddress = -1); @@ -182,7 +182,7 @@ private: std::map gridCoarse5d_; unsigned int nd_; // random number generator - RngPt rng4d_; + RngPt rng4d_{nullptr}; // object store std::vector object_; std::map objectAddress_; diff --git a/Hadrons/Modules/MAction/DWF.hpp b/Hadrons/Modules/MAction/DWF.hpp index bac70ae8..21024a4f 100644 --- a/Hadrons/Modules/MAction/DWF.hpp +++ b/Hadrons/Modules/MAction/DWF.hpp @@ -112,8 +112,6 @@ void TDWF::setup(void) << par().mass << ", M5= " << par().M5 << " and Ls= " << par().Ls << " using gauge field '" << par().gauge << "'" << std::endl; - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; auto &U = envGet(GaugeField, par().gauge); auto &g4 = *envGetGrid(FermionField); @@ -121,8 +119,26 @@ void TDWF::setup(void) auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); typename DomainWallFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, DomainWallFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, implParams); } diff --git a/Hadrons/Modules/MAction/MobiusDWF.hpp b/Hadrons/Modules/MAction/MobiusDWF.hpp index 6b960e2e..3ce02b91 100644 --- a/Hadrons/Modules/MAction/MobiusDWF.hpp +++ b/Hadrons/Modules/MAction/MobiusDWF.hpp @@ -112,17 +112,33 @@ void TMobiusDWF::setup(void) << ", b= " << par().b << ", c= " << par().c << " using gauge field '" << par().gauge << "'" << std::endl; - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; - + auto &U = envGet(GaugeField, par().gauge); auto &g4 = *envGetGrid(FermionField); auto &grb4 = *envGetRbGrid(FermionField); auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); typename MobiusFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, MobiusFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, par().b, par().c, implParams); diff --git a/Hadrons/Modules/MAction/ScaledDWF.hpp b/Hadrons/Modules/MAction/ScaledDWF.hpp index 8742a820..3e052225 100644 --- a/Hadrons/Modules/MAction/ScaledDWF.hpp +++ b/Hadrons/Modules/MAction/ScaledDWF.hpp @@ -111,8 +111,6 @@ void TScaledDWF::setup(void) << ", scale= " << par().scale << " using gauge field '" << par().gauge << "'" << std::endl; - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; auto &U = envGet(GaugeField, par().gauge); auto &g4 = *envGetGrid(FermionField); @@ -120,8 +118,26 @@ void TScaledDWF::setup(void) auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); typename ScaledShamirFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, ScaledShamirFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, par().scale, implParams); diff --git a/Hadrons/Modules/MAction/Wilson.hpp b/Hadrons/Modules/MAction/Wilson.hpp index 3ce5569e..11df89c3 100644 --- a/Hadrons/Modules/MAction/Wilson.hpp +++ b/Hadrons/Modules/MAction/Wilson.hpp @@ -109,15 +109,31 @@ void TWilson::setup(void) { LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass << " using gauge field '" << par().gauge << "'" << std::endl; - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; auto &U = envGet(GaugeField, par().gauge); auto &grid = *envGetGrid(FermionField); auto &gridRb = *envGetRbGrid(FermionField); typename WilsonFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, WilsonFermion, getName(), 1, U, grid, gridRb, par().mass, implParams); } diff --git a/Hadrons/Modules/MAction/WilsonClover.hpp b/Hadrons/Modules/MAction/WilsonClover.hpp index a9ea8de5..d80bccdf 100644 --- a/Hadrons/Modules/MAction/WilsonClover.hpp +++ b/Hadrons/Modules/MAction/WilsonClover.hpp @@ -112,17 +112,34 @@ void TWilsonClover::setup(void) { LOG(Message) << "Setting up Wilson clover fermion matrix with m= " << par().mass << " using gauge field '" << par().gauge << "'" << std::endl; - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; LOG(Message) << "Clover term csw_r: " << par().csw_r << " csw_t: " << par().csw_t << std::endl; + auto &U = envGet(GaugeField, par().gauge); auto &grid = *envGetGrid(FermionField); auto &gridRb = *envGetRbGrid(FermionField); typename WilsonCloverFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, WilsonCloverFermion, getName(), 1, U, grid, gridRb, par().mass, par().csw_r, par().csw_t, par().clover_anisotropy, implParams); diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/Hadrons/Modules/MAction/ZMobiusDWF.hpp index e780bd73..56efe1a2 100644 --- a/Hadrons/Modules/MAction/ZMobiusDWF.hpp +++ b/Hadrons/Modules/MAction/ZMobiusDWF.hpp @@ -118,10 +118,7 @@ void TZMobiusDWF::setup(void) { LOG(Message) << " omega[" << i << "]= " << par().omega[i] << std::endl; } - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; - env().createGrid(par().Ls); auto &U = envGet(GaugeField, par().gauge); auto &g4 = *envGetGrid(FermionField); auto &grb4 = *envGetRbGrid(FermionField); @@ -129,8 +126,26 @@ void TZMobiusDWF::setup(void) auto &grb5 = *envGetRbGrid(FermionField, par().Ls); auto omega = par().omega; typename ZMobiusFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, ZMobiusFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, omega, par().b, par().c, implParams); diff --git a/Hadrons/Modules/MGauge/GaugeFix.hpp b/Hadrons/Modules/MGauge/GaugeFix.hpp index 87822a7e..57b1e082 100644 --- a/Hadrons/Modules/MGauge/GaugeFix.hpp +++ b/Hadrons/Modules/MGauge/GaugeFix.hpp @@ -7,6 +7,7 @@ Source file: Hadrons/Modules/MGauge/GaugeFix.hpp Copyright (C) 2015-2019 Author: Antonin Portelli +Author: Nils Asmussen Author: Peter Boyle This program is free software; you can redistribute it and/or modify @@ -58,6 +59,7 @@ class TGaugeFix: public Module { public: GAUGE_TYPE_ALIASES(GImpl,); + typedef typename GImpl::GaugeLinkField GaugeMat; public: // constructor TGaugeFix(const std::string name); @@ -94,7 +96,7 @@ std::vector TGaugeFix::getInput(void) template std::vector TGaugeFix::getOutput(void) { - std::vector out = {getName()}; + std::vector out = {getName(), getName()+"_xform"}; return out; } @@ -103,6 +105,7 @@ template void TGaugeFix::setup(void) { envCreateLat(GaugeField, getName()); + envCreateLat(GaugeMat, getName()+"_xform"); } @@ -116,6 +119,7 @@ void TGaugeFix::execute(void) LOG(Message) << par().gauge << std::endl; auto &U = envGet(GaugeField, par().gauge); auto &Umu = envGet(GaugeField, getName()); + auto &xform = envGet(GaugeMat, getName()+"_xform"); LOG(Message) << "Gauge Field fetched" << std::endl; //do we allow maxiter etc to be user set? Real alpha = par().alpha; @@ -123,8 +127,8 @@ void TGaugeFix::execute(void) Real Omega_tol = par().Omega_tol; Real Phi_tol = par().Phi_tol; bool Fourier = par().Fourier; - FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(U,alpha,maxiter,Omega_tol,Phi_tol,Fourier); Umu = U; + FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier); LOG(Message) << "Gauge Fixed" << std::endl; } diff --git a/Hadrons/Utilities/HadronsXmlValidate.cc b/Hadrons/Utilities/HadronsXmlValidate.cc new file mode 100644 index 00000000..73cf3139 --- /dev/null +++ b/Hadrons/Utilities/HadronsXmlValidate.cc @@ -0,0 +1,64 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Utilities/HadronsXmlValidate.cc + +Copyright (C) 2015-2019 + +Author: Antonin Portelli + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace QCD; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // parse command line + std::string parameterFileName; + + if (argc != 2) + { + std::cerr << "usage: " << argv[0] << " "; + std::cerr << std::endl; + std::exit(EXIT_FAILURE); + } + parameterFileName = argv[1]; + + try + { + Application application(parameterFileName); + + application.parseParameterFile(parameterFileName); + auto &vm = VirtualMachine::getInstance(); + vm.getModuleGraph(); + LOG(Message) << "Application valid (check XML warnings though)" + << std::endl; + } + catch (const std::exception& e) + { + Exceptions::abort(e); + } + + return EXIT_SUCCESS; +} diff --git a/Hadrons/Utilities/Makefile.am b/Hadrons/Utilities/Makefile.am index 4f324d6d..6c48d966 100644 --- a/Hadrons/Utilities/Makefile.am +++ b/Hadrons/Utilities/Makefile.am @@ -1,8 +1,11 @@ -bin_PROGRAMS = HadronsXmlRun HadronsFermionEP64To32 HadronsContractor HadronsContractorBenchmark +bin_PROGRAMS = HadronsXmlRun HadronsXmlValidate HadronsFermionEP64To32 HadronsContractor HadronsContractorBenchmark HadronsXmlRun_SOURCES = HadronsXmlRun.cc HadronsXmlRun_LDADD = ../libHadrons.a ../../Grid/libGrid.a +HadronsXmlValidate_SOURCES = HadronsXmlValidate.cc +HadronsXmlValidate_LDADD = ../libHadrons.a ../../Grid/libGrid.a + HadronsFermionEP64To32_SOURCES = EigenPackCast.cc HadronsFermionEP64To32_CXXFLAGS = $(AM_CXXFLAGS) -DFIN=WilsonImplD::FermionField -DFOUT=WilsonImplF::FermionField HadronsFermionEP64To32_LDADD = ../libHadrons.a ../../Grid/libGrid.a diff --git a/documentation/GridXcode/GridXcFig1.png b/documentation/GridXcode/GridXcFig1.png new file mode 100644 index 00000000..73fd0625 Binary files /dev/null and b/documentation/GridXcode/GridXcFig1.png differ diff --git a/documentation/GridXcode/GridXcFig2.png b/documentation/GridXcode/GridXcFig2.png new file mode 100644 index 00000000..9b169491 Binary files /dev/null and b/documentation/GridXcode/GridXcFig2.png differ diff --git a/documentation/GridXcode/GridXcFig3.png b/documentation/GridXcode/GridXcFig3.png new file mode 100644 index 00000000..02595468 Binary files /dev/null and b/documentation/GridXcode/GridXcFig3.png differ diff --git a/documentation/GridXcode/GridXcFig4.png b/documentation/GridXcode/GridXcFig4.png new file mode 100644 index 00000000..5abd150b Binary files /dev/null and b/documentation/GridXcode/GridXcFig4.png differ diff --git a/documentation/GridXcode/readme.md b/documentation/GridXcode/readme.md new file mode 100644 index 00000000..b75d9323 --- /dev/null +++ b/documentation/GridXcode/readme.md @@ -0,0 +1,422 @@ +# Using Xcode with Grid on Mac OS + +This guide explains how to use Xcode as an IDE on Mac OS. + +# Initial setup + +For first time setup of the Xcode and Grid build environment on Mac OS, you will need to do the following, in this order: + +1. Install Xcode and the Xcode command-line utilities +2. Set Grid environment variables +3. Install and build Open MPI ***optional*** +4. Install and build Grid pre-requisites +5. Install, Configure and Build Grid + +Apple's [Xcode website][Xcode] is the go-to reference for 1, and the definitive reference for 4 and 5 is the [Grid Documentation][GridDoc]. + +[Xcode]: https://developer.apple.com/xcode/ +[GridDoc]: https://github.com/paboyle/Grid/blob/develop/documentation/Grid.pdf + +The following sections explain these steps in more detail + +## 1. Install Xcode and the Xcode command-line utilities + +See Apple's [Xcode website][Xcode] for instructions on installing Xcode. + +Once Xcode is installed, install the Xcode command-line utilities using: + + xcode-select --install + +*NB: the screenshots from this guide were generated from Xcode 10.1.* + +## 2. Set Grid environment variables + +To make sure we can share Xcode projects via git and have them work without requiring modification, we will define Grid environment variables. To make sure these environment variables will be available to the Xcode build system, issue the following command: + + defaults write com.apple.dt.Xcode UseSanitizedBuildSystemEnvironment -bool NO + +These are the environment variables we will define for Grid: + +Variable | Typical Value | Use +--- | --- | --- +`Grid` | `/Users/user_id/src/Grid` | Path to grid source +`GridPre` | `/Users/user_id/bin` | Path to install directory containing grid pre-requisites built from source +`GridPkg` | **MacPorts**=`/opt/local`, **Homebrew**=`/usr/local` | Path to package manager install directory + +Choose either of the following ways to do this, and when you're done, log out and in again. To check these have been set: + + printenv|grep -i grid + +### Method 1 -- Apple Script + +* Start *Script Editor* (cmd-space, *script editor*) +* Click on *New Document*. Paste the following into the new script, editing the paths appropriately (just replace `user_id` with your *user_id* if you are unsure): + +```apple script +do shell script "launchctl setenv Grid $HOME/src/Grid +launchctl setenv GridPre $HOME/bin +launchctl setenv GridPkg /opt/local" +``` + +* Save the script inside `~/Applications` and give it the name `GridEnv.app`. +* Open `System Preferences`, `Users & Groups` +* Click on `Login Items` +* Click the plus sign to add a new login item +* Select the `~/Applications` folder and select `GridEnv.app` + +Log out and in again. + +### Method 2 -- `environment.plist` + +Make the file `environment.plist` in `~/Library/LaunchAgents` with the following contents, editing the paths appropriately (just replace `user_id` with your *user_id* if you are unsure): + +```html + + + + +Label +Grid.startup +ProgramArguments + +sh +-c +launchctl setenv Grid $HOME/src/Grid +launchctl setenv GridPre $HOME/bin +launchctl setenv GridPkg /opt/local + + +RunAtLoad + + + +``` + +## 3. Install and build Open MPI -- ***optional*** + +Download the latest version of [Open MPI][OMPI] version 3.1 (I used 3.1.3). + +NB: Grid does not have any dependencies on fortran, however many standard scientific packages do, so you may as well download GNU fortran (e.g. MacPorts ``gfortran`` package) and build Open MPI like so: + +[OMPI]: https://www.open-mpi.org/software/ompi/v3.1/ + + ../configure CC=clang CXX=clang++ F77=gfortran FC=gfortran CXXFLAGS=-g --prefix=$GridPre/openmpi-3.1.3 + make -j 4 all install + +(If you don't want to bother with fortran bindings, just don't include the F77 and FC flags) + +## 4. Install and build Grid pre-requisites + +To simplify the installation of **Grid pre-requisites**, you can use your favourite package manager, e.g.: + +### 1. [MacPorts][MacPorts] + +[MacPorts]: https://www.macports.org "MacPorts package manager" + +Install [MacPorts][MacPorts] if you haven't done so already, and then install packages with: + + sudo port install + +These are the `portname`s for mandatory Grid libraries: + +* git +* gmp +* mpfr + +and these are the `portname`s for optional Grid libraries: + +* fftw-3 +* hdf5 +* lapack +* doxygen +* OpenBLAS + +***Please update this list with any packages I've missed! ... and double-check whether OpenBLAS is really for Grid*** + +### 2. [Homebrew][Homebrew] + +[Homebrew]: https://brew.sh "Homebrew package manager" + +Install [Homebrew][Homebrew] if you haven't done so already, and then install packages with: + + sudo brew install + +The same packages are available as from MacPorts. + +### Install LIME ***optional*** + +There isn't currently a port for [C-LIME][C-LIME], so download the source and then build it: + +[C-LIME]: https://usqcd-software.github.io/c-lime/ "C-language API for Lattice QCD Interchange Message Encapsulation / Large Internet Message Encapsulation" + + ../configure --prefix=$GridPre/lime-1.3.2 CC=clang + make -j 4 + make install + +## 5. Install, Configure and Build Grid + +### 5.1 Install Grid + +[Grid]: https://github.com/paboyle/Grid + +Start by cloning [Grid (from GitHub)][Grid] ***into the directory you specified in*** `$Grid`. Bear in mind that git will create the `Grid` subdirectory to place Grid in, so for example if `$Grid` is set to `~/src/Grid` then install Grid with: + + cd ~/src + +followed by either: + + git clone git@github.com:paboyle/Grid.git + +or + + git clone https://github.com/paboyle/Grid.git + +depending on whether you are using https or ssh. + +### 5.2 Configure Grid + +The Xcode build system supports multiple configurations for each project, by default: `Debug` and `Release`, but more configurations can be defined. We will create separate Grid build directories for each configuration, using the Grid **Autoconf** build system to make each configuration. NB: it is **not** necessary to run `make install` on them once they are built (IDE features such as *jump to definition* will work better of you don't). + +Below are shown the `configure` script invocations for three recommended configurations. You are free to define more, less or different configurations, but as a minimum, be sure to build a `Debug` configuration. + +#### 1. `Debug` + +This is the build for every day developing and debugging with Xcode. It uses the Xcode clang c++ compiler, without MPI, and defaults to double-precision. Xcode builds the `Debug` configuration with debug symbols for full debugging: + + ../configure --with-hdf5=$GridPkg --with-gmp=$GridPkg --with-mpfr=$GridPkg --with-fftw=$GridPkg --with-lime=$GridPre/lime-1.3.2 --enable-simd=GEN --enable-precision=double CXX=clang++ --prefix=$GridPre/GridDebug --enable-comms=none --enable-doxygen-doc + +#### 2. `Release` + +Since Grid itself doesn't really have debug configurations, the release build is recommended to be the same as `Debug`, except using single-precision (handy for validation): + + ../configure --with-hdf5=$GridPkg --with-gmp=$GridPkg --with-mpfr=$GridPkg --with-fftw=$GridPkg --with-lime=$GridPre/lime-1.3.2 --enable-simd=GEN --enable-precision=single CXX=clang++ --prefix=$GridPre/GridRelease --enable-comms=none --enable-doxygen-doc + +#### 3. `MPIDebug` + +Debug configuration with MPI: + + ../configure --with-hdf5=$GridPkg --with-gmp=$GridPkg --with-mpfr=$GridPkg --with-fftw=$GridPkg --with-lime=$GridPre/lime-1.3.2 --enable-simd=GEN --enable-precision=double CXX=clang++ --prefix=$GridPre/GridMPIDebug --enable-comms=mpi-auto MPICXX=$GridPre/openmpi-3.1.3/bin/mpicxx --enable-doxygen-doc + +### 5.3 Build Grid + +Each configuration must be built before they can be used. You can either: + +1. Use automake and the Grid Makefile with `make -j 4` (NB: you **do not** need to run `make install` for these to work with Xcode) +2. Build `Grid` and `Hadrons` under Xcode (see below) + +# Make a new application which links to Grid / Hadrons + +Making an Xcode project which links to Grid / Hadrons is straightforward: + +* Make a new application (in the usual way) +* Configure your application to use Grid (via three project settings:) + 1. `HEADER_SEARCH_PATHS` + 2. `LIBRARY_SEARCH_PATHS` + 3. `OTHER_LDFLAGS` +* Make additional configurations, e.g. `MPIDebug` (NB Xcode will make `Debug` and `Release` by default) + +Detailed instructions follow, but instead of following the instructions in this section, you can clone `HelloGrid` from the [University of Edinburgh GitLab site][HelloGrid]. + +[HelloGrid]: https://git.ecdf.ed.ac.uk/s1786208/HelloGrid + +## Make a new application + +To make a hello world application for Grid: + +* Start Xcode +* Click 'Create a new project' +* Click ‘macOS’, then in the ‘Application’ section choose ‘Command Line Tool’, then click ‘Next’ +* Choose options for your new project: + * Product Name: HelloGrid + * Team: None + * Organisation Name: sopa + * Organisation Identifier: uk.ac.ed.ph + * Language: C++ + * ... then click ‘Next’ +* Choose a location for your project, e.g. `$HOME/src`. NB: The project and all it’s files will be created inside `$HOME/src/HelloGrid`. If you are using Git, you can put the new project under Git source control immediately, if you like. Now click ‘Create’. + +## Configure your new application to use Grid + +Click the project name (`HelloGrid`) in the project navigator pane on the left (command-1 if it's not visible), then click the project name (`HelloGrid`) under `PROJECT` in the second pane. Click the `Build Settings` tab on the right, then under that click `All` and `Combined`. You should see: + +![Project settings](GridXcFig1.png) + +We now need to make changes to two sections (these are listed in alphabetical order), bearing in mind that if you are not using MPI (or you gave your build directories different names) replace `build_mpidebug` and `build_mpirelease` with the directory names you used. + +### 1. Search Paths + +#### HEADER_SEARCH_PATHS + +Obtain a list of header locations required by Grid by running the following from your Grid build directory (choose an MPI configuration if you built one, e.g. `MPIDebug`): + + ./grid-config --cxxflags + +Output should look similar to: + + -I$GridPre/openmpi-3.1.3/include -I$GridPkg/include -I$GridPre/lime-1.3.2/include -I$GridPkg/include -I$GridPkg/include -I$GridPkg/include -O3 -g -std=c++11 + +The header locations follow the `-I` switches. You can ignore the other switches, and you can ignore duplicate entries, which just mean that your package manager has installed multiple packages in the same location. + +*Note: `grid-config` will output absolute paths. Make sure to replace absolute paths with environment variables (such as `$GridPre`) in your settings, so that the project will work unmodified for other collaborators downloading the same project from git.* + +Set HEADER_SEARCH_PATHS to: + + $Grid/build$(CONFIGURATION)/Grid + $Grid + $Grid/Grid + +followed by (***the order is important***) the locations reported by `grid-config --cxxflags`, ignoring duplicates, e.g.: + + $GridPre/openmpi-3.1.3/include + $GridPkg/include + $GridPre/lime-1.3.2/include + +**Note: the easiest way to set this value is to put it all on one line, space separated, and edit the text to the right of `HEADER_SEARCH_PATHS`**, i.e.: + + $Grid/build$(CONFIGURATION)/Grid $Grid $Grid/Grid $GridPre/openmpi-3.1.3/include $GridPkg/include $GridPre/lime-1.3.2/include + +#### LIBRARY_SEARCH_PATHS + +Obtain a list of library locations required by Grid by running the following from your Grid build directory (again, choose an MPI configuration if you built one, e.g. `MPIDebug`): + + ./grid-config --ldflags + +Output should look similar to: + + -L$GridPre/openmpi-3.1.3/lib -L$GridPkg/lib -L$GridPre/lime-1.3.2/lib -L$GridPkg/lib -L$GridPkg/lib -L$GridPkg/lib + +Paste the output ***with `$Grid/build$(CONFIGURATION)/Grid $Grid/build$(CONFIGURATION)/Hadrons ` prepended*** into `LIBRARY_SEARCH_PATHS`: + + $Grid/build$(CONFIGURATION)/Grid $Grid/build$(CONFIGURATION)/Hadrons $GridPre/openmpi-3.1.3/lib $GridPkg/lib $GridPre/lime-1.3.2/lib + +### 2. Linking + +#### OTHER_LDFLAGS + +The easiest way to link to all required libraries is to obtain a list of all libraries required by Grid by running the following from your Grid build directory: + + ./grid-config --libs + +and pasting the output ***with `-lGrid -lHadrons ` prepended*** (including the `-l` switches) directly into `OTHER_LDFLAGS`, e.g.: + + -lGrid -lHadrons -lmpi -lhdf5_cpp -lz -lcrypto -llime -lfftw3f -lfftw3 -lmpfr -lgmp -lstdc++ -lm -lz -lhdf5 + +## Make additional configurations + +On the project settings, `Info` tab, click the plus sign underneath configurations: + +![Add configurations](GridXcFig4.png) + +Choose `Duplicate "Debug" Configuration` (you can choose `Release` if you prefer) and give the new configuration a name, e.g. `MPIDebug`. + +## Edit your source code + +A hello world for grid is: + +```c++ +#include +using namespace Grid; + +int main(int argc, char * argv[]) { + Grid_init(&argc,&argv); + std::cout << GridLogMessage << "Hello Grid" << std::endl; + Grid_finalize(); + return 0; +} +``` + +## Create a `.gitignore` file for Xcode + +You can create an up-to-date .gitignore file to ignore all the Xcode temporary build files using [gitignore.io][GIO]. + +[GIO]: https://www.gitignore.io/api/xcode + +NB: If you let Xcode add your project to git when you created it, you probably want to remove your personal scheme selection from git: + + git rm --cached HelloGrid.xcodeproj/xcuserdata/$USER.xcuserdatad/xcschemes/xcschememanagement.plist + +## Run your program under the Xcode debugger + +First, specify command-line arguments. From the menu, select `Product`, then `Scheme`, then `Edit Scheme`. Select `Run` on the left, then select the `Arguments` tab on the right. Add the following to `Arguments passed on Launch`: + + --grid 4.4.4.8 + +If your program will be manipulating files, it's a good idea to specify the working directory on the `Options` tab under `Use Custom Working Directory` (by default, Xcode launches the program inside the Xcode build folder). + +Then click `Close`. + +Let's set a breakpoint by clicking on: + + Grid_finalize(); + +then from the menu selecting `Debug`, then `Breakpoints`, then `Add Breakpoint at Current Line`. + +Now click on the `Play` button (the right pointing triangle just to the right of the maximise button) to run your program under the debugger. (You may see dialog boxes the first couple of times asking whether to allow MPI to receive network requests - say yes to these.) + +The debug output pane opens at the bottom of Xcode, with output on the right (ending with `Hello Grid`) and local variables on the left i.e.: + +![Running under the debugger](GridXcFig2.png) + +See the Xcode documentation to learn about the debugger. When you're done, press `ctl-cmd-Y` to let the program run to completion. + +# Debugging multiple MPI processes under Xcode + +You could tell Xcode to use mpirun to launch multiple copies of a target executable, however if you do this the debugger will attach to mpirun - not your target process. + +Instead: + +1. Set a breakpoint just inside `main()` (otherwise your programs may complete before you attach to them all) + +2. From the `Debug` menu, select `Attach to Process by PID or Name ...`. In the `PID or Process Name` field, enter the name of your target. Then click `Attach`. + +3. From a terminal session, locate and run your executable using `mpirun` (*the mangled name of the project build products will not be exactly the same as this example*): + + `$GridPre/openmpi-3.1.3/bin/mpirun -np 2 ~/Library/Developer/Xcode/DerivedData/HelloGrid-fiyyuveptaqelbbvllomcgjyvghr/Build/Products/Debug/HelloGrid --grid 4.4.4.8 --mpi 1.1.1.2` + + The Xcode debugger will attach to the first process. + +4. From the `Debug` menu in Xcode, select `Attach to Process`, and other running instances of your application will appear at the top of the list. Attach to as many instances as you wish to debug. + +5. Click on the first process (which should have stopped at the breakpoint) and restart it with ctl-cmd-y + +You are now debugging multiple MPI instances, and the Xcode debugger should look similar to this: + +![Debugging multiple MPI instances under the Xcode debugger](GridXcFig3.png) + +# Build `Grid` and `Hadrons` libraries under Xcode + +If you want to build `Grid` and `Hadrons` libraries using Xcode, you will need to: + +1. Make new library targets for `Grid` and `Hadrons` + +2. Add Grid source folders to your project: + a. Right click project then `Add files to "project" ...` + b. Choose `$Grid/Grid` folder + c. Select `Create groups` (`folder references` doesn't work) + d. Select `Grid` (containing just the Grid sources, not the entire Git repository) as your target + e. Click `Add` + +3. Add Hadrons source folders to your project + a. As per `Grid`, but change the target to `Hadrons` + b. For each source file in `Archive`, remove them from the target (option-command-1, then untick source files) + +4. Set the following values for each target in `Build Settings` + + Group | Variable | Value + --- | --- | --- + `Deployment` | `DSTROOT` | `$Grid/build$(CONFIGURATION)` *(do this for the entire project)* + `Deployment` | `DEPLOYMENT_LOCATION` | `Yes` + `Deployment` | `INSTALL_PATH` | `$(PRODUCT_NAME)/` + `Deployment` | `SKIP_INSTALL` | `No` + `Linking` | `OTHER_LDFLAGS` | remove `-lGrid -lHadrons` from the list + + This ensures that the libraries are copied back into the build folders when they are made (removing the need to run `make -j 4`) + +5. For `Grid`, in `Build Settings` in the `Build Options` group, set: + + Variable | Configuration | Value + --- | --- | --- + `EXCLUDED_SOURCE_FILE_NAMES` | Non-MPI configurations (`Debug` and `Release`) | `$(Grid)/Grid/communicator/Communicator_mpi3.cc $(Grid)/Grid/communicator/SharedMemoryMPI.cc` + `EXCLUDED_SOURCE_FILE_NAMES` | MPI configurations (`MPIDebug`) | `$(Grid)/Grid/communicator/Communicator_none.cc $(Grid)/Grid/communicator/SharedMemoryNone.cc` + +You should now be able to build and debug any configuration. diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index 15602e93..d1ec1309 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -4,11 +4,12 @@ Source file: ./tests/Test_serialisation.cc - Copyright (C) 2015-2016 + Copyright (C) 2015-2019 Author: Guido Cossu Author: Antonin Portelli Author: Peter Boyle +Author: Michael Marshall This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -28,6 +29,7 @@ Author: Peter Boyle *************************************************************************************/ /* END LEGAL */ #include +#include using namespace Grid; using namespace Grid::QCD; @@ -80,26 +82,159 @@ double d = 2*M_PI; bool b = false; template -void ioTest(const std::string &filename, const O &object, const std::string &name) +void ioTest(const std::string &filename, const O &object, const std::string &name, + const char * tag = "testobject", unsigned short Precision = 0 ) { + std::cout << "IO test: " << name << " -> " << filename << " ..."; // writer needs to be destroyed so that writing physically happens { W writer(filename); - - write(writer, "testobject", object); + if( Precision ) + writer.setPrecision(Precision); + write(writer, tag , object); } + std::cout << " done. reading ..."; R reader(filename); - O buf; - bool good; + std::unique_ptr buf( new O ); // In case object too big for stack - read(reader, "testobject", buf); - good = (object == buf); - std::cout << name << " IO test: " << (good ? "success" : "failure"); - std::cout << std::endl; - if (!good) exit(EXIT_FAILURE); + read(reader, tag, *buf); + bool good = Serializable::CompareMember(object, *buf); + if (!good) { + std::cout << " failure!" << std::endl; + exit(EXIT_FAILURE); + } + std::cout << " done." << std::endl; } +// The only way I could get these iterators to work is to put the begin() and end() functions in the Eigen namespace +// So if Eigen ever defines these, we'll have a conflict and have to change this +namespace Eigen { + template + inline typename std::enable_if::value, typename EigenIO::Traits::scalar_type *>::type + begin( ET & et ) { return reinterpret_cast::scalar_type *>(et.data()); } + template + inline typename std::enable_if::value, typename EigenIO::Traits::scalar_type *>::type + end( ET & et ) { return begin(et) + et.size() * EigenIO::Traits::count; } +} + +// Perform I/O tests on a range of tensor types +// Test coverage: scalars, complex and GridVectors in single, double and default precision +class TensorIO : public Serializable { + using TestScalar = ComplexD; + using SR3 = Eigen::Sizes<9,4,2>; + using SR5 = Eigen::Sizes<5,4,3,2,1>; + using ESO = Eigen::StorageOptions; + using TensorRank3 = Eigen::Tensor; + using TensorR5 = Eigen::TensorFixedSize; + using TensorR5Alt = Eigen::TensorFixedSize; + using Tensor942 = Eigen::TensorFixedSize; + using aTensor942 = std::vector; + using Perambulator = Eigen::Tensor; + using LSCTensor = Eigen::TensorFixedSize>; + + static const Real FlagR; + static const Complex Flag; + static const ComplexF FlagF; + static const TestScalar FlagTS; + static const char * const pszFilePrefix; + + void Init(unsigned short Precision) + { + for( auto &s : Perambulator1 ) s = Flag; + for( auto &s : Perambulator2 ) s = Flag; + for( auto &s : tensorR5 ) s = FlagR; + for( auto &s : tensorRank3 ) s = FlagF; + for( auto &s : tensor_9_4_2 ) s = FlagTS; + for( auto &t : atensor_9_4_2 ) + for( auto &s : t ) s = FlagTS; + for( auto &s : MyLSCTensor ) s = Flag; + } + + // Perform an I/O test for a single Eigen tensor (of any type) + template + static void TestOne(const char * MyTypeName, unsigned short Precision, std::string &filename, + const char * pszExtension, unsigned int &TestNum, + typename EigenIO::Traits::scalar_type Flag, IndexTypes... otherDims) + { + using Traits = EigenIO::Traits; + using scalar_type = typename Traits::scalar_type; + std::unique_ptr pTensor{new T(otherDims...)}; + for( auto &s : * pTensor ) s = Flag; + filename = pszFilePrefix + std::to_string(++TestNum) + "_" + MyTypeName + pszExtension; + ioTest(filename, * pTensor, MyTypeName, MyTypeName); + } + +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TensorIO + , SpinColourVector, spinColourVector + , SpinColourMatrix, spinColourMatrix + , std::vector, DistilParameterNames + , std::vector, DistilParameterValues + , Perambulator, Perambulator1 + , Perambulator, Perambulator2 + , TensorR5, tensorR5 + , TensorRank3, tensorRank3 + , Tensor942, tensor_9_4_2 + , aTensor942, atensor_9_4_2 + , LSCTensor, MyLSCTensor + ); + TensorIO() + : DistilParameterNames {"do", "androids", "dream", "of", "electric", "sheep?"} + , DistilParameterValues{2,3,1,4,5,1} + , Perambulator1(2,3,1,4,5,1) + , Perambulator2(7,1,6,1,5,1) + , tensorRank3(7,3,2) + , atensor_9_4_2(3) {} + +#define TEST_PARAMS( T ) #T, Precision, filename, pszExtension, TestNum + + // Perform a series of I/O tests for Eigen tensors, including a serialisable object + template + static void Test(const char * pszExtension, unsigned short Precision = 0) + { + // Perform a series of tests on progressively more complex tensors + unsigned int TestNum = 0; + std::string filename; + // Rank 1 tensor containing a single integer + using TensorSingle = Eigen::TensorFixedSize>; + TestOne( TEST_PARAMS( TensorSingle ), 7 ); // lucky! + // Rather convoluted way of defining four complex numbers + using TensorSimple = Eigen::Tensor, 6>; + using I = typename TensorSimple::Index; // NB: Never specified, so same for all my test tensors + // Try progressively more complicated tensors + TestOne( TEST_PARAMS( TensorSimple ), FlagTS, 1,1,1,1,1,1 ); + TestOne( TEST_PARAMS( TensorRank3 ), FlagF, 6, 3, 2 ); + TestOne(TEST_PARAMS( Tensor942 ), FlagTS); + TestOne(TEST_PARAMS( LSCTensor ), Flag ); + TestOne(TEST_PARAMS( TensorR5 ), FlagR); + // Now test a serialisable object containing a number of tensors + { + static const char MyTypeName[] = "TensorIO"; + filename = pszFilePrefix + std::to_string(++TestNum) + "_" + MyTypeName + pszExtension; + std::unique_ptr pObj{new TensorIO()}; + pObj->Init(Precision); + ioTest(filename, * pObj, MyTypeName, MyTypeName, Precision); + } + // Stress test. Too large for the XML or text readers and writers! +#ifdef STRESS_TEST + const std::type_info &tw = typeid( WTR_ ); + if( tw == typeid( Hdf5Writer ) || tw == typeid( BinaryWriter ) ) { + using LCMTensor=Eigen::TensorFixedSize,2>,7>,3>, + Eigen::Sizes<2,4,11,10,9>, Eigen::StorageOptions::RowMajor>; + std::cout << "sizeof( LCMTensor ) = " << sizeof( LCMTensor ) / 1024 / 1024 << " MB" << std::endl; + TestOne(TEST_PARAMS( LCMTensor ), Flag); + } +#endif + } +}; + +const Real TensorIO::FlagR {1}; +const Complex TensorIO::Flag {1,-1}; +const ComplexF TensorIO::FlagF {1,-1}; +const TensorIO::TestScalar TensorIO::FlagTS{1,-1}; +const char * const TensorIO::pszFilePrefix = "tensor_"; + template void tensorConvTestFn(GridSerialRNG &rng, const std::string label) { @@ -121,12 +256,12 @@ 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 + GridSerialRNG rng; rng.SeedFixedIntegers(std::vector({42,10,81,9})); - - std::cout << "==== basic IO" << std::endl; + XmlWriter WR("bother.xml"); // test basic type writing @@ -146,7 +281,6 @@ int main(int argc,char **argv) // test serializable class writing myclass obj(1234); // non-trivial constructor std::vector vec; - std::pair pair; std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl; write(WR,"obj",obj); @@ -154,15 +288,15 @@ int main(int argc,char **argv) vec.push_back(obj); vec.push_back(myclass(5678)); vec.push_back(myclass(3838)); - pair = std::make_pair(myenum::red, myenum::blue); write(WR, "objvec", vec); 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::pair pair = std::make_pair(myenum::red, myenum::blue); std::cout << pair << std::endl; // read tests @@ -184,7 +318,15 @@ int main(int argc,char **argv) #ifdef HAVE_HDF5 ioTest("iotest.h5", obj, "HDF5 (object) "); ioTest("iotest.h5", vec, "HDF5 (vector of objects)"); + std::cout << "\n==== detailed Hdf5 tensor tests (Grid::EigenIO)" << std::endl; + TensorIO::Test(".h5"); #endif + std::cout << "\n==== detailed binary tensor tests (Grid::EigenIO)" << std::endl; + TensorIO::Test(".bin"); + std::cout << "\n==== detailed xml tensor tests (Grid::EigenIO)" << std::endl; + TensorIO::Test(".xml", 6); + std::cout << "\n==== detailed text tensor tests (Grid::EigenIO)" << std::endl; + TensorIO::Test(".dat", 5); std::cout << "\n==== vector flattening/reconstruction" << std::endl; typedef std::vector>> vec3d; @@ -227,4 +369,6 @@ int main(int argc,char **argv) tensorConvTest(rng, ColourVector); tensorConvTest(rng, SpinMatrix); tensorConvTest(rng, SpinVector); + + Grid_finalize(); } diff --git a/tests/smearing/Test_smearing.cc b/tests/smearing/Test_smearing.cc new file mode 100644 index 00000000..996671f0 --- /dev/null +++ b/tests/smearing/Test_smearing.cc @@ -0,0 +1,73 @@ +/* + * + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_cg_prec.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); + + LatticeFermion src(&Grid); random(pRNG,src); + RealD nrm = norm2(src); + LatticeFermion result(&Grid); result=zero; + LatticeGaugeField Umu(&Grid); + // SU3::HotConfiguration(pRNG,Umu); + SU3::ColdConfiguration(Umu); + std::vector U(4,&Grid); + + for(int mu=0;mu(Umu,mu); + } + + std::vector site({4,4,0,0}); + src=zero; + SpinColourVector scv; + scv=zero; + scv()(0)(0) = 1.0; + pokeSite(scv,src,site); + + CovariantSmearing::GaussianSmear(U, src, 2.0, 50, Tdir); + + std::cout << src < + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ +/* */ + + +#include +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main(int argc, char **argv) { + + + Grid_init(&argc, &argv); + + // clang-format off + GridCartesian *FGrid_d = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi()); + GridCartesian *FGrid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *FrbGrid_d = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_d); + GridRedBlackCartesian *FrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_f); + // clang-format on + + std::vector fSeeds({1, 2, 3, 4}); + GridParallelRNG fPRNG(FGrid_d); + fPRNG.SeedFixedIntegers(fSeeds); + + // clang-format off + LatticeFermionD src_d(FGrid_d); gaussian(fPRNG, src_d); + LatticeFermionD resultMGD_d(FGrid_d); resultMGD_d = zero; + LatticeFermionD resultMGF_d(FGrid_d); resultMGF_d = zero; + LatticeGaugeFieldD Umu_d(FGrid_d); + +#if 0 + { + FieldMetaData header; + std::string file("./qcdsf.769.00399.lime"); + std::cout < MdagMOpDwc_d(Dwc_d); + MdagMLinearOperator MdagMOpDwc_f(Dwc_f); + + std::cout << GridLogMessage << "**************************************************" << std::endl; + std::cout << GridLogMessage << "Testing single-precision Multigrid for Wilson Clover" << std::endl; + std::cout << GridLogMessage << "**************************************************" << std::endl; + + auto MGPreconDwc_f = createMGInstance(mgParams, levelInfo_f, Dwc_f, Dwc_f); + + MGPreconDwc_f->setup(); + + if(GridCmdOptionExists(argv, argv + argc, "--runchecks")) { + MGPreconDwc_f->runChecks(1e-6); + } + + MixedPrecisionFlexibleGeneralisedMinimalResidual MPFGMRESPREC( + 1.0e-12, 50000, FGrid_f, *MGPreconDwc_f, 100, false); + + std::cout << std::endl << "Starting with a new solver" << std::endl; + MPFGMRESPREC(MdagMOpDwc_d, src_d, resultMGF_d); + + MGPreconDwc_f->reportTimings(); + + if(GridCmdOptionExists(argv, argv + argc, "--docomparison")) { + + std::cout << GridLogMessage << "**************************************************" << std::endl; + std::cout << GridLogMessage << "Testing double-precision Multigrid for Wilson Clover" << std::endl; + std::cout << GridLogMessage << "**************************************************" << std::endl; + + auto MGPreconDwc_d = createMGInstance(mgParams, levelInfo_d, Dwc_d, Dwc_d); + + MGPreconDwc_d->setup(); + + if(GridCmdOptionExists(argv, argv + argc, "--runchecks")) { + MGPreconDwc_d->runChecks(1e-13); + } + + FlexibleGeneralisedMinimalResidual FGMRESPREC(1.0e-12, 50000, *MGPreconDwc_d, 100, false); + + std::cout << std::endl << "Starting with a new solver" << std::endl; + FGMRESPREC(MdagMOpDwc_d, src_d, resultMGD_d); + + MGPreconDwc_d->reportTimings(); + + std::cout << GridLogMessage << "**************************************************" << std::endl; + std::cout << GridLogMessage << "Comparing single-precision Multigrid with double-precision one for Wilson Clover" << std::endl; + std::cout << GridLogMessage << "**************************************************" << std::endl; + + LatticeFermionD diffFullSolver(FGrid_d); + + RealD deviationFullSolver = axpy_norm(diffFullSolver, -1.0, resultMGF_d, resultMGD_d); + + // clang-format off + LatticeFermionF src_f(FGrid_f); precisionChange(src_f, src_d); + LatticeFermionF resMGF_f(FGrid_f); resMGF_f = zero; + LatticeFermionD resMGD_d(FGrid_d); resMGD_d = zero; + // clang-format on + + (*MGPreconDwc_f)(src_f, resMGF_f); + (*MGPreconDwc_d)(src_d, resMGD_d); + + LatticeFermionD diffOnlyMG(FGrid_d); + LatticeFermionD resMGF_d(FGrid_d); + precisionChange(resMGF_d, resMGF_f); + + RealD deviationOnlyPrec = axpy_norm(diffOnlyMG, -1.0, resMGF_d, resMGD_d); + + // clang-format off + std::cout << GridLogMessage << "Absolute difference between FGMRES preconditioned by double and single precicision MG: " << deviationFullSolver << std::endl; + std::cout << GridLogMessage << "Relative deviation between FGMRES preconditioned by double and single precicision MG: " << deviationFullSolver / norm2(resultMGD_d) << std::endl; + std::cout << GridLogMessage << "Absolute difference between one iteration of MG Prec in double and single precision: " << deviationOnlyPrec << std::endl; + std::cout << GridLogMessage << "Relative deviation between one iteration of MG Prec in double and single precision: " << deviationOnlyPrec / norm2(resMGD_d) << std::endl; + // clang-format on + } + + Grid_finalize(); +}