diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp new file mode 100644 index 00000000..c54ad6a6 --- /dev/null +++ b/Hadrons/A2AMatrix.hpp @@ -0,0 +1,97 @@ +#ifndef A2A_Matrix_hpp_ +#define A2A_Matrix_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE + +template +class A2AMatrixIo +{ +public: + A2AMatrixIo(void) = default; + A2AMatrixIo(std::string filename, std::string dataname, + const unsigned int nt, const unsigned int ni, + const unsigned int nj, const unsigned int blockSize); + ~A2AMatrixIo(void) = default; + void initFile(const MetadataType &d); + void saveBlock(const T *data, const unsigned int i, const unsigned int j); +private: + std::string filename_, dataname_; + unsigned int nt_, ni_, nj_, blockSize_; +}; + +template +A2AMatrixIo::A2AMatrixIo(std::string filename, + std::string dataname, + const unsigned int nt, + const unsigned int ni, + const unsigned int nj, + const unsigned int blockSize) +: filename_(filename), dataname_(dataname) +, nt_(nt), ni_(ni), nj_(nj), blockSize_(blockSize) +{} + +template +void A2AMatrixIo::initFile(const MetadataType &d) +{ +#ifdef HAVE_HDF5 + std::vector dim = {static_cast(nt_), + static_cast(ni_), + static_cast(nj_)}, + chunk = {static_cast(nt_), + static_cast(blockSize_), + static_cast(blockSize_)}; + H5NS::DataSpace dataspace(dim.size(), dim.data()); + H5NS::DataSet dataset; + H5NS::DSetCreatPropList plist; + + // create empty file just with metadata + { + Hdf5Writer writer(filename_); + write(writer, dataname_, d); + } + + // create the dataset + Hdf5Writer writer(filename_); + + push(writer, dataname_); + auto &group = writer.getGroup(); + plist.setChunk(chunk.size(), chunk.data()); + dataset = group.createDataSet("data", Hdf5Type::type(), dataspace, plist); +#else + HADRONS_ERROR(Implementation, "all-to-all matrix I/O needs HDF5 library"); +#endif +} + +template +void A2AMatrixIo::saveBlock(const T *data, + const unsigned int i, + const unsigned int j) +{ +#ifdef HAVE_HDF5 + Hdf5Reader reader(filename_); + std::vector count = {nt_, blockSize_, blockSize_}, + offset = {0, static_cast(i), + static_cast(j)}, + stride = {1, 1, 1}, + block = {1, 1, 1}; + H5NS::DataSpace memspace(count.size(), count.data()), dataspace; + H5NS::DataSet dataset; + size_t shift; + + push(reader, dataname_); + auto &group = reader.getGroup(); + dataset = group.openDataSet("data"); + dataspace = dataset.getSpace(); + dataspace.selectHyperslab(H5S_SELECT_SET, count.data(), offset.data(), + stride.data(), block.data()); + dataset.write(data, Hdf5Type::type(), memspace, dataspace); +#else + HADRONS_ERROR(Implementation, "all-to-all matrix I/O needs HDF5 library"); +#endif +} + +END_HADRONS_NAMESPACE + +#endif // A2A_Matrix_hpp_ \ No newline at end of file diff --git a/Hadrons/Makefile.am b/Hadrons/Makefile.am index 93251dca..f10a9bca 100644 --- a/Hadrons/Makefile.am +++ b/Hadrons/Makefile.am @@ -15,7 +15,9 @@ libHadrons_adir = $(includedir)/Hadrons nobase_libHadrons_a_HEADERS = \ $(modules_hpp) \ A2AVectors.hpp \ + A2AMatrix.hpp \ Application.hpp \ + DilutedNoise.hpp \ EigenPack.hpp \ Environment.hpp \ Exceptions.hpp \ diff --git a/Hadrons/Modules/MContraction/A2AMesonField.hpp b/Hadrons/Modules/MContraction/A2AMesonField.hpp index 4ca538e1..2af842dc 100644 --- a/Hadrons/Modules/MContraction/A2AMesonField.hpp +++ b/Hadrons/Modules/MContraction/A2AMesonField.hpp @@ -34,6 +34,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include @@ -62,6 +63,14 @@ class A2AMesonFieldPar: Serializable std::vector, mom); }; +class A2AMesonFieldMetadata: Serializable +{ + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMesonFieldMetadata, + std::vector, momentum, + Gamma::Algebra, gamma); +}; + template class TA2AMesonField : public Module { @@ -70,6 +79,13 @@ public: SOLVER_TYPE_ALIASES(FImpl,); typedef Eigen::TensorMap> MesonField; typedef Eigen::TensorMap> MesonFieldIo; + typedef A2AMatrixIo MatrixIo; + struct IoHelper + { + MatrixIo io; + A2AMesonFieldMetadata metadata; + size_t offset; + }; public: // constructor TA2AMesonField(const std::string name); @@ -86,16 +102,13 @@ private: // IO std::string ioname(unsigned int m, unsigned int g) const; std::string filename(unsigned int m, unsigned int g) const; - void initFile(unsigned int m, unsigned int g); - void saveBlock(const MesonFieldIo &mf, - unsigned int m, unsigned int g, - unsigned int i, unsigned int j); + void saveBlock(MF_IO_TYPE *data, IoHelper &h, unsigned int i, unsigned int j); private: bool hasPhase_{false}; std::string momphName_; std::vector gamma_; std::vector> mom_; - std::vector> nodeFile_; + std::vector nodeIo_; }; MODULE_REGISTER(A2AMesonField, ARG(TA2AMesonField), MContraction); @@ -331,51 +344,42 @@ void TA2AMesonField::execute(void) makeFileDir(filename(0, 0), env().getGrid()); #ifdef MF_PARALLEL_IO env().getGrid()->Barrier(); - nodeFile_.clear(); + nodeIo_.clear(); for(int f = myRank; f < nmom*ngamma; f += nRank) { - std::pair file; + const unsigned int m = f/ngamma, g = f % ngamma; + IoHelper h; - file.first = f/ngamma; - file.second = f % ngamma; - nodeFile_.push_back(file); + h.io = MatrixIo(filename(m, g), ioname(m, g), nt, N_i, N_j, block); + for (auto pmu: mom_[m]) + { + h.metadata.momentum.push_back(pmu); + } + h.metadata.gamma = gamma_[g]; + h.offset = (m*ngamma + g)*nt*block*block; + nodeIo_.push_back(h); } // parallel IO - for (auto &f: nodeFile_) + for (auto &h: nodeIo_) { - auto m = f.first, g = f.second; - - if ((i == 0) and (j == 0)) - { - startTimer("IO: file creation"); - initFile(m, g); - stopTimer("IO: file creation"); - } - startTimer("IO: write block"); - saveBlock(mfBlock, m, g, i, j); - stopTimer("IO: write block"); + saveBlock(mfBlock.data(), h, i, j); } env().getGrid()->Barrier(); #else // serial IO for(int m = 0; m < nmom; m++) for(int g = 0; g < ngamma; g++) - { - if ((i == 0) and (j == 0)) + { + IoHelper h; + + h.io = MatrixIo(filename(m, g), ioname(m, g), nt, N_i, N_j, block); + for (auto pmu: mom_[m]) { - startTimer("IO: file creation"); - if (env().getGrid()->IsBoss()) - { - initFile(m, g); - } - stopTimer("IO: file creation"); + h.metadata.momentum.push_back(pmu); } - startTimer("IO: write block"); - if (env().getGrid()->IsBoss()) - { - saveBlock(mfBlock, m, g, i, j); - } - stopTimer("IO: write block"); + h.metadata.gamma = gamma_[g]; + h.offset = (m*ngamma + g)*nt*block*block; + saveBlock(mfBlock.data(), h, i, j); } #endif stopTimer("IO: total"); @@ -412,71 +416,18 @@ std::string TA2AMesonField::filename(unsigned int m, unsigned int g) cons } template -void TA2AMesonField::initFile(unsigned int m, unsigned int g) -{ -#ifdef HAVE_HDF5 - std::string f = filename(m, g); - auto &v = envGet(std::vector, par().v); - auto &w = envGet(std::vector, par().w); - int nt = env().getDim().back(); - int N_i = w.size(); - int N_j = v.size(); - - Hdf5Writer writer(f); - std::vector dim = {static_cast(nt), - static_cast(N_i), - static_cast(N_j)}, - chunk = {static_cast(nt), - static_cast(par().block), - static_cast(par().block)}; - H5NS::DataSpace dataspace(dim.size(), dim.data()); - H5NS::DataSet dataset; - H5NS::DSetCreatPropList plist; - - push(writer, ioname(m, g)); - write(writer, "momentum", mom_[m]); - write(writer, "gamma", gamma_[g]); - auto &group = writer.getGroup(); - plist.setChunk(chunk.size(), chunk.data()); - dataset = group.createDataSet("mesonField", Hdf5Type::type(), - dataspace, plist); -#else - HADRONS_ERROR(Implementation, "meson field I/O needs HDF5 library"); -#endif -} - -template -void TA2AMesonField::saveBlock(const MesonFieldIo &mf, - unsigned int m, unsigned int g, +void TA2AMesonField::saveBlock(MF_IO_TYPE *data, IoHelper &h, unsigned int i, unsigned int j) { -#ifdef HAVE_HDF5 - std::string f = filename(m, g); - Hdf5Reader reader(f); - hsize_t nt = mf.dimension(2), - Ni = mf.dimension(3), - Nj = mf.dimension(4); - std::vector count = {nt, Ni, Nj}, - offset = {0, static_cast(i), - static_cast(j)}, - stride = {1, 1, 1}, - block = {1, 1, 1}; - H5NS::DataSpace memspace(count.size(), count.data()), dataspace; - H5NS::DataSet dataset; - size_t shift; - - push(reader, ioname(m, g)); - auto &group = reader.getGroup(); - dataset = group.openDataSet("mesonField"); - dataspace = dataset.getSpace(); - dataspace.selectHyperslab(H5S_SELECT_SET, count.data(), offset.data(), - stride.data(), block.data()); - shift = (m*mf.dimension(1) + g)*nt*Ni*Nj; - dataset.write(mf.data() + shift, Hdf5Type::type(), memspace, - dataspace); -#else - HADRONS_ERROR(Implementation, "meson field I/O needs HDF5 library"); -#endif + if ((i == 0) and (j == 0)) + { + startTimer("IO: file creation"); + h.io.initFile(h.metadata); + stopTimer("IO: file creation"); + } + startTimer("IO: write block"); + h.io.saveBlock(data + h.offset, i, j); + stopTimer("IO: write block"); } END_MODULE_NAMESPACE