diff --git a/Hadrons/Modules/MDistil/Distil.hpp b/Hadrons/Modules/MDistil/Distil.hpp index ae68ebcd..4071d5bd 100644 --- a/Hadrons/Modules/MDistil/Distil.hpp +++ b/Hadrons/Modules/MDistil/Distil.hpp @@ -238,7 +238,7 @@ inline GridCartesian * MakeLowerDimGrid( GridCartesian * gridHD ) Perambulator object ******************************************************************************/ -template +template class NamedTensor : public Eigen::Tensor { public: @@ -261,20 +261,26 @@ public: // load and save - not virtual - probably all changes inline void load(const std::string filename); inline void save(const std::string filename) const; - inline void ReadTemporary(const std::string filename); - inline void WriteTemporary(const std::string filename) const; + inline void ReadBinary(const std::string filename); + inline void WriteBinary(const std::string filename); }; /****************************************************************************** - Save NamedTensor binary format + Save NamedTensor binary format (NB: On-disk format is Big Endian) ******************************************************************************/ -template -void NamedTensor::WriteTemporary(const std::string filename) const { +template +void NamedTensor::WriteBinary(const std::string filename) { LOG(Message) << "Writing NamedTensor to \"" << filename << "\"" << std::endl; std::ofstream w(filename, std::ios::binary); + // Number of Scalar_Unit objects per Scalar_ + constexpr unsigned int Scalar_Unit_Size{sizeof(Scalar_Unit)}; + assert((Scalar_Unit_Size == 2 || Scalar_Unit_Size == 4 || Scalar_Unit_Size == 8 ) + && "Scalar_Unit_Size should be 2, 4 or 8"); + assert((sizeof(Scalar_) % Scalar_Unit_Size) == 0 && "Scalar_ is not composed of Scalar_Unit_Size" ); // Size of the data (in bytes) - const std::streamsize TotalDataSize{static_cast(this->size() * sizeof(Scalar_))}; + const auto NumElements{this->size()}; + const std::streamsize TotalDataSize{static_cast(NumElements * sizeof(Scalar_))}; uint64_t u64 = htobe64(static_cast(TotalDataSize)); w.write(reinterpret_cast(&u64), sizeof(u64)); // number of dimensions which aren't 1 @@ -300,7 +306,29 @@ void NamedTensor::WriteTemporary(const std::string filenam d++; } // Actual data - w.write(reinterpret_cast(this->data()), TotalDataSize); + char * const pStart{reinterpret_cast(this->data())}; + // Swap to network byte order in place (alternative is to copy memory - still slow) + void * const pEnd{pStart + TotalDataSize}; + if(Scalar_Unit_Size == 8) + for(uint64_t * p = reinterpret_cast(pStart) ; p < pEnd ; p++ ) + * p = htobe64( * p ); + else if(Scalar_Unit_Size == 4) + for(uint32_t * p = reinterpret_cast(pStart) ; p < pEnd ; p++ ) + * p = htobe32( * p ); + else if(Scalar_Unit_Size == 2) + for(uint16_t * p = reinterpret_cast(pStart) ; p < pEnd ; p++ ) + * p = htobe16( * p ); + w.write(pStart, TotalDataSize); + // Swap back from network byte order + if(Scalar_Unit_Size == 8) + for(uint64_t * p = reinterpret_cast(pStart) ; p < pEnd ; p++ ) + * p = be64toh( * p ); + else if(Scalar_Unit_Size == 4) + for(uint32_t * p = reinterpret_cast(pStart) ; p < pEnd ; p++ ) + * p = be32toh( * p ); + else if(Scalar_Unit_Size == 2) + for(uint16_t * p = reinterpret_cast(pStart) ; p < pEnd ; p++ ) + * p = be16toh( * p ); // checksum #ifdef USE_IPP uint32_t u32 = htobe32(GridChecksum::crc32c(this->data(), TotalDataSize)); @@ -311,15 +339,21 @@ void NamedTensor::WriteTemporary(const std::string filenam } /****************************************************************************** - Load NamedTensor binary format + Load NamedTensor binary format (NB: On-disk format is Big Endian) ******************************************************************************/ -template -void NamedTensor::ReadTemporary(const std::string filename) { +template +void NamedTensor::ReadBinary(const std::string filename) { LOG(Message) << "Reading NamedTensor from \"" << filename << "\"" << std::endl; std::ifstream r(filename, std::ios::binary); + // Number of Scalar_Unit objects per Scalar_ + constexpr unsigned int Scalar_Unit_Size{sizeof(Scalar_Unit)}; + assert((Scalar_Unit_Size == 2 || Scalar_Unit_Size == 4 || Scalar_Unit_Size == 8 ) + && "Scalar_Unit_Size should be 2, 4 or 8"); + assert((sizeof(Scalar_) % Scalar_Unit_Size) == 0 && "Scalar_ is not composed of Scalar_Unit_Size" ); // Size of the data in bytes - const std::streamsize TotalDataSize{static_cast(this->size() * sizeof(Scalar_))}; + const auto NumElements{this->size()}; + const std::streamsize TotalDataSize{static_cast(NumElements * sizeof(Scalar_))}; uint64_t u64; r.read(reinterpret_cast(&u64), sizeof(u64)); assert( TotalDataSize == be64toh( u64 ) && "Error: Size of the data in bytes" ); @@ -350,7 +384,19 @@ void NamedTensor::ReadTemporary(const std::string filename d++; } // Actual data - r.read(reinterpret_cast(this->data()),TotalDataSize); + char * const pStart{reinterpret_cast(this->data())}; + void * const pEnd{pStart + TotalDataSize}; + r.read(pStart,TotalDataSize); + // Swap back from network byte order + if(Scalar_Unit_Size == 8) + for(uint64_t * p = reinterpret_cast(pStart) ; p < pEnd ; p++ ) + * p = be64toh( * p ); + else if(Scalar_Unit_Size == 4) + for(uint32_t * p = reinterpret_cast(pStart) ; p < pEnd ; p++ ) + * p = be32toh( * p ); + else if(Scalar_Unit_Size == 2) + for(uint16_t * p = reinterpret_cast(pStart) ; p < pEnd ; p++ ) + * p = be16toh( * p ); // checksum uint32_t u32; r.read(reinterpret_cast(&u32), sizeof(u32)); @@ -360,15 +406,15 @@ void NamedTensor::ReadTemporary(const std::string filename #else u32 -= GridChecksum::crc32(this->data(), TotalDataSize); #endif - assert( u32 == 0 && "checksum"); + assert( u32 == 0 && "Perambulator checksum invalid"); } /****************************************************************************** Save NamedTensor Hdf5 format ******************************************************************************/ -template -void NamedTensor::save(const std::string filename) const { +template +void NamedTensor::save(const std::string filename) const { LOG(Message) << "Writing NamedTensor to \"" << filename << "\"" << std::endl; #ifndef HAVE_HDF5 LOG(Message) << "Error: I/O for NamedTensor requires HDF5" << std::endl; @@ -382,8 +428,8 @@ void NamedTensor::save(const std::string filename) const { Load NamedTensor Hdf5 format ******************************************************************************/ -template -void NamedTensor::load(const std::string filename) { +template +void NamedTensor::load(const std::string filename) { LOG(Message) << "Reading NamedTensor from \"" << filename << "\"" << std::endl; #ifndef HAVE_HDF5 LOG(Message) << "Error: I/O for NamedTensor requires HDF5" << std::endl; @@ -400,8 +446,8 @@ void NamedTensor::load(const std::string filename) { Perambulator object ******************************************************************************/ -template -using Perambulator = NamedTensor; +template +using Perambulator = NamedTensor; /************************************************************************************* diff --git a/Hadrons/Modules/MDistil/DistilVectors.hpp b/Hadrons/Modules/MDistil/DistilVectors.hpp index c5526bc5..39ba4d09 100644 --- a/Hadrons/Modules/MDistil/DistilVectors.hpp +++ b/Hadrons/Modules/MDistil/DistilVectors.hpp @@ -134,7 +134,7 @@ void TDistilVectors::execute(void) //auto &noise = envGet(std::vector>>, par().noise); auto &noise = envGet(std::vector, par().noise); - auto &perambulator = envGet(Perambulator, par().perambulator); + auto &perambulator = envGet(Perambulator, par().perambulator); auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); auto &rho = envGet(std::vector, getName() + "_rho"); auto &phi = envGet(std::vector, getName() + "_phi"); diff --git a/Hadrons/Modules/MDistil/PerambLight.hpp b/Hadrons/Modules/MDistil/PerambLight.hpp index 3471ce38..7e710d51 100644 --- a/Hadrons/Modules/MDistil/PerambLight.hpp +++ b/Hadrons/Modules/MDistil/PerambLight.hpp @@ -139,7 +139,7 @@ void TPerambLight::setup(void) //envCreate(std::complex, getName() + "_debug_delete_me_3", 1, z); //envCreate(std::complex, getName() + "_debug_delete_me_4", 1, {0.6 COMMA -3.1}); //envCreate(std::array, getName() + "_debug_delete_me_5", 1, {"One" COMMA "Two" COMMA "Three"}); - envCreate(Perambulator, getName() + "_perambulator_light", 1, + envCreate(Perambulator, getName() + "_perambulator_light", 1, sIndexNames,Distil.Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI); envCreate(std::vector, getName() + "_noise", 1, nvec*Distil.Ns*Distil.Nt*Distil.nnoise); @@ -195,7 +195,7 @@ void TPerambLight::execute(void) //auto &noise = envGet(std::vector>>, par().noise); auto &noise = envGet(std::vector, getName() + "_noise"); - auto &perambulator = envGet(Perambulator, getName() + "_perambulator_light"); + auto &perambulator = envGet(Perambulator, getName() + "_perambulator_light"); auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); auto &unsmeared_sink = envGet(std::vector, getName() + "_unsmeared_sink"); @@ -268,7 +268,7 @@ void TPerambLight::execute(void) bExists = true; } if( bExists ) { - perambulator.ReadTemporary(PerambFileName); + perambulator.ReadBinary(PerambFileName); return; } } @@ -368,7 +368,7 @@ void TPerambLight::execute(void) // THIS IS WHERE WE WANT TO SAVE THE PERAMBULATORS TO DISK if(PerambFileName.length()) - perambulator.WriteTemporary(PerambFileName); + perambulator.WriteBinary(PerambFileName); } END_MODULE_NAMESPACE diff --git a/tests/hadrons/Test_hadrons_distil.cc b/tests/hadrons/Test_hadrons_distil.cc index eaa2ef19..34c27df3 100644 --- a/tests/hadrons/Test_hadrons_distil.cc +++ b/tests/hadrons/Test_hadrons_distil.cc @@ -254,7 +254,7 @@ bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true ) #ifdef DEBUG -typedef Grid::Hadrons::MDistil::NamedTensor MyTensor; +typedef Grid::Hadrons::MDistil::NamedTensor MyTensor; void DebugShowTensor(MyTensor &x, const char * n) { @@ -270,7 +270,7 @@ void DebugShowTensor(MyTensor &x, const char * n) MyTensor::Index SizeCalculated{1}; std::cout << "Dimensions again"; for(int i=0 ; i < d.size() ; i++ ) { - std::cout << " : [" << i << "]=" << d[i]; + std::cout << " : [" << i << ", " << x.IndexNames[i] << "]=" << d[i]; SizeCalculated *= d[i]; } std::cout << std::endl; @@ -281,7 +281,7 @@ void DebugShowTensor(MyTensor &x, const char * n) for( int i = 0 ; i < d[0] ; i++ ) for( int j = 0 ; j < d[1] ; j++ ) for( int k = 0 ; k < d[2] ; k++ ) { - x(i,j,k) = std::complex(SizeCalculated, SizeCalculated); + x(i,j,k) = std::complex(SizeCalculated, -SizeCalculated); SizeCalculated--; } // Show raw data @@ -300,19 +300,20 @@ bool DebugEigenTest() std::array as={"Alpha", "Beta", "Gamma"}; MyTensor x(as, 2,1,4); DebugShowTensor(x, "x"); - x.WriteTemporary(pszTestFileName); + x.WriteBinary(pszTestFileName); + DebugShowTensor(x, "x"); // Test initialisation of an array of strings for( auto a : as ) std::cout << a << std::endl; - Grid::Hadrons::MDistil::Perambulator p{as,2,7,2}; + Grid::Hadrons::MDistil::Perambulator p{as,2,7,2}; DebugShowTensor(p, "p"); std::cout << "p.IndexNames follow" << std::endl; for( auto a : p.IndexNames ) std::cout << a << std::endl; // Now see whether we can read a tensor back - std::array a2={"Alpha", "Delta", "Gamma"}; - MyTensor y(a2, 2,1,4); - y.ReadTemporary(pszTestFileName); + std::array a2={"Alpha", "Gamma", "Delta"}; + MyTensor y(a2, 2,4,1); + y.ReadBinary(pszTestFileName); DebugShowTensor(y, "y"); return true; }