diff --git a/Grid/Makefile.am b/Grid/Makefile.am index 7c3c151b..aff5b9d0 100644 --- a/Grid/Makefile.am +++ b/Grid/Makefile.am @@ -70,7 +70,7 @@ endif lib_LIBRARIES = libGrid.a CCFILES += $(extra_sources) -HFILES += $(extra_headers) Config.h Version.h +HFILES += $(extra_headers) libGrid_a_SOURCES = $(CCFILES) libGrid_adir = $(includedir)/Grid diff --git a/Grid/parallelIO/IldgIO.h b/Grid/parallelIO/IldgIO.h index 576f08bd..80d135d2 100644 --- a/Grid/parallelIO/IldgIO.h +++ b/Grid/parallelIO/IldgIO.h @@ -31,6 +31,7 @@ directory #include #include #include +#include #include #include @@ -654,7 +655,8 @@ class IldgWriter : public ScidacWriter { // Fill ILDG header data struct ////////////////////////////////////////////////////// ildgFormat ildgfmt ; - ildgfmt.field = std::string("su3gauge"); + const std::string stNC = std::to_string( Nc ) ; + ildgfmt.field = std::string("su"+stNC+"gauge"); if ( format == std::string("IEEE32BIG") ) { ildgfmt.precision = 32; @@ -871,7 +873,8 @@ class IldgReader : public GridLimeReader { } else { assert(found_ildgFormat); - assert ( ildgFormat_.field == std::string("su3gauge") ); + const std::string stNC = std::to_string( Nc ) ; + assert ( ildgFormat_.field == std::string("su"+stNC+"gauge") ); /////////////////////////////////////////////////////////////////////////////////////// // Populate our Grid metadata as best we can @@ -879,7 +882,7 @@ class IldgReader : public GridLimeReader { std::ostringstream vers; vers << ildgFormat_.version; FieldMetaData_.hdr_version = vers.str(); - FieldMetaData_.data_type = std::string("4D_SU3_GAUGE_3X3"); + FieldMetaData_.data_type = std::string("4D_SU"+stNC+"_GAUGE_"+stNC+"x"+stNC); FieldMetaData_.nd=4; FieldMetaData_.dimension.resize(4); diff --git a/Grid/parallelIO/MetaData.h b/Grid/parallelIO/MetaData.h index af8b3f76..6b9d8708 100644 --- a/Grid/parallelIO/MetaData.h +++ b/Grid/parallelIO/MetaData.h @@ -6,8 +6,8 @@ Copyright (C) 2015 - Author: Peter Boyle + Author: Jamie Hudspith 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 @@ -182,8 +182,8 @@ class GaugeStatistics public: void operator()(Lattice & data,FieldMetaData &header) { - header.link_trace=WilsonLoops::linkTrace(data); - header.plaquette =WilsonLoops::avgPlaquette(data); + header.link_trace = WilsonLoops::linkTrace(data); + header.plaquette = WilsonLoops::avgPlaquette(data); } }; typedef GaugeStatistics PeriodicGaugeStatistics; @@ -203,20 +203,24 @@ template<> inline void PrepareMetaData(Lattice 1 ) ; for(int mu=0;mu using iLorentzColour2x3 = iVector, 2>, Nd >; +template using iLorentzColour2x3 = iVector, Nc-1>, Nd >; typedef iLorentzColour2x3 LorentzColour2x3; typedef iLorentzColour2x3 LorentzColour2x3F; @@ -278,7 +282,6 @@ struct GaugeSimpleMunger{ template struct GaugeSimpleUnmunger { - void operator()(sobj &in, fobj &out) { for (int mu = 0; mu < Nd; mu++) { for (int i = 0; i < Nc; i++) { @@ -317,8 +320,8 @@ template struct Gauge3x2munger{ void operator() (fobj &in,sobj &out){ for(int mu=0;mu struct Gauge3x2unmunger{ void operator() (sobj &in,fobj &out){ for(int mu=0;mu Author: Peter Boyle Author: paboyle + Author: Jamie Hudspith 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 @@ -30,6 +31,8 @@ #ifndef GRID_NERSC_IO_H #define GRID_NERSC_IO_H +#include + NAMESPACE_BEGIN(Grid); using namespace Grid; @@ -145,15 +148,17 @@ public: std::string format(header.floating_point); - int ieee32big = (format == std::string("IEEE32BIG")); - int ieee32 = (format == std::string("IEEE32")); - int ieee64big = (format == std::string("IEEE64BIG")); - int ieee64 = (format == std::string("IEEE64") || format == std::string("IEEE64LITTLE")); + const int ieee32big = (format == std::string("IEEE32BIG")); + const int ieee32 = (format == std::string("IEEE32")); + const int ieee64big = (format == std::string("IEEE64BIG")); + const int ieee64 = (format == std::string("IEEE64") || \ + format == std::string("IEEE64LITTLE")); uint32_t nersc_csum,scidac_csuma,scidac_csumb; // depending on datatype, set up munger; // munger is a function of - if ( header.data_type == std::string("4D_SU3_GAUGE") ) { + const std::string stNC = std::to_string( Nc ) ; + if ( header.data_type == std::string("4D_SU"+stNC+"_GAUGE") ) { if ( ieee32 || ieee32big ) { BinaryIO::readLatticeObject (Umu,file,Gauge3x2munger(), offset,format, @@ -164,7 +169,7 @@ public: (Umu,file,Gauge3x2munger(),offset,format, nersc_csum,scidac_csuma,scidac_csumb); } - } else if ( header.data_type == std::string("4D_SU3_GAUGE_3x3") ) { + } else if ( header.data_type == std::string("4D_SU"+stNC+"_GAUGE_"+stNC+"x"+stNC) ) { if ( ieee32 || ieee32big ) { BinaryIO::readLatticeObject (Umu,file,GaugeSimpleMunger(),offset,format, @@ -231,6 +236,7 @@ public: header.sequence_number = sequence_number; header.ensemble_id = ens_id; header.ensemble_label = ens_label; + header.hdr_version = "1.0" ; typedef LorentzColourMatrixD fobj3D; typedef LorentzColour2x3D fobj2D; @@ -244,10 +250,14 @@ public: uint64_t offset; - // Sod it -- always write 3x3 double - header.floating_point = std::string("IEEE64BIG"); - header.data_type = std::string("4D_SU3_GAUGE_3x3"); - GaugeSimpleUnmunger munge; + // Sod it -- always write NcxNc double + header.floating_point = std::string("IEEE64BIG"); + const std::string stNC = std::to_string( Nc ) ; + if( two_row ) { + header.data_type = std::string("4D_SU" + stNC + "_GAUGE" ); + } else { + header.data_type = std::string("4D_SU" + stNC + "_GAUGE_" + stNC + "x" + stNC ); + } if ( grid->IsBoss() ) { truncate(file); offset = writeHeader(header,file); @@ -255,8 +265,15 @@ public: grid->Broadcast(0,(void *)&offset,sizeof(offset)); uint32_t nersc_csum,scidac_csuma,scidac_csumb; - BinaryIO::writeLatticeObject(Umu,file,munge,offset,header.floating_point, - nersc_csum,scidac_csuma,scidac_csumb); + if( two_row ) { + Gauge3x2unmunger munge; + BinaryIO::writeLatticeObject(Umu,file,munge,offset,header.floating_point, + nersc_csum,scidac_csuma,scidac_csumb); + } else { + GaugeSimpleUnmunger munge; + BinaryIO::writeLatticeObject(Umu,file,munge,offset,header.floating_point, + nersc_csum,scidac_csuma,scidac_csumb); + } header.checksum = nersc_csum; if ( grid->IsBoss() ) { writeHeader(header,file); @@ -288,8 +305,7 @@ public: header.plaquette=0.0; MachineCharacteristics(header); - uint64_t offset; - + uint64_t offset; #ifdef RNG_RANLUX header.floating_point = std::string("UINT64"); header.data_type = std::string("RANLUX48"); @@ -329,7 +345,7 @@ public: GridBase *grid = parallel.Grid(); - uint64_t offset = readHeader(file,grid,header); + uint64_t offset = readHeader(file,grid,header); FieldMetaData clone(header); diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h index 0367c9fa..e002e3d5 100644 --- a/Grid/qcd/utils/WilsonLoops.h +++ b/Grid/qcd/utils/WilsonLoops.h @@ -125,7 +125,6 @@ public: return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME } - ////////////////////////////////////////////////// // average over all x,y,z the temporal loop ////////////////////////////////////////////////// @@ -165,7 +164,7 @@ public: double vol = Umu.Grid()->gSites(); - return p.real() / vol / 4.0 / 3.0; + return p.real() / vol / (4.0 * Nc ) ; }; ////////////////////////////////////////////////// diff --git a/configure.ac b/configure.ac index 406b0b74..9ab0595a 100644 --- a/configure.ac +++ b/configure.ac @@ -159,7 +159,7 @@ case ${ac_ZMOBIUS} in esac ############### Nc AC_ARG_ENABLE([Nc], - [AC_HELP_STRING([--enable-Nc=2|3|4], [enable number of colours])], + [AC_HELP_STRING([--enable-Nc=2|3|4|5], [enable number of colours])], [ac_Nc=${enable_Nc}], [ac_Nc=3]) case ${ac_Nc} in diff --git a/tests/IO/Test_nersc_io.cc b/tests/IO/Test_nersc_io.cc index c15c320e..ae5b9c0d 100644 --- a/tests/IO/Test_nersc_io.cc +++ b/tests/IO/Test_nersc_io.cc @@ -147,7 +147,7 @@ int main (int argc, char ** argv) Complex p = TensorRemove(Tp); std::cout< Author: Peter Boyle Author: Guido Cossu +Author: Jamie Hudspith 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,14 +43,14 @@ directory using namespace std; using namespace Grid; - ; +; int main(int argc, char** argv) { Grid_init(&argc, &argv); std::vector latt({4, 4, 4, 8}); GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( - latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); + latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid); @@ -60,15 +61,19 @@ int main(int argc, char** argv) { << std::endl; SU2::printGenerators(); std::cout << "Dimension of adjoint representation: "<< SU2Adjoint::Dimension << std::endl; + + // guard as this code fails to compile for Nc != 3 +#if (Nc == 3) + SU2Adjoint::printGenerators(); SU2::testGenerators(); SU2Adjoint::testGenerators(); - + std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "* Generators for SU(Nc" << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; SU3::printGenerators(); std::cout << "Dimension of adjoint representation: "<< SU3Adjoint::Dimension << std::endl; SU3Adjoint::printGenerators(); @@ -111,12 +116,10 @@ int main(int argc, char** argv) { // AdjointRepresentation has the predefined number of colours Nc // Representations RepresentationTypes(grid); - - LatticeGaugeField U(grid), V(grid); SU3::HotConfiguration(gridRNG, U); SU3::HotConfiguration(gridRNG, V); - + // Adjoint representation // Test group structure // (U_f * V_f)_r = U_r * V_r @@ -127,17 +130,17 @@ int main(int argc, char** argv) { SU3::LatticeMatrix Vmu = peekLorentz(V,mu); pokeLorentz(UV,Umu*Vmu, mu); } - + AdjRep.update_representation(UV); typename AdjointRep::LatticeField UVr = AdjRep.U; // (U_f * V_f)_r - - + + AdjRep.update_representation(U); typename AdjointRep::LatticeField Ur = AdjRep.U; // U_r - + AdjRep.update_representation(V); typename AdjointRep::LatticeField Vr = AdjRep.U; // V_r - + typename AdjointRep::LatticeField UrVr(grid); UrVr = Zero(); for (int mu = 0; mu < Nd; mu++) { @@ -145,10 +148,10 @@ int main(int argc, char** argv) { typename AdjointRep::LatticeMatrix Vrmu = peekLorentz(Vr,mu); pokeLorentz(UrVr,Urmu*Vrmu, mu); } - + typename AdjointRep::LatticeField Diff_check = UVr - UrVr; std::cout << GridLogMessage << "Group structure SU("<::AdjointLieAlgebraMatrix(h_adj,Ar); - + // Re-extract h_adj SU3::LatticeAlgebraVector h_adj2(grid); SU_Adjoint::projectOnAlgebra(h_adj2, Ar); SU3::LatticeAlgebraVector h_diff = h_adj - h_adj2; std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl; - + // Exponentiate typename AdjointRep::LatticeMatrix Uadj(grid); Uadj = expMat(Ar, 1.0, 16); - + typename AdjointRep::LatticeMatrix uno(grid); uno = 1.0; // Check matrix Uadj, must be real orthogonal typename AdjointRep::LatticeMatrix Ucheck = Uadj - conjugate(Uadj); std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck) - << std::endl; - + << std::endl; + Ucheck = Uadj * adj(Uadj) - uno; std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck) - << std::endl; + << std::endl; Ucheck = adj(Uadj) * Uadj - uno; std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck) - << std::endl; - - + << std::endl; + // Construct the fundamental matrix in the group SU3::LatticeMatrix Af(grid); SU3::FundamentalLieAlgebraMatrix(h_adj,Af); @@ -193,72 +195,65 @@ int main(int argc, char** argv) { SU3::LatticeMatrix UnitCheck(grid); UnitCheck = Ufund * adj(Ufund) - uno_f; std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck) - << std::endl; + << std::endl; UnitCheck = adj(Ufund) * Ufund - uno_f; std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck) - << std::endl; - + << std::endl; // Tranform to the adjoint representation U = Zero(); // fill this with only one direction pokeLorentz(U,Ufund,0); // the representation transf acts on full gauge fields - + AdjRep.update_representation(U); Ur = AdjRep.U; // U_r typename AdjointRep::LatticeMatrix Ur0 = peekLorentz(Ur,0); // this should be the same as Uadj - + typename AdjointRep::LatticeMatrix Diff_check_mat = Ur0 - Uadj; std::cout << GridLogMessage << "Projections structure check group difference : " << norm2(Diff_check_mat) << std::endl; - - - // TwoIndexRep tests - std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "* eS^{ij} base for SU(2)" << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "Dimension of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl; SU2TwoIndexSymm::printBase(); - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "Generators of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl; + std::cout << GridLogMessage << "*********************************************" + << std::endl; + std::cout << GridLogMessage << "Generators of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl; SU2TwoIndexSymm::printGenerators(); - std::cout << GridLogMessage << "Test of Two Index Symmetric Generators: "<< SU2TwoIndexSymm::Dimension << std::endl; + std::cout << GridLogMessage << "Test of Two Index Symmetric Generators: "<< SU2TwoIndexSymm::Dimension << std::endl; SU2TwoIndexSymm::testGenerators(); std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "* eAS^{ij} base for SU(2)" << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl; SU2TwoIndexAntiSymm::printBase(); - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl; + std::cout << GridLogMessage << "*********************************************" + << std::endl; + std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl; SU2TwoIndexAntiSymm::printGenerators(); std::cout << GridLogMessage << "Test of Two Index anti-Symmetric Generators: "<< SU2TwoIndexAntiSymm::Dimension << std::endl; SU2TwoIndexAntiSymm::testGenerators(); - - std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "Test for the Two Index Symmetric projectors" - << std::endl; + << std::endl; // Projectors SU3TwoIndexSymm::LatticeTwoIndexMatrix Gauss2(grid); random(gridRNG,Gauss2); @@ -276,13 +271,13 @@ int main(int argc, char** argv) { SU3::LatticeAlgebraVector diff2 = ha - hb; std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; - std::cout << GridLogMessage << "*********************************************" - << std::endl; + std::cout << GridLogMessage << "*********************************************" + << std::endl; std::cout << GridLogMessage << "Test for the Two index anti-Symmetric projectors" - << std::endl; + << std::endl; // Projectors SU3TwoIndexAntiSymm::LatticeTwoIndexMatrix Gauss2a(grid); random(gridRNG,Gauss2a); @@ -300,11 +295,11 @@ int main(int argc, char** argv) { SU3::LatticeAlgebraVector diff2a = ha - hb; std::cout << GridLogMessage << "Difference: " << norm2(diff2a) << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "Two index Symmetric: Checking Group Structure" - << std::endl; + << std::endl; // Testing HMC representation classes TwoIndexRep< Nc, Symmetric > TIndexRep(grid); @@ -313,7 +308,7 @@ int main(int argc, char** argv) { LatticeGaugeField U2(grid), V2(grid); SU3::HotConfiguration(gridRNG, U2); SU3::HotConfiguration(gridRNG, V2); - + LatticeGaugeField UV2(grid); UV2 = Zero(); for (int mu = 0; mu < Nd; mu++) { @@ -321,16 +316,16 @@ int main(int argc, char** argv) { SU3::LatticeMatrix Vmu2 = peekLorentz(V2,mu); pokeLorentz(UV2,Umu2*Vmu2, mu); } - + TIndexRep.update_representation(UV2); typename TwoIndexRep< Nc, Symmetric >::LatticeField UVr2 = TIndexRep.U; // (U_f * V_f)_r - + TIndexRep.update_representation(U2); typename TwoIndexRep< Nc, Symmetric >::LatticeField Ur2 = TIndexRep.U; // U_r - + TIndexRep.update_representation(V2); typename TwoIndexRep< Nc, Symmetric >::LatticeField Vr2 = TIndexRep.U; // V_r - + typename TwoIndexRep< Nc, Symmetric >::LatticeField Ur2Vr2(grid); Ur2Vr2 = Zero(); for (int mu = 0; mu < Nd; mu++) { @@ -338,11 +333,11 @@ int main(int argc, char** argv) { typename TwoIndexRep< Nc, Symmetric >::LatticeMatrix Vrmu2 = peekLorentz(Vr2,mu); pokeLorentz(Ur2Vr2,Urmu2*Vrmu2, mu); } - + typename TwoIndexRep< Nc, Symmetric >::LatticeField Diff_check2 = UVr2 - Ur2Vr2; std::cout << GridLogMessage << "Group structure SU("<::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym); - + // Re-extract h_sym SU3::LatticeAlgebraVector h_sym2(grid); SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym); SU3::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2; std::cout << GridLogMessage << "Projections structure check vector difference (Two Index Symmetric): " << norm2(h_diff_sym) << std::endl; - - + // Exponentiate typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix U2iS(grid); U2iS = expMat(Ar_sym, 1.0, 16); - + typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix uno2iS(grid); uno2iS = 1.0; // Check matrix U2iS, must be real orthogonal typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ucheck2iS = U2iS - conjugate(U2iS); std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iS) - << std::endl; - + << std::endl; + Ucheck2iS = U2iS * adj(U2iS) - uno2iS; std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iS) - << std::endl; + << std::endl; Ucheck2iS = adj(U2iS) * U2iS - uno2iS; std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iS) - << std::endl; - - - + << std::endl; + // Construct the fundamental matrix in the group SU3::LatticeMatrix Af_sym(grid); SU3::FundamentalLieAlgebraMatrix(h_sym,Af_sym); @@ -386,147 +378,137 @@ int main(int argc, char** argv) { SU3::LatticeMatrix UnitCheck2(grid); UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f; std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2) - << std::endl; + << std::endl; UnitCheck2 = adj(Ufund2) * Ufund2 - uno_f; std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2) - << std::endl; - - + << std::endl; + + // Tranform to the 2Index Sym representation U = Zero(); // fill this with only one direction pokeLorentz(U,Ufund2,0); // the representation transf acts on full gauge fields - + TIndexRep.update_representation(U); Ur2 = TIndexRep.U; // U_r typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ur02 = peekLorentz(Ur2,0); // this should be the same as U2iS - + typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Diff_check_mat2 = Ur02 - U2iS; std::cout << GridLogMessage << "Projections structure check group difference (Two Index Symmetric): " << norm2(Diff_check_mat2) << std::endl; - - - - + if (TwoIndexRep::Dimension != 1){ - std::cout << GridLogMessage << "*********************************************" - << std::endl; + std::cout << GridLogMessage << "*********************************************" + << std::endl; - std::cout << GridLogMessage << "Two Index anti-Symmetric: Check Group Structure" - << std::endl; - // Testing HMC representation classes - TwoIndexRep< Nc, AntiSymmetric > TIndexRepA(grid); + std::cout << GridLogMessage << "Two Index anti-Symmetric: Check Group Structure" + << std::endl; + // Testing HMC representation classes + TwoIndexRep< Nc, AntiSymmetric > TIndexRepA(grid); - // Test group structure - // (U_f * V_f)_r = U_r * V_r - LatticeGaugeField U2A(grid), V2A(grid); - SU3::HotConfiguration(gridRNG, U2A); - SU3::HotConfiguration(gridRNG, V2A); + // Test group structure + // (U_f * V_f)_r = U_r * V_r + LatticeGaugeField U2A(grid), V2A(grid); + SU3::HotConfiguration(gridRNG, U2A); + SU3::HotConfiguration(gridRNG, V2A); - LatticeGaugeField UV2A(grid); - UV2A = Zero(); - for (int mu = 0; mu < Nd; mu++) { - SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu); - SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu); - pokeLorentz(UV2A,Umu2A*Vmu2A, mu); + LatticeGaugeField UV2A(grid); + UV2A = Zero(); + for (int mu = 0; mu < Nd; mu++) { + SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu); + SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu); + pokeLorentz(UV2A,Umu2A*Vmu2A, mu); + } + + TIndexRep.update_representation(UV2A); + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField UVr2A = TIndexRepA.U; // (U_f * V_f)_r + + TIndexRep.update_representation(U2A); + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2A = TIndexRepA.U; // U_r + + TIndexRep.update_representation(V2A); + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Vr2A = TIndexRepA.U; // V_r + + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2Vr2A(grid); + Ur2Vr2A = Zero(); + for (int mu = 0; mu < Nd; mu++) { + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Urmu2A = peekLorentz(Ur2A,mu); + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Vrmu2A = peekLorentz(Vr2A,mu); + pokeLorentz(Ur2Vr2A,Urmu2A*Vrmu2A, mu); + } + + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Diff_check2A = UVr2A - Ur2Vr2A; + std::cout << GridLogMessage << "Group structure SU("<::LatticeMatrix Ar_Asym(grid); + random(gridRNG,h_Asym); + h_Asym = real(h_Asym); + SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym); + + // Re-extract h_sym + SU3::LatticeAlgebraVector h_Asym2(grid); + SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym); + SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2; + std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl; + + + // Exponentiate + typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix U2iAS(grid); + U2iAS = expMat(Ar_Asym, 1.0, 16); + + typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix uno2iAS(grid); + uno2iAS = 1.0; + // Check matrix U2iS, must be real orthogonal + typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ucheck2iAS = U2iAS - conjugate(U2iAS); + std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iAS) + << std::endl; + + Ucheck2iAS = U2iAS * adj(U2iAS) - uno2iAS; + std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iAS) + << std::endl; + Ucheck2iAS = adj(U2iAS) * U2iAS - uno2iAS; + std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iAS) + << std::endl; + + + + // Construct the fundamental matrix in the group + SU3::LatticeMatrix Af_Asym(grid); + SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym); + SU3::LatticeMatrix Ufund2A(grid); + Ufund2A = expMat(Af_Asym, 1.0, 16); + SU3::LatticeMatrix UnitCheck2A(grid); + UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f; + std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A) + << std::endl; + UnitCheck2A = adj(Ufund2A) * Ufund2A - uno_f; + std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2A) + << std::endl; + + + // Tranform to the 2Index Sym representation + U = Zero(); // fill this with only one direction + pokeLorentz(U,Ufund2A,0); // the representation transf acts on full gauge fields + + TIndexRepA.update_representation(U); + Ur2A = TIndexRepA.U; // U_r + typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ur02A = peekLorentz(Ur2A,0); // this should be the same as U2iS + + typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Diff_check_mat2A = Ur02A - U2iAS; + std::cout << GridLogMessage << "Projections structure check group difference (Two Index anti-Symmetric): " << norm2(Diff_check_mat2A) << std::endl; + + } else { + std::cout << GridLogMessage << "Skipping Two Index anti-Symmetric tests " + "because representation is trivial (dim = 1)" + << std::endl; } - - TIndexRep.update_representation(UV2A); - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField UVr2A = TIndexRepA.U; // (U_f * V_f)_r - - TIndexRep.update_representation(U2A); - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2A = TIndexRepA.U; // U_r - - TIndexRep.update_representation(V2A); - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Vr2A = TIndexRepA.U; // V_r - - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2Vr2A(grid); - Ur2Vr2A = Zero(); - for (int mu = 0; mu < Nd; mu++) { - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Urmu2A = peekLorentz(Ur2A,mu); - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Vrmu2A = peekLorentz(Vr2A,mu); - pokeLorentz(Ur2Vr2A,Urmu2A*Vrmu2A, mu); - } - - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Diff_check2A = UVr2A - Ur2Vr2A; - std::cout << GridLogMessage << "Group structure SU("<::LatticeMatrix Ar_Asym(grid); - random(gridRNG,h_Asym); - h_Asym = real(h_Asym); - SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym); - - // Re-extract h_sym - SU3::LatticeAlgebraVector h_Asym2(grid); - SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym); - SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2; - std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl; - - - // Exponentiate - typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix U2iAS(grid); - U2iAS = expMat(Ar_Asym, 1.0, 16); - - typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix uno2iAS(grid); - uno2iAS = 1.0; - // Check matrix U2iS, must be real orthogonal - typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ucheck2iAS = U2iAS - conjugate(U2iAS); - std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iAS) - << std::endl; - - Ucheck2iAS = U2iAS * adj(U2iAS) - uno2iAS; - std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iAS) - << std::endl; - Ucheck2iAS = adj(U2iAS) * U2iAS - uno2iAS; - std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iAS) - << std::endl; - - - - // Construct the fundamental matrix in the group - SU3::LatticeMatrix Af_Asym(grid); - SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym); - SU3::LatticeMatrix Ufund2A(grid); - Ufund2A = expMat(Af_Asym, 1.0, 16); - SU3::LatticeMatrix UnitCheck2A(grid); - UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f; - std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A) - << std::endl; - UnitCheck2A = adj(Ufund2A) * Ufund2A - uno_f; - std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2A) - << std::endl; - - - // Tranform to the 2Index Sym representation - U = Zero(); // fill this with only one direction - pokeLorentz(U,Ufund2A,0); // the representation transf acts on full gauge fields - - TIndexRepA.update_representation(U); - Ur2A = TIndexRepA.U; // U_r - typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ur02A = peekLorentz(Ur2A,0); // this should be the same as U2iS - - typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Diff_check_mat2A = Ur02A - U2iAS; - std::cout << GridLogMessage << "Projections structure check group difference (Two Index anti-Symmetric): " << norm2(Diff_check_mat2A) << std::endl; - -} else { - std::cout << GridLogMessage << "Skipping Two Index anti-Symmetric tests " - "because representation is trivial (dim = 1)" - << std::endl; -} - - - - - - - - +#endif Grid_finalize(); } diff --git a/tests/core/Test_reunitarise.cc b/tests/core/Test_reunitarise.cc index 6644be1a..10423159 100644 --- a/tests/core/Test_reunitarise.cc +++ b/tests/core/Test_reunitarise.cc @@ -122,14 +122,15 @@ int main (int argc, char ** argv) std::cout << "Determinant defect before projection " < using namespace std; using namespace Grid; +#ifdef HAVE_LIME + template class LocalCoherenceLanczosScidac : public LocalCoherenceLanczos { @@ -249,3 +251,11 @@ int main (int argc, char ** argv) { Grid_finalize(); } +#else + +int main( void ) +{ + return 0 ; +} + +#endif // HAVE_LIME_H diff --git a/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc b/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc index ecd0e629..669a7b6d 100644 --- a/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc +++ b/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc @@ -36,7 +36,8 @@ Author: Peter Boyle using namespace std; using namespace Grid; - ; + +#ifdef HAVE_LIME template class LocalCoherenceLanczosScidac : public LocalCoherenceLanczos @@ -250,3 +251,11 @@ int main (int argc, char ** argv) { Grid_finalize(); } +#else + +int main( void ) +{ + return 0 ; +} + +#endif diff --git a/tests/solver/Test_coarse_even_odd.cc b/tests/solver/Test_coarse_even_odd.cc index dfbab747..c7127121 100644 --- a/tests/solver/Test_coarse_even_odd.cc +++ b/tests/solver/Test_coarse_even_odd.cc @@ -93,8 +93,16 @@ int main(int argc, char** argv) { // Setup of Dirac Matrix and Operator // ///////////////////////////////////////////////////////////////////////////// - LatticeGaugeField Umu(Grid_f); SU3::HotConfiguration(pRNG_f, Umu); - + LatticeGaugeField Umu(Grid_f); +#if (Nc==2) + SU2::HotConfiguration(pRNG_f, Umu); +#elif (defined Nc==3) + SU3::HotConfiguration(pRNG_f, Umu); +#elif (defined Nc==4) + SU4::HotConfiguration(pRNG_f, Umu); +#elif (defined Nc==5) + SU5::HotConfiguration(pRNG_f, Umu); +#endif RealD checkTolerance = (getPrecision::value == 1) ? 1e-7 : 1e-15; RealD mass = -0.30; diff --git a/tests/solver/Test_dwf_mrhs_cg.cc b/tests/solver/Test_dwf_mrhs_cg.cc index b912ba4f..4242f64b 100644 --- a/tests/solver/Test_dwf_mrhs_cg.cc +++ b/tests/solver/Test_dwf_mrhs_cg.cc @@ -34,6 +34,7 @@ using namespace Grid; int main (int argc, char ** argv) { +#ifdef HAVE_LIME typedef typename DomainWallFermionR::FermionField FermionField; typedef typename DomainWallFermionR::ComplexField ComplexField; typename DomainWallFermionR::ImplParams params; @@ -237,4 +238,5 @@ int main (int argc, char ** argv) } Grid_finalize(); +#endif // HAVE_LIME }