diff --git a/Grid/parallelIO/IldgIO.h b/Grid/parallelIO/IldgIO.h index 688d1ac9b8..fc47fc4e52 100644 --- a/Grid/parallelIO/IldgIO.h +++ b/Grid/parallelIO/IldgIO.h @@ -4,7 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid Source file: ./lib/parallelIO/IldgIO.h -Copyright (C) 2015 +Copyright (C) 2015, 2026 + +Author: Peter Boyle +Author: Guido Cossu +Author: Gaurav Ray 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 @@ -413,7 +417,7 @@ class GridLimeWriter : public BinaryIO // This routine is Collectively called by all nodes // in communicator used by the field.Grid() //////////////////////////////////////////////////// - template + template void writeLimeLatticeBinaryObject(Lattice &field,std::string record_name,int control=BINARYIO_LEXICOGRAPHIC) { //////////////////////////////////////////////////////////////////// @@ -438,9 +442,9 @@ class GridLimeWriter : public BinaryIO //////////////////////////////////////////// // Create record header //////////////////////////////////////////// - typedef typename vobj::scalar_object sobj; int err; uint32_t nersc_csum,scidac_csuma,scidac_csumb; + typedef typename GaugeUnMunger::out_type sobj; uint64_t PayloadSize = sizeof(sobj) * grid->_gsites; if ( boss_node ) { createLimeRecordHeader(record_name, 0, 0, PayloadSize); @@ -463,9 +467,11 @@ class GridLimeWriter : public BinaryIO /////////////////////////////////////////// // The above is collective. Write by other means into the binary record /////////////////////////////////////////// - std::string format = getFormatString(); - BinarySimpleMunger munge; - BinaryIO::writeLatticeObject(field, filename, munge, offset1, format,nersc_csum,scidac_csuma,scidac_csumb,control); + std::string format = getFormatString(); + + GaugeUnMunger unmunger; + + BinaryIO::writeLatticeObject(field, filename, unmunger, offset1, format, nersc_csum, scidac_csuma, scidac_csumb, control); /////////////////////////////////////////// // Wind forward and close the record @@ -484,7 +490,7 @@ class GridLimeWriter : public BinaryIO err=limeWriterCloseRecord(LimeW); assert(err>=0); } //////////////////////////////////////// - // Write checksum element, propagaing forward from the BinaryIO + // Write checksum element, propagating forward from the BinaryIO // Always pair a checksum with a binary object, and close message //////////////////////////////////////// scidacChecksum checksum; @@ -620,7 +626,7 @@ class IldgWriter : public ScidacWriter { { uint64_t PayloadSize = LFN.size(); int err; - createLimeRecordHeader(ILDG_DATA_LFN, 0 , 0, PayloadSize); + createLimeRecordHeader(ILDG_DATA_LFN, 1 , 1, PayloadSize); err=limeWriteRecordData(const_cast(LFN.c_str()), &PayloadSize,LimeW); assert(err>=0); err=limeWriterCloseRecord(LimeW); assert(err>=0); } @@ -630,13 +636,13 @@ class IldgWriter : public ScidacWriter { // Don't require scidac records EXCEPT checksum // Use Grid MetaData object if present. //////////////////////////////////////////////////////////////// - template - void writeConfiguration(Lattice &Umu,int sequence,std::string LFN,std::string description) + template + void writeConfiguration(Lattice &Umu, int sequence, std::string LFN, std::string description) { GridBase * grid = Umu.Grid(); - typedef Lattice GaugeField; - typedef vLorentzColourMatrixD vobj; - typedef typename vobj::scalar_object sobj; + + // catch malformed (wrt group_name) template instantiations at compile-time + static_assert( std::is_same_v || (std::is_same_v && Nc%2==0), "IldgWriter supports SU(Nc) and Sp(Nc=2k) lattices. For Sp fields Nc must be even" ); //////////////////////////////////////// // fill the Grid header @@ -649,29 +655,51 @@ class IldgWriter : public ScidacWriter { stats Stats; Stats(Umu,header); - - std::string format = header.floating_point; + header.ensemble_id = description; header.ensemble_label = description; header.sequence_number = sequence; header.ildg_lfn = LFN; - - assert ( (format == std::string("IEEE32BIG")) - ||(format == std::string("IEEE64BIG")) ); - ////////////////////////////////////////////////////// // Fill ILDG header data struct ////////////////////////////////////////////////////// ildgFormat ildgfmt ; const std::string stNC = std::to_string( Nc ) ; - ildgfmt.field = std::string("su"+stNC+"gauge"); - if ( format == std::string("IEEE32BIG") ) { + // use the gauge group to populate the 'field' element of ildg header + if constexpr ( is_su::value ) { + std::cout << GridLogMessage << "writing SU(" << stNC << ") field" << std::endl; + ildgfmt.field = std::string("su"+stNC+"gauge"); + } else if constexpr ( is_sp::value ) { + std::cout << GridLogMessage << "writing Sp(" << stNC << ") field" << std::endl; + ildgfmt.field = std::string("sp"+stNC+"gauge"); + } else { + static_assert(1, "Unrecognised group; unable to determine field tag"); + } + + // populate 'rows' element of ildg header + if constexpr( matrix_fmt==MatrixFormat::REDUCED && is_su::value ) { + ildgfmt.rows = Nc-1 ; + } else if constexpr( matrix_fmt==MatrixFormat::REDUCED && is_sp::value ) { + ildgfmt.rows = Nc/2 ; + } else if constexpr( matrix_fmt==MatrixFormat::FULL ) { + ildgfmt.rows = Nc; + } else { + static_assert(1, "Unknown MatrixFormat specified"); + } + + // set fmt string for single/double precision + if constexpr( fp_fmt == FloatingPointFormat::IEEE32BIG ) { + header.floating_point = std::string("IEEE32BIG"); ildgfmt.precision = 32; - } else { + } else if constexpr ( fp_fmt == FloatingPointFormat::IEEE64BIG ) { + header.floating_point = std::string("IEEE64BIG"); ildgfmt.precision = 64; + } else { + static_assert(1, "Unknown FloatingPointFormat specified"); } - ildgfmt.version = 1.0; + + ildgfmt.version = 1.2; ildgfmt.lx = header.dimension[0]; ildgfmt.ly = header.dimension[1]; ildgfmt.lz = header.dimension[2]; @@ -704,8 +732,8 @@ class IldgWriter : public ScidacWriter { writeLimeObject(1,0,_scidacRecord,_scidacRecord.SerialisableClassName(),std::string(SCIDAC_PRIVATE_RECORD_XML)); writeLimeObject(0,0,info,info.SerialisableClassName(),std::string(SCIDAC_RECORD_XML)); writeLimeObject(0,0,ildgfmt,std::string("ildgFormat") ,std::string(ILDG_FORMAT)); // rec + writeLimeLatticeBinaryObject(Umu,std::string(ILDG_BINARY_DATA)); // Closes message with checksum writeLimeIldgLFN(header.ildg_lfn); // rec - writeLimeLatticeBinaryObject(Umu,std::string(ILDG_BINARY_DATA)); // Closes message with checksum // limeDestroyWriter(LimeW); } }; @@ -713,6 +741,62 @@ class IldgWriter : public ScidacWriter { class IldgReader : public GridLimeReader { public: + ////////////////////////////////////////////////////////////////// + // this helper function wraps the logic for choosing + // the correct munger and intermediate lattice data type + // when reading a lattice field from disk. + // This is a runtime choice as Grid has to first read the + // header of the lattice cfg. Therefore templating this function on + // FloatingPointFormat etc. does not work. It is a rather + // cumbersome if statement with 6 branches, + // but the benefit of organising IldgReader this way is it makes + // readConfiguration clearer and more readable. + ////////////////////////////////////////////////////////////////// + template + void readLatticeBinaryObject(Lattice &Umu, std::string filename, FloatingPointFormat fp_fmt, MatrixFormat matrix_fmt, bool is_grp_su, bool is_grp_sp, uint64_t &offset, uint32_t &nersc_csum, uint32_t &scidac_csuma, uint32_t &scidac_csumb) + { + + // need all the types we could possibly read from + // including the intermediate data types for + // reduced format lattices and single/double precision + typedef typename vobj::scalar_object sobj; + typedef LorentzColourMatrixF fobj; + typedef LorentzColourMatrixD dobj; + typedef LorentzColour2x3F fobjsuR; + typedef LorentzColour2x3D dobjsuR; + typedef LorentzColourNx2NF fobjspR; + typedef LorentzColourNx2ND dobjspR; + + std::string format; + + if ( fp_fmt == FloatingPointFormat::IEEE64BIG ) { + format = std::string("IEEE64BIG"); + if ( is_grp_su && matrix_fmt==MatrixFormat::REDUCED ) { + GaugeSUmunger munge; + BinaryIO::readLatticeObject(Umu, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); + } else if ( is_grp_sp && matrix_fmt==MatrixFormat::REDUCED ) { + GaugeSpmunger munge; + BinaryIO::readLatticeObject(Umu, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); + } else { + GaugeSimpleMunger munge; + BinaryIO::readLatticeObject(Umu, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); + } + } else if ( fp_fmt == FloatingPointFormat::IEEE32BIG) { + format = std::string("IEEE32BIG"); + if ( is_grp_su && matrix_fmt==MatrixFormat::REDUCED ) { + GaugeSUmunger munge; + BinaryIO::readLatticeObject(Umu, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); + } else if ( is_grp_sp && matrix_fmt==MatrixFormat::REDUCED ) { + GaugeSpmunger munge; + BinaryIO::readLatticeObject(Umu, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); + } else { + GaugeSimpleMunger munge; + BinaryIO::readLatticeObject(Umu, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); + } + } else { assert("Unable to determine which readLatticeObject function template to instantiate."); } + + } + //////////////////////////////////////////////////////////////// // Read either Grid/SciDAC/ILDG configuration // Don't require scidac records EXCEPT checksum @@ -720,15 +804,8 @@ class IldgReader : public GridLimeReader { // Else use ILDG MetaData object if present. // Else use SciDAC MetaData object if present. //////////////////////////////////////////////////////////////// - template - void readConfiguration(Lattice &Umu, FieldMetaData &FieldMetaData_) { - - typedef Lattice GaugeField; - typedef typename GaugeField::vector_object vobj; - typedef typename vobj::scalar_object sobj; - - typedef LorentzColourMatrixF fobj; - typedef LorentzColourMatrixD dobj; + template + void readConfiguration(Lattice &Umu, FieldMetaData &FieldMetaData_) { GridBase *grid = Umu.Grid(); @@ -755,8 +832,15 @@ class IldgReader : public GridLimeReader { uint32_t scidac_csuma; uint32_t scidac_csumb; + // these variables store information about the lattice that is read + // from its ildg-format header. if matrix_fmt==MatrixFormat::REDUCED + // then Grid will reconstruct the full matrix using the appropriate munger. + bool is_grp_su = false; + bool is_grp_sp = false; + MatrixFormat matrix_fmt; // Binary format std::string format; + FloatingPointFormat fp_fmt; ////////////////////////////////////////////////////////////////////////// // Loop over all records @@ -787,11 +871,52 @@ class IldgReader : public GridLimeReader { std::string xmlstring(&xmlc[0]); if ( !strncmp(limeReaderType(LimeR), ILDG_FORMAT,strlen(ILDG_FORMAT)) ) { + // 'older' format ildg lattices don't have a element in their header. + // Grid's XML parsing doesn't support optional parameters with default values. + // Rather than rewrite that, here we explicitly add the element with default + // value = Nc to the ildg header if it is missing, before passing it to XmlReader." + pugi::xml_document doc; + doc.load_string(xmlstring.c_str()); + + if(doc.child("ildgFormat").child("rows")) { + std::string rows = doc.child("ildgFormat").child("rows").child_value(); + std::cout << GridLogMessage << " element present = " << rows << ". So this lattice might be ildg 1.2 compliant." << std::endl; + } else { + std::cout << GridLogMessage << " element not present - adding it after ." << std::endl; + pugi::xml_node ildgfmt = doc.child("ildgFormat"); + const std::string stNC = std::to_string(Nc); + ildgfmt.insert_child_after("rows", ildgfmt.child("field")).text().set(stNC.c_str()); + + // write result back into xmlstring + std::ostringstream ss; + doc.save(ss); + xmlstring = ss.str(); + } + XmlReader RD(xmlstring, true, ""); read(RD,"ildgFormat",ildgFormat_); - if ( ildgFormat_.precision == 64 ) format = std::string("IEEE64BIG"); - if ( ildgFormat_.precision == 32 ) format = std::string("IEEE32BIG"); + if ( ildgFormat_.precision == 64 ) { fp_fmt = FloatingPointFormat::IEEE64BIG; } + if ( ildgFormat_.precision == 32 ) { fp_fmt = FloatingPointFormat::IEEE32BIG; } + + // set vars to tell Grid if it needs to reconstruct the full field. + std::cout << GridLogMessage << "ildgFormat_.rows is " << ildgFormat_.rows << std::endl; + matrix_fmt = (ildgFormat_.rows < Nc) ? MatrixFormat::REDUCED : MatrixFormat::FULL; + if( !strncmp(ildgFormat_.field.c_str(),"su",2) ) { is_grp_su = true; } + if( !strncmp(ildgFormat_.field.c_str(),"sp",2) ) { is_grp_sp = true; } + // check if field element corresponds to either su or sp + assert( is_grp_su || is_grp_sp ); + + // check rows and Nc are related in the way we expect for each gauge group. + if (matrix_fmt==MatrixFormat::REDUCED && is_grp_su) { + assert( ildgFormat_.rows == Nc-1 ); + } + if (matrix_fmt==MatrixFormat::REDUCED && is_grp_sp) { + assert( ildgFormat_.rows == Nc/2 ); + } + if (matrix_fmt==MatrixFormat::FULL) { + assert( ildgFormat_.rows == Nc ); + } assert( ildgFormat_.lx == dims[0]); assert( ildgFormat_.ly == dims[1]); @@ -848,14 +973,10 @@ class IldgReader : public GridLimeReader { // Binary data ///////////////////////////////// // std::cout << GridLogMessage << "ILDG Binary record found : " ILDG_BINARY_DATA << std::endl; + uint64_t offset= ftello(File); - if ( format == std::string("IEEE64BIG") ) { - GaugeSimpleMunger munge; - BinaryIO::readLatticeObject< vobj, dobj >(Umu, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); - } else { - GaugeSimpleMunger munge; - BinaryIO::readLatticeObject< vobj, fobj >(Umu, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); - } + + readLatticeBinaryObject(Umu, filename, fp_fmt, matrix_fmt, is_grp_su, is_grp_sp, offset, nersc_csum, scidac_csuma, scidac_csumb); found_ildgBinary = 1; } @@ -876,7 +997,7 @@ class IldgReader : public GridLimeReader { if ( found_FieldMetaData ) { - std::cout << GridLogMessage<<"Grid MetaData was record found: configuration was probably written by Grid ! Yay ! "< Author: Jamie Hudspith + Author: Gaurav Ray 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 @@ -209,7 +210,7 @@ inline void reconstruct3(LorentzColourMatrix & cm) cm(mu)()(1,0) = -adj(cm(mu)()(0,y)) ; cm(mu)()(1,1) = adj(cm(mu)()(0,x)) ; #else - const int x=0 , y=1 , z=2 ; // a little disinenuous labelling + const int x=0 , y=1 , z=2 ; // a little disingenuous labelling cm(mu)()(2,x) = adj(cm(mu)()(0,y)*cm(mu)()(1,z)-cm(mu)()(0,z)*cm(mu)()(1,y)); //x= yz-zy cm(mu)()(2,y) = adj(cm(mu)()(0,z)*cm(mu)()(1,x)-cm(mu)()(0,x)*cm(mu)()(1,z)); //y= zx-xz cm(mu)()(2,z) = adj(cm(mu)()(0,x)*cm(mu)()(1,y)-cm(mu)()(0,y)*cm(mu)()(1,x)); //z= xy-yx @@ -217,15 +218,143 @@ inline void reconstruct3(LorentzColourMatrix & cm) } } +// the elements of the final row of an SU(N) matrix +// can be obtained from the determinants of the N N-1 by N-1 matrices +// formed by deleting each column and the Nth row of the SU(N) matrix. +// 1) create a ColourSubMatrix object to store the N-1 dim matrix. +// 2) peek the colour matrix in a particular direction. +// 3) fill the ColourSubMatrix object from the peeked matrix. +// 4) take the Determinant and fill the corresponding element +// 5) repeat for each Lorentz index +inline void reconstructSU(LorentzColourMatrix &cm) +{ + using ColourSubMatrix = iScalar > > ; + + ColourSubMatrix Su; // for the Nc-1 by Nc-1 matrix + ColourMatrix SU; // for the Nc-1 by N matrix read from disk + + for(int mu=0; mu(cm,mu); + for(int k=0; k(cm,SU,mu); + } +} + +// this function determines if a sequence of integers {i0...ik} +// is an even or odd parity permutation. even/odd if the number of +// inversions needed to get back to the lexicographic 1st sequence is +// even or odd. ex: {1,2,0} --> {1,0,2} --> {0,1,2} implies {1,2,0} is even. +inline bool is_perm_even(std::vector &v) { + + int n = v.size(); + std::vector a(n,0); + int c = 0; + + for(int j=0; j indices(Nc); + std::iota(indices.begin(), indices.end(), 0); // {0,1,...,Nc-1} + + for(int mu=0; mu cols = indices; + cols.erase(cols.begin()+k); // {0,...,k-1,k+1,...Nc-1} + cm(mu)()(Nc-1,k) = 0; + // construct sum of products + // as we cycle through permutations of cols + do + { + std::vector perm = cols; + perm.emplace_back(k); + ComplexD prod(1,0); + for(int i=0;i using iLorentzColour2x3 = iVector, Nc-1>, Nd >; - typedef iLorentzColour2x3 LorentzColour2x3; typedef iLorentzColour2x3 LorentzColour2x3F; typedef iLorentzColour2x3 LorentzColour2x3D; +// intermediate types for Sp(2N) fields +template using iLorentzColourNx2N = iVector, Nc/2>, Nd >; +typedef iLorentzColourNx2N LorentzColourNx2NF; +typedef iLorentzColourNx2N LorentzColourNx2ND; + ///////////////////////////////////////////////////////////////////////////////// // Simple classes for precision conversion ///////////////////////////////////////////////////////////////////////////////// @@ -273,8 +402,8 @@ struct GaugeSimpleMunger{ void operator()(fobj &in, sobj &out) { for (int mu = 0; mu < Nd; mu++) { for (int i = 0; i < Nc; i++) { - for (int j = 0; j < Nc; j++) { - out(mu)()(i, j) = in(mu)()(i, j); + for (int j = 0; j < Nc; j++) { + out(mu)()(i, j) = in(mu)()(i, j); }} } }; @@ -341,5 +470,148 @@ struct Gauge3x2unmunger{ } }; +template +struct GaugeSUmunger; +// two specialisations for two ways to reconstruct +// NcxNc matrix from (Nc-1)xNc matrix + +template +struct GaugeSUmunger{ + void operator() (fobj &in,sobj &out){ + for(int mu=0;mu +struct GaugeSUmunger{ + void operator() (fobj &in,sobj &out){ + for(int mu=0;mu +struct GaugeSUunmunger{ + void operator() (sobj &in,fobj &out){ + for(int mu=0;mu +struct GaugeSpmunger{ + void operator() (fobj &in,sobj &out){ + for(int mu=0;mu +struct GaugeSpunmunger{ + void operator() (sobj &in,fobj &out){ + for(int mu=0;mu no group reduction - treat SU and Sp the same +// > group reduction for SU fields +// > group reduction for Sp fields +// It also exposes the intermediate out_type for use in +// IldgWriter.writeConfiguration +///////////////////////////////////////////////////////// +template +struct GaugeUnMunger; + +// no group reduction +template +struct GaugeUnMunger +{ + using in_type = typename vobj::scalar_object; + using out_type = typename std::tuple_element_t(fp_fmt), std::tuple>; + BinarySimpleUnmunger unmunger; + + void operator() (in_type &in, out_type &out){ + unmunger(in,out); + } +}; + +//template specialisation for SU +template +struct GaugeUnMunger +{ + using in_type = typename vobj::scalar_object; + using tmp_type = typename std::tuple_element_t(fp_fmt), std::tuple>; + using out_type = typename std::tuple_element_t(fp_fmt), std::tuple>; + + BinarySimpleUnmunger binary_unmunger; + GaugeSUunmunger gauge_unmunger; + + void operator() (in_type &in, out_type &out){ + tmp_type tmp; + binary_unmunger(in, tmp); + gauge_unmunger(tmp,out); + } +}; + +//template specialisation for Sp +template +struct GaugeUnMunger +{ + using in_type = typename vobj::scalar_object; + using tmp_type = typename std::tuple_element_t(fp_fmt), std::tuple>; + using out_type = typename std::tuple_element_t(fp_fmt), std::tuple>; + + BinarySimpleUnmunger binary_unmunger; + GaugeSpunmunger gauge_unmunger; + + void operator() (in_type &in, out_type &out){ + tmp_type tmp; + binary_unmunger(in, tmp); + gauge_unmunger(tmp, out); + } +}; + + NAMESPACE_END(Grid); diff --git a/Grid/qcd/hmc/checkpointers/BaseCheckpointer.h b/Grid/qcd/hmc/checkpointers/BaseCheckpointer.h index 72f5e39bab..03972a8561 100644 --- a/Grid/qcd/hmc/checkpointers/BaseCheckpointer.h +++ b/Grid/qcd/hmc/checkpointers/BaseCheckpointer.h @@ -39,15 +39,21 @@ class CheckpointerParameters : Serializable { std::string, rng_prefix, int, saveInterval, bool, saveSmeared, - std::string, format, ); + std::string, format, + std::string, group, + bool, reduced_matrix, + bool, unique_su, ); - CheckpointerParameters(std::string cf = "cfg", std::string sf="cfg_smr" , std::string rn = "rng", - int savemodulo = 1, const std::string &f = "IEEE64BIG") + CheckpointerParameters(std::string cf = "cfg", std::string sf="cfg_smr" , std::string rn = "rng", int savemodulo = 1, bool save_smr = false, const std::string& f = "IEEE64BIG", std::string grp = "su", bool rdc_mat = false, bool unq_su = false) : config_prefix(cf), smeared_prefix(sf), rng_prefix(rn), saveInterval(savemodulo), - format(f){}; + saveSmeared(save_smr), + format(f), + group(grp), + reduced_matrix(rdc_mat), + unique_su(unq_su) {}; template diff --git a/Grid/qcd/hmc/checkpointers/ILDGCheckpointer.h b/Grid/qcd/hmc/checkpointers/ILDGCheckpointer.h index 50cca63d17..80e413ba11 100644 --- a/Grid/qcd/hmc/checkpointers/ILDGCheckpointer.h +++ b/Grid/qcd/hmc/checkpointers/ILDGCheckpointer.h @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: ./lib/qcd/hmc/ILDGCheckpointer.h -Copyright (C) 2016 +Copyright (C) 2016, 2026 Author: Guido Cossu +Author: Gaurav Ray 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 @@ -54,20 +55,75 @@ class ILDGHmcCheckpointer : public BaseHmcCheckpointer { // check here that the format is valid int ieee32big = (Params.format == std::string("IEEE32BIG")); - int ieee32 = (Params.format == std::string("IEEE32")); int ieee64big = (Params.format == std::string("IEEE64BIG")); - int ieee64 = (Params.format == std::string("IEEE64")); - if (!(ieee64big || ieee32 || ieee32big || ieee64)) { + if (!(ieee64big || ieee32big)) { std::cout << GridLogError << "Unrecognized file format " << Params.format << std::endl; std::cout << GridLogError - << "Allowed: IEEE32BIG | IEEE32 | IEEE64BIG | IEEE64" + << "Allowed: IEEE32BIG | IEEE64BIG" << std::endl; + exit(1); + } + if ( !((Params.group == std::string("su")) || (Params.group == std::string("sp"))) ) { + std::cout << GridLogError << "Unrecognized gauge group " << Params.group << std::endl; + std::cout << GridLogError << "Allowed: su | sp" << std::endl; exit(1); } + + if ( Params.group == std::string("sp") && Nc%2!=0 ) { + std::cout << GridLogError << "Nc=" << Nc; + std::cout << ", Sp fields require even Nc" << std::endl; + exit(1); + } + + } + + // choose appropriate template instantiation here since the + // non-const checkpointer parameters cannot be used directly as + // template arguments to IldgWriter + void chooseIldgWriter( std::string format, std::string group, bool reduced_matrix, + std::string lat_obj, int traj, + ConfigurationBase &SmartConfig) { + + GridBase *grid = SmartConfig.get_U(false).Grid(); + + IldgWriter _IldgWriter(grid->IsBoss()); + _IldgWriter.open(lat_obj); + + if(format=="IEEE64BIG") { + if(group=="su" && reduced_matrix) { + _IldgWriter.writeConfiguration(SmartConfig.get_U(false), traj, lat_obj, lat_obj); + } + else if (group=="su" && !reduced_matrix) { + _IldgWriter.writeConfiguration(SmartConfig.get_U(false), traj, lat_obj, lat_obj); + } + else if (group=="sp" && reduced_matrix) { + _IldgWriter.writeConfiguration(SmartConfig.get_U(false), traj, lat_obj, lat_obj); + } + else if (group=="sp" && !reduced_matrix) { + _IldgWriter.writeConfiguration(SmartConfig.get_U(false), traj, lat_obj, lat_obj); + } + } + else if (format=="IEEE32BIG") { + if(group=="su" && reduced_matrix) { + _IldgWriter.writeConfiguration(SmartConfig.get_U(false), traj, lat_obj, lat_obj); + } + else if (group=="su" && !reduced_matrix) { + _IldgWriter.writeConfiguration(SmartConfig.get_U(false), traj, lat_obj, lat_obj); + } + else if (group=="sp" && reduced_matrix) { + _IldgWriter.writeConfiguration(SmartConfig.get_U(false), traj, lat_obj, lat_obj); + } + else if (group=="sp" && !reduced_matrix) { + _IldgWriter.writeConfiguration(SmartConfig.get_U(false), traj, lat_obj, lat_obj); + } + } + + _IldgWriter.close(); } + void TrajectoryComplete(int traj, ConfigurationBase &SmartConfig, @@ -76,7 +132,6 @@ class ILDGHmcCheckpointer : public BaseHmcCheckpointer { if ((traj % Params.saveInterval) == 0) { std::string config, rng, smr; this->build_filenames(traj, Params, config, smr, rng); - GridBase *grid = SmartConfig.get_U(false).Grid(); uint32_t nersc_csum,scidac_csuma,scidac_csumb; BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb); std::cout << GridLogMessage << "Written BINARY RNG " << rng @@ -86,11 +141,8 @@ class ILDGHmcCheckpointer : public BaseHmcCheckpointer { << scidac_csumb << std::dec << std::endl; - - IldgWriter _IldgWriter(grid->IsBoss()); - _IldgWriter.open(config); - _IldgWriter.writeConfiguration(SmartConfig.get_U(false), traj, config, config); - _IldgWriter.close(); + chooseIldgWriter( Params.format, Params.group, Params.reduced_matrix, config, traj, + SmartConfig ); std::cout << GridLogMessage << "Written ILDG Configuration on " << config << " checksum " << std::hex @@ -100,17 +152,15 @@ class ILDGHmcCheckpointer : public BaseHmcCheckpointer { << std::dec << std::endl; if ( Params.saveSmeared ) { - IldgWriter _IldgWriter(grid->IsBoss()); - _IldgWriter.open(smr); - _IldgWriter.writeConfiguration(SmartConfig.get_U(true), traj, smr, smr); - _IldgWriter.close(); - - std::cout << GridLogMessage << "Written ILDG Configuration on " << smr - << " checksum " << std::hex - << nersc_csum<<"/" - << scidac_csuma<<"/" - << scidac_csumb - << std::dec << std::endl; + chooseIldgWriter( Params.format, Params.group, Params.reduced_matrix, smr, traj, + SmartConfig ); + + std::cout << GridLogMessage << "Written ILDG Configuration on " << smr + << " checksum " << std::hex + << nersc_csum<<"/" + << scidac_csuma<<"/" + << scidac_csumb + << std::dec << std::endl; } } @@ -123,15 +173,18 @@ class ILDGHmcCheckpointer : public BaseHmcCheckpointer { this->check_filename(rng); this->check_filename(config); - - uint32_t nersc_csum,scidac_csuma,scidac_csumb; BinaryIO::readRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb); FieldMetaData header; IldgReader _IldgReader; _IldgReader.open(config); - _IldgReader.readConfiguration(U,header); // format from the header + + if(Params.unique_su) { + _IldgReader.readConfiguration(U,header);// format from the header + } + else { _IldgReader.readConfiguration(U,header); } + _IldgReader.close(); std::cout << GridLogMessage << "Read ILDG Configuration from " << config diff --git a/documentation/manual.rst b/documentation/manual.rst index 0c463672c1..b42bd42ce6 100644 --- a/documentation/manual.rst +++ b/documentation/manual.rst @@ -2104,21 +2104,24 @@ These are specialised to SciDAC writers, introducing facilities for generating t void readScidacFieldRecord(Lattice &field,userRecord &_userRecord); }; -They are also specialised to ILDG format writers, available and defined only for Gauge configurations:: +They are also specialised to ILDG format writers, available and defined only for Gauge configurations. These ``IldgWriter`` and ``IldgReader`` template classes can read and write ILDG-compliant SU(N) and Sp(2N) lattice configurations as full NxN and 2Nx2N matrix fields respectively. They also support reduced format read/write with (N-1)xN matrix fields for SU(N) and Nx2N matrix fields for Sp(2N), as specified by `Rev. 1.2 `_ of the ILDG Binary File Format.:: class IldgWriter : public ScidacWriter { - template - void writeConfiguration(Lattice > &Umu,int sequence,std::string LFN,std::string description); + template + void writeConfiguration(Lattice &Umu, int sequence, std::string LFN, std::string description) }; class IldgReader : public GridLimeReader { - template - void readConfiguration(Lattice > &Umu, FieldMetaData &FieldMetaData_) ; + + template + void readConfiguration(Lattice &Umu, FieldMetaData &FieldMetaData_); }; +``writeConfiguration``'s template parameters; ``group_name``, ``MatrixFormat``, and ``FloatingPointFormat`` specify the format the lattices are to be written to disk with. At present there are two options for the group, SU and Sp, two options for the matrix format, reduced and not reduced, and two options for the floating point format, 32 bit and 64 bit precision. Reduced format Sp gauge fields take up half as much storage as non-reduced fields while SU(N) fields take up N-1/N as much storage as their non-reduced counterparts. When reading reduced format lattices from disk there are multiple ways the full SU(N) field can be recovered. The ``unique_su`` template parameter of ``readConfiguration`` controls which of two alternative functions is used. Set ``unique_su=true`` if you want ``Grid`` to follow the prescription laid out in the ILDG Binary File Format specification when reconstructing reduced format lattices. + **Implementation detail** The lattice field routines internally use the above Binary routines to write the bulk data at an offset diff --git a/tests/IO/Test_IldgWriter.cc b/tests/IO/Test_IldgWriter.cc new file mode 100644 index 0000000000..003bb79717 --- /dev/null +++ b/tests/IO/Test_IldgWriter.cc @@ -0,0 +1,74 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./Test_IldgWriter.cc + + Copyright (C) 2026 + +Author: Azusa Yamaguchi +Author: Peter Boyle +Author: Gaurav Ray + + 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; + +// Unit tests to check the IldgWriter class. + +void checkWriteLimeIldgLFN(std::string &test_string); + +// write a lime record and then read it back and check it matches. +void checkWriteLimeIldgLFN(std::string &test_string) { + auto simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); + auto mpi_layout = GridDefaultMpi(); + Coordinate latt_size ({8,8,8,16}); + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + + std::string test_file = "test_file.lime"; + std::string check_string; + // open lime file + IldgWriter _ildgWriter(Grid.IsBoss()); + _ildgWriter.open(test_file); + _ildgWriter.writeLimeIldgLFN( test_string ); + _ildgWriter.close(); + // read back record + IldgReader _ildgReader; + _ildgReader.open(test_file); + _ildgReader.readLimeObject(check_string, ILDG_DATA_LFN); + std::cout < +Author: Peter Boyle +Author: Gaurav Ray + + 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; + +// this test demonstrates and checks IldgWriter/Readers' ability to +// read/write ildg 1.2 compliant SU and Sp lattices, including in reduced format. + +////////////////////////////////////////////// +// this template function returns a +// LatticeGaugeField of the chosen gauge +// group (passed as a template argument) +////////////////////////////////////////////// +template +LatticeGaugeField generateHotFieldConfiguration( GridCartesian &Grid, std::vector seed ) { + + GridParallelRNG pRNGa(&Grid); + std::cout < ) { + SU::HotConfiguration(pRNGa,Umu); + } else if constexpr ( std::is_same_v ) { + Sp::HotConfiguration(pRNGa,Umu); + } else { static_assert(true, "Grid does not recognise the gauge group"); } + + return Umu; +} + +///////////////////////////////////////////////////////////// +// this template function writes a lattice +// field of a given gaugeGroup to disk. It can write in +// reduced format and single/double precision depending on +// the values of matrix_fmt and fp_fmt. unique_su toggles +// between different functions when reading reduced fmt lattices. +///////////////////////////////////////////////////////////// +template +void writeReadIldgConfiguration( LatticeGaugeField &Umu, GridCartesian &Grid, std::string file) { + + if constexpr( std::is_same_v && N%2==1) { + std::cout <= 4" << std::endl; + std::cout <(Umu,4000,std::string("dummy_ildg_LFN"),std::string("dummy_config")); + _IldgWriter.close(); + + LatticeGaugeField Umu_saved(&Grid); + + FieldMetaData header; + + std::cout <(Umu_saved, header); + _IldgReader.close(); + std::cout < seed0 = {15,16,23,42}; + std::vector seed1 = {95,63,71,66}; + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + + const bool not_unique_su = false; // use reconstructSU + const bool unique_su = true; // use unique_reconstructSU + + LatticeGaugeField Umu = generateHotFieldConfiguration(Grid, seed0); + LatticeGaugeField UmuSp = generateHotFieldConfiguration(Grid, seed1); + + // write and read back SU lattices + writeReadIldgConfiguration(Umu, Grid, "./ckpoint_su64_"+std::to_string(Nc)+"x"+std::to_string(Nc)+".4000"); + writeReadIldgConfiguration(Umu, Grid, "./ckpoint_su64_"+std::to_string(Nc-1)+"x"+std::to_string(Nc)+".4000"); + writeReadIldgConfiguration(Umu, Grid, "./ckpoint_su64_unique_"+std::to_string(Nc-1)+"x"+std::to_string(Nc)+".4000"); + + writeReadIldgConfiguration(Umu, Grid, "./ckpoint_su32_"+std::to_string(Nc)+"x"+std::to_string(Nc)+".4000"); + writeReadIldgConfiguration(Umu, Grid, "./ckpoint_su32_"+std::to_string(Nc-1)+"x"+std::to_string(Nc)+".4000"); + writeReadIldgConfiguration(Umu, Grid, "./ckpoint_su32_unique_"+std::to_string(Nc-1)+"x"+std::to_string(Nc)+".4000"); + + // write and read Sp lattices + writeReadIldgConfiguration(UmuSp, Grid, "./ckpoint_sp64_"+std::to_string(Nc)+"x"+std::to_string(Nc)+".4000"); + writeReadIldgConfiguration(UmuSp, Grid, "./ckpoint_sp64_"+std::to_string(Nc/2)+"x"+std::to_string(Nc)+".4000"); + + writeReadIldgConfiguration(UmuSp, Grid, "./ckpoint_sp32_"+std::to_string(Nc)+"x"+std::to_string(Nc)+".4000"); + writeReadIldgConfiguration(UmuSp, Grid, "./ckpoint_sp32_"+std::to_string(Nc/2)+"x"+std::to_string(Nc)+".4000"); + + Grid_finalize(); +#endif +} diff --git a/tests/IO/Test_io_mungers.cc b/tests/IO/Test_io_mungers.cc new file mode 100644 index 0000000000..0e13946ded --- /dev/null +++ b/tests/IO/Test_io_mungers.cc @@ -0,0 +1,659 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./Test_io_mungers.cc + + Copyright (C) 2026 + +Author: Azusa Yamaguchi +Author: Peter Boyle +Author: Gaurav Ray + + 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; + +// tests to check the various (un)munger classes +// and functions in parallelIO/Metadata.h +void check_reconstruct3(); +void check_reconstructSU(); +void check_is_perm_even(); +void check_unique_reconstructSU(); +void checkGauge3x2mungers(); +void checkGaugeSUmungers(); +void check_reconstructSp(); +void checkGaugeSpmungers(); +void checkBinarySimpleMungers(); +void checkGaugeSimpleMungers(); +void checkGaugeDoubleStoredMungers(); + +void mock_SU_field(); +void mock_Sp_field(); + +////////////////////////////////////////////// +// helper functions for generating (mock) SU/Sp matrices +////////////////////////////////////////////// +void mock_SU_field( ColourMatrixD &cm, std::vector seed ) { + + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seed); + + ColourMatrixD lie, la, ta; + ComplexD ca; + + lie = Zero(); + + for (int a = 0; a < Nc*Nc-1; a++) { + random(sRNG, ca); + + ca = (ca + conjugate(ca)) * 0.5; + ca = ca - 0.5; + + SU::generator(a, ta); + + la = timesI(ca * ta); + + lie = lie + la; // e^{i la ta} + } + + cm = Exponentiate(lie, 2.0); + +} + +void mock_Sp_field( ColourMatrixD &cm, std::vector seed ) { + + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seed); + + std::vector elem(Nc*Nc/2); + for(auto&& e: elem) { + random(sRNG,e); + } + + // fill upper left and upper right blocks only + for(int i=0;i({30})); + mock_SU_field(SU3_yfield, std::vector({95})); + mock_SU_field(SU3_zfield, std::vector({7})); + mock_SU_field(SU3_tfield, std::vector({10})); + + //set last row equal to zero + for(int j=0;j(scalar, mu); + auto det = Determinant(new_cm); + assert( abs(norm2(det)-1.0) < 1e-12 ); + assert( norm2( new_cm*adj(new_cm)-1.0 ) < 1e-12 ); + } + +} + +void check_reconstructSU() { + + ColourMatrixD SUN_xfield, SUN_yfield, SUN_zfield, SUN_tfield; + + mock_SU_field(SUN_xfield, std::vector({4})); + mock_SU_field(SUN_yfield, std::vector({8})); + mock_SU_field(SUN_zfield, std::vector({15})); + mock_SU_field(SUN_tfield, std::vector({16})); + + //set last row equal to zero + for(int j=0;j(scalar, mu); + auto det = Determinant(new_cm); + assert( abs( norm2(det)-1.0 ) < 1e-12 ); + assert( norm2( new_cm*adj(new_cm)-1.0 ) < 1e-12 ); + } + +} + +/////////////////////////////////////////////////////////////////////////// +// test is_perm_even via the Johnson-Trotter algorithm which +// exhausts **ALL** the permutations of {0...Nc-1}. +// J-T generates permutations by swapping one adjacent pair of elements +// at a time. each step flips the parity - this is the basis of the test. +// implementation derived from: +// https://github.com/Lampjaw/Johnson-Trotter/blob/master/main.cpp +// ref: https://www.nist.gov/dads/HTML/johnsonTrotter.html +/////////////////////////////////////////////////////////////////////////// +void check_is_perm_even() { + bool isMobile = true; + int count = 1; + + // Initialization and establishment of number and direction lists + std::vector elements(Nc); + std::iota(elements.begin(), elements.end(), 0); + std::vector dirs(Nc, "<"); + + // start with an even permutation + assert( is_perm_even(elements) == true ); + + // Generate permutations as long as it's possible + while(isMobile) { + isMobile = false; + int mobileIndex = Nc; + elements[Nc] = 0; + + // Find the greatest mobile integer + for(int i = 0; i < Nc; i ++) { + // If integer is not on the edge of the list and pointing away from an integer + if(dirs[i] == "<" && i > 0){ + // If Integer is greater than the one it's pointing to and + if(elements[i] > elements[i-1] && elements[i] > elements[mobileIndex]) { + isMobile = true; // greater than the previously found mobile integer + mobileIndex = i; + } + } + else if(dirs[i] == ">" && i < Nc-1) { + if(elements[i] > elements[i+1] && elements[i] > elements[mobileIndex]) { + isMobile = true; + mobileIndex = i; + } + } + } + // If there is a valid mobile integer swap with the one it's pointing to + if(isMobile) { + if(dirs[mobileIndex] == "<") { // Integers pointing left + std::iter_swap(elements.begin()+mobileIndex-1, elements.begin()+mobileIndex); + std::iter_swap(dirs.begin()+mobileIndex-1, dirs.begin()+mobileIndex); + // Change the direction of every integer greater than the mobile integer + for(int i = 0; i < Nc; i++) { + if(elements[i] > elements[mobileIndex - 1]){ + if(dirs[i] == "<") { dirs[i] = ">"; } + else { dirs[i] = "<"; } + } + } + } + else { // Integers pointing right + std::iter_swap(elements.begin()+mobileIndex, elements.begin()+mobileIndex+1); + std::iter_swap(dirs.begin()+mobileIndex, dirs.begin()+mobileIndex+1); + // Change the direction of every integer greater than the mobile integer + for(int i = 0; i < Nc; i++) { + if(elements[i] > elements[mobileIndex + 1]){ + if(dirs[i] == "<") { dirs[i] = ">"; } + else { dirs[i] = "<"; } + } + } + } + // test if permutation is even or odd + if(count%2==0) { assert(is_perm_even(elements) == true); } + else { assert(is_perm_even(elements) == false); } + count++; + } + } +} + + +void check_unique_reconstructSU() { + + ColourMatrixD SUN_xfield, SUN_yfield, SUN_zfield, SUN_tfield; + + mock_SU_field(SUN_xfield, std::vector({2})); + mock_SU_field(SUN_yfield, std::vector({40})); + mock_SU_field(SUN_zfield, std::vector({28})); + mock_SU_field(SUN_tfield, std::vector({25})); + + //set last row equal to zero + for(int j=0;j(scalar, mu); + auto det = Determinant(new_cm); + assert( abs( norm2(det)-1.0 ) < 1e-12 ); + assert( norm2( new_cm*adj(new_cm)-1.0 ) < 1e-12 ); + } + +} + +void checkGauge3x2mungers() { + + ColourMatrixD SU3_xfield, SU3_yfield, SU3_zfield, SU3_tfield; + + mock_SU_field(SU3_xfield, std::vector({66})); + mock_SU_field(SU3_yfield, std::vector({63})); + mock_SU_field(SU3_zfield, std::vector({71})); + mock_SU_field(SU3_tfield, std::vector({98})); + + LorentzColourMatrixD scalar = Zero(); + + pokeLorentz(scalar, SU3_xfield, 0); + pokeLorentz(scalar, SU3_yfield, 1); + pokeLorentz(scalar, SU3_zfield, 2); + pokeLorentz(scalar, SU3_tfield, 3); + + // reduce field. Gauge3x2unmunger + Gauge3x2unmunger unmunger; + LorentzColour2x3D scalar_rdc; + unmunger(scalar,scalar_rdc); + + // reconstruct full field. Gauge3x2munger + Gauge3x2munger munger; + LorentzColourMatrixD scalar_recon; + munger(scalar_rdc,scalar_recon); + + // round-trip test + assert(norm2(scalar_recon-scalar)<1e-10); + +} + +template +void checkGaugeSUmungers() { + + ColourMatrixD SUN_xfield, SUN_yfield, SUN_zfield, SUN_tfield; + + mock_SU_field(SUN_xfield, std::vector({360})); + mock_SU_field(SUN_yfield, std::vector({804})); + mock_SU_field(SUN_zfield, std::vector({77})); + mock_SU_field(SUN_tfield, std::vector({65})); + + LorentzColourMatrixD scalar = Zero(); + + pokeLorentz(scalar, SUN_xfield, 0); + pokeLorentz(scalar, SUN_yfield, 1); + pokeLorentz(scalar, SUN_zfield, 2); + pokeLorentz(scalar, SUN_tfield, 3); + + // reduce field. GaugeSUunmunger + GaugeSUunmunger unmunger; + LorentzColour2x3D scalar_rdc = Zero(); + unmunger(scalar,scalar_rdc); + + // reconstruct full field. GaugeSUmunger + GaugeSUmunger munger; + LorentzColourMatrixD scalar_recon; + munger(scalar_rdc, scalar_recon); + + // round-trip test + assert(norm2(scalar_recon-scalar)<1e-12); + +} + + +void check_reconstructSp() { + + ColourMatrixD Sp_xfield, Sp_yfield, Sp_zfield, Sp_tfield; + Sp_xfield = Sp_yfield = Sp_zfield = Sp_tfield = Zero(); + + mock_Sp_field(Sp_xfield, std::vector({1066})); + mock_Sp_field(Sp_yfield, std::vector({1995})); + mock_Sp_field(Sp_zfield, std::vector({997})); + mock_Sp_field(Sp_tfield, std::vector({58})); + + LorentzColourMatrixD scalar = Zero(); + + pokeLorentz(scalar, Sp_xfield, 0); + pokeLorentz(scalar, Sp_yfield, 1); + pokeLorentz(scalar, Sp_zfield, 2); + pokeLorentz(scalar, Sp_tfield, 3); + + reconstructSp(scalar); + + // test Sp block structure + // A B + // -B* A* + for(int mu=0; mu(scalar,mu); + for(int i=0;i({23})); + mock_Sp_field(Sp_yfield, std::vector({42})); + mock_Sp_field(Sp_zfield, std::vector({93})); + mock_Sp_field(Sp_tfield, std::vector({49})); + + LorentzColourMatrixD scalar = Zero(); + + pokeLorentz(scalar, Sp_xfield, 0); + pokeLorentz(scalar, Sp_yfield, 1); + pokeLorentz(scalar, Sp_zfield, 2); + pokeLorentz(scalar, Sp_tfield, 3); + + reconstructSp(scalar); + + // reduce field. GaugeSpunmunger + GaugeSpunmunger unmunger; + LorentzColourNx2ND scalar_rdc = Zero(); + unmunger(scalar,scalar_rdc); + + // reconstruct full field. GaugeSpmunger + GaugeSpmunger munger; + LorentzColourMatrixD scalar_recon; + munger(scalar_rdc,scalar_recon); + + // round-trip test + assert(scalar==scalar_recon); + +} + + +void checkBinarySimpleMungers() { + + LorentzColourMatrixF in_scalar_objectF; // single precision + LorentzColourMatrixF out_scalar_objectF = Zero(); + LorentzColourMatrixD in_scalar_objectD; // double precision + LorentzColourMatrixD out_scalar_objectD = Zero(); + + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(std::vector({117})); + + random(sRNG,in_scalar_objectF); + random(sRNG,in_scalar_objectD); + + // BinarySimpleUnmunger + BinarySimpleUnmunger unmungerFF; + BinarySimpleUnmunger unmungerFD; + BinarySimpleUnmunger unmungerDD; + BinarySimpleUnmunger unmungerDF; + + unmungerFF(in_scalar_objectF, out_scalar_objectF); + unmungerDD(in_scalar_objectD, out_scalar_objectD); + assert(in_scalar_objectF==out_scalar_objectF); + assert(in_scalar_objectD==out_scalar_objectD); + + random(sRNG,in_scalar_objectF); + unmungerFD(in_scalar_objectF, out_scalar_objectD); + unmungerDF(out_scalar_objectD, out_scalar_objectF); + assert(in_scalar_objectF==out_scalar_objectF); + + random(sRNG,in_scalar_objectD); + unmungerDF(in_scalar_objectD, out_scalar_objectF); + unmungerFD(out_scalar_objectF, out_scalar_objectD); + // lose exactness when going from double-->single precision + assert(norm2( in_scalar_objectD - out_scalar_objectD ) < 1e-10); + + // BinarySimpleMunger + BinarySimpleMunger mungerFF; + BinarySimpleMunger mungerDF; + BinarySimpleMunger mungerDD; + BinarySimpleMunger mungerFD; + + random(sRNG,in_scalar_objectF); + random(sRNG,in_scalar_objectD); + + mungerFF(in_scalar_objectF, out_scalar_objectF); + mungerDD(in_scalar_objectD, out_scalar_objectD); + assert(in_scalar_objectD==out_scalar_objectD); + assert(in_scalar_objectF==out_scalar_objectF); + + random(sRNG,in_scalar_objectF); + mungerFD(in_scalar_objectF, out_scalar_objectD); + mungerDF(out_scalar_objectD, out_scalar_objectF); + assert(in_scalar_objectF==out_scalar_objectF); + + random(sRNG,in_scalar_objectD); + mungerDF(in_scalar_objectD, out_scalar_objectF); + mungerFD(out_scalar_objectF, out_scalar_objectD); + // lose exactness when going from double-->single precision + assert(norm2( in_scalar_objectD - out_scalar_objectD ) < 1e-10); + +} + +// these are used in NerscIO.h +void checkGaugeSimpleMungers() { + + LorentzColourMatrixF in_scalar_objectF; // single precision + LorentzColourMatrixF out_scalar_objectF = Zero(); + LorentzColourMatrixD in_scalar_objectD; // double precision + LorentzColourMatrixD out_scalar_objectD = Zero(); + + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(std::vector({57})); + + random(sRNG,in_scalar_objectF); + random(sRNG,in_scalar_objectD); + + // GaugeSimpleUnmunger + GaugeSimpleUnmunger unmungerFF; + GaugeSimpleUnmunger unmungerFD; + GaugeSimpleUnmunger unmungerDD; + GaugeSimpleUnmunger unmungerDF; + + unmungerFF(in_scalar_objectF, out_scalar_objectF); + unmungerDD(in_scalar_objectD, out_scalar_objectD); + assert(in_scalar_objectF==out_scalar_objectF); + assert(in_scalar_objectD==out_scalar_objectD); + + random(sRNG,in_scalar_objectF); + unmungerFD(in_scalar_objectF, out_scalar_objectD); + unmungerDF(out_scalar_objectD, out_scalar_objectF); + assert(in_scalar_objectF==out_scalar_objectF); + + random(sRNG,in_scalar_objectD); + unmungerDF(in_scalar_objectD, out_scalar_objectF); + unmungerFD(out_scalar_objectF, out_scalar_objectD); + // lose exactness when going from double-->single precision + assert(norm2( in_scalar_objectD - out_scalar_objectD ) < 1e-10); + + // GaugeSimpleMunger + GaugeSimpleMunger mungerFF; + GaugeSimpleMunger mungerDF; + GaugeSimpleMunger mungerDD; + GaugeSimpleMunger mungerFD; + + random(sRNG,in_scalar_objectF); + random(sRNG,in_scalar_objectD); + + mungerFF(in_scalar_objectF, out_scalar_objectF); + mungerDD(in_scalar_objectD, out_scalar_objectD); + assert(in_scalar_objectD==out_scalar_objectD); + assert(in_scalar_objectF==out_scalar_objectF); + + random(sRNG,in_scalar_objectF); + mungerFD(in_scalar_objectF, out_scalar_objectD); + mungerDF(out_scalar_objectD, out_scalar_objectF); + assert(in_scalar_objectF==out_scalar_objectF); + + random(sRNG,in_scalar_objectD); + mungerDF(in_scalar_objectD, out_scalar_objectF); + mungerFD(out_scalar_objectF, out_scalar_objectD); + // lose exactness when going from double-->single precision + assert(norm2( in_scalar_objectD - out_scalar_objectD ) < 1e-10); + +} + +void checkGaugeDoubleStoredMungers() { + + DoubleStoredColourMatrixF in_scalar_objectF; // single precision + DoubleStoredColourMatrixF out_scalar_objectF = Zero(); + DoubleStoredColourMatrixD in_scalar_objectD; // double precision + DoubleStoredColourMatrixD out_scalar_objectD = Zero(); + + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(std::vector({169})); + + random(sRNG,in_scalar_objectF); + random(sRNG,in_scalar_objectD); + + // GaugeDoubleStoredUnmunger + GaugeDoubleStoredUnmunger unmungerFF; + GaugeDoubleStoredUnmunger unmungerFD; + GaugeDoubleStoredUnmunger unmungerDD; + GaugeDoubleStoredUnmunger unmungerDF; + + unmungerFF(in_scalar_objectF, out_scalar_objectF); + unmungerDD(in_scalar_objectD, out_scalar_objectD); + assert(in_scalar_objectF==out_scalar_objectF); + assert(in_scalar_objectD==out_scalar_objectD); + + random(sRNG,in_scalar_objectF); + unmungerFD(in_scalar_objectF, out_scalar_objectD); + unmungerDF(out_scalar_objectD, out_scalar_objectF); + assert(in_scalar_objectF==out_scalar_objectF); + + random(sRNG,in_scalar_objectD); + unmungerDF(in_scalar_objectD, out_scalar_objectF); + unmungerFD(out_scalar_objectF, out_scalar_objectD); + // lose exactness when going from double-->single precision + assert(norm2( in_scalar_objectD - out_scalar_objectD ) < 1e-10); + + // GaugeDoubleStoredMunger + GaugeDoubleStoredMunger mungerFF; + GaugeDoubleStoredMunger mungerDF; + GaugeDoubleStoredMunger mungerDD; + GaugeDoubleStoredMunger mungerFD; + + random(sRNG,in_scalar_objectF); + random(sRNG,in_scalar_objectD); + + mungerFF(in_scalar_objectF, out_scalar_objectF); + mungerDD(in_scalar_objectD, out_scalar_objectD); + assert(in_scalar_objectD==out_scalar_objectD); + assert(in_scalar_objectF==out_scalar_objectF); + + random(sRNG,in_scalar_objectF); + mungerFD(in_scalar_objectF, out_scalar_objectD); + mungerDF(out_scalar_objectD, out_scalar_objectF); + assert(in_scalar_objectF==out_scalar_objectF); + + random(sRNG,in_scalar_objectD); + mungerDF(in_scalar_objectD, out_scalar_objectF); + mungerFD(out_scalar_objectF, out_scalar_objectD); + // lose exactness when going from double-->single precision + assert(norm2( in_scalar_objectD - out_scalar_objectD ) < 1e-10); + +} + + +int main (int argc, char ** argv) +{ +#ifdef HAVE_LIME + Grid_init(&argc,&argv); + + std::cout <1 && Nc<4) { + check_reconstruct3(); + std::cout << GridLogMessage << "reconstruct3: PASS" << std::endl; + + checkGauge3x2mungers(); + std::cout << GridLogMessage << "Gauge3x2mungers: PASS" << std::endl; + } + + check_reconstructSU(); + std::cout << GridLogMessage << "reconstructSU: PASS" << std::endl; + + checkGaugeSUmungers(); + std::cout << GridLogMessage << "(unique_su=false) GaugeSUmungers: PASS" << std::endl; + + check_unique_reconstructSU(); + std::cout << GridLogMessage << "unique_reconstructSU: PASS" << std::endl; + + checkGaugeSUmungers(); + std::cout << GridLogMessage << "(unique_su=true) GaugeSUmungers: PASS" << std::endl; + + if constexpr(Nc>2 && Nc%2==0) { + std::cout < +Author: neo +Author: Guido Cossu +Author: Gaurav Ray + +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 + +int main(int argc, char **argv) +{ + using namespace Grid; + + Grid_init(&argc, &argv); + GridLogLayout(); + + typedef GenericSpHMCRunner HMCWrapper; + HMCWrapper TheHMC; + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint/lat"; + CPparams.rng_prefix = "ckpoint/rng"; + CPparams.smeared_prefix = "ckpoint/smr"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE32BIG"; + CPparams.saveSmeared = true; + + // Ildg specific parameters + CPparams.group = "sp"; + CPparams.reduced_matrix = true; + //CPparams.unique_su = true; // only applies to su fields + + TheHMC.Resources.LoadILDGCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "55 1 37 3 97 8"; + RNGpar.parallel_seeds = "4 8 15 16 23 42"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef TopologicalChargeMod QObs; + TopologyObsParameters TopParams; + TopParams.interval = 5; + TopParams.do_smearing = true; + TopParams.Smearing.init_step_size = 0.01; + TopParams.Smearing.tolerance = 1e-5; + TopParams.Smearing.meas_interval = 50; + TopParams.Smearing.maxTau = 2.0; + TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + RealD beta = 8.0 ; + SpWilsonGaugeActionR Waction(beta); + + ActionLevel Level1(1); + Level1.push_back(&Waction); + TheHMC.TheAction.push_back(Level1); + ///////////////////////////////////////////////////////////// + + // HMC parameters are serialisable + TheHMC.Parameters.MD.MDsteps = 10; + TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + TheHMC.Run(); + + Grid_finalize(); + +} // main