diff --git a/Grid/qcd/utils/Sp2n.h b/Grid/qcd/utils/Sp2n.h index 1dafd771..ad6933e3 100644 --- a/Grid/qcd/utils/Sp2n.h +++ b/Grid/qcd/utils/Sp2n.h @@ -433,7 +433,85 @@ public: ColdConfiguration(out); } + static void OmegaInvariance(ColourMatrix &in) + { + + ColourMatrix Omega; + Omega = Zero(); + + for (int i = 0; i < ncolour; i++) + { + Omega()()(i, 2*ncolour-1-i) = 1.; + Omega()()(2*ncolour-1-i, i) = -1; + } + + auto diff = Omega - (in * Omega * transpose(in) ); + auto sdiff = norm2(diff); + if (norm2(sdiff) < 1e-8) + { + std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl; + } else { + std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl; + } + + } + template + static void OmegaInvariance(GaugeField &in) + { + typedef typename GaugeField::vector_type vector_type; + typedef iSp2nMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType U(in.Grid()); + LatticeMatrixType OmegaLatt(in.Grid()); + LatticeMatrixType identity(in.Grid()); + RealD vol = in.Grid()->gSites(); + ColourMatrix Omega; + + OmegaLatt = Zero(); + Omega = Zero(); + identity = 1.; + + U = PeekIndex(in,1); + + OmegaInvariance(U); + + } + + static void OmegaInvariance(LatticeColourMatrixD &in) + { + + LatticeColourMatrixD OmegaLatt(in.Grid()); + LatticeColourMatrixD identity(in.Grid()); + RealD vol = in.Grid()->gSites(); + ColourMatrix Omega; + + OmegaLatt = Zero(); + Omega = Zero(); + identity = 1.; + + for (int i = 0; i < ncolour; i++) + { + Omega()()(i, 2*ncolour-1-i) = 1.; + Omega()()(2*ncolour-1-i, i) = -1; + } + + OmegaLatt = OmegaLatt + (identity*Omega); + + auto diff = OmegaLatt - (in*Omega*transpose(in)); + auto sdiff = sum(diff); + //assert( norm2(sdiff) < 1e-8 ); + if (norm2(sdiff) < 1e-8) + { + std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl; + } else { + std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl; + } + + + } + }; // end of class Sp diff --git a/configure.ac b/configure.ac index caadca41..c6a85588 100644 --- a/configure.ac +++ b/configure.ac @@ -787,7 +787,6 @@ os (target) : $target_os compiler vendor : ${ax_cv_cxx_compiler_vendor} compiler version : ${ax_cv_gxx_version} ----- BUILD OPTIONS ----------------------------------- -#gauge group : `if test "${ac_ENABLE_SP}x" == "yesx"; then echo Sp2n; else echo SUn; fi` Nc : ${ac_Nc} SIMD : ${ac_SIMD}${SIMD_GEN_WIDTH_MSG} Threading : ${ac_openmp} diff --git a/tests/sp2n/Test_Sp_start.cc b/tests/sp2n/Test_Sp_start.cc index eaeb59f4..da19779f 100644 --- a/tests/sp2n/Test_Sp_start.cc +++ b/tests/sp2n/Test_Sp_start.cc @@ -14,54 +14,79 @@ int main (int argc, char **argv) GridRedBlackCartesian RBGrid(&Grid); LatticeGaugeField Umu(&Grid); - LatticeGaugeField identity(&Grid); LatticeColourMatrixD U(&Grid); - LatticeColourMatrixD Cidentity(&Grid); LatticeColourMatrixD aux(&Grid); - - identity = 1.0; - Cidentity = 1.0; - Complex i(0., 1.); + LatticeGaugeField identity(&Grid); + LatticeColourMatrixD Cidentity(&Grid); double vol = Umu.Grid()->gSites(); - + const int nsp = Nc / 2; + identity = 1.; + Cidentity = 1.; + std::vector pseeds({1,2,3,4,5}); std::vector sseeds({6,7,8,9,10}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(pseeds); GridSerialRNG sRNG; sRNG.SeedFixedIntegers(sseeds); std::cout << GridLogMessage << std::endl; - - + std::cout << GridLogMessage << "-------" << std::endl; + std::cout << GridLogMessage << "Checking Cold Configuration " << std::endl; + std::cout << GridLogMessage << std::endl; Sp::ColdConfiguration(pRNG,Umu); + Sp::OmegaInvariance(Umu); + std::cout << GridLogMessage << std::endl; + Umu = Umu - identity; U = PeekIndex(Umu,1); std::cout << GridLogMessage << "U - 1 = " << norm2(sum(U)) << std::endl; assert ( norm2(sum(U)) < 1e-6); std::cout << GridLogMessage << std::endl; - - std::cout << GridLogMessage << "Hot Start " << std::endl; + std::cout << GridLogMessage << "-------" << std::endl; + std::cout << GridLogMessage << "Checking Hot Configuration" << std::endl; + std::cout << GridLogMessage << std::endl; Sp::HotConfiguration(pRNG,Umu); + Sp::OmegaInvariance(Umu); + std::cout << GridLogMessage << std::endl; + U = PeekIndex(Umu,1); + std::cout << GridLogMessage << "Checking unitarity " << std::endl; + aux = U*adj(U) - Cidentity; + assert( norm2(aux) < 1e-6 ); + std::cout << GridLogMessage << "ok" << std::endl; + std::cout << GridLogMessage << std::endl; + std::cout << GridLogMessage << "Checking determinant " << std::endl; + std::cout << GridLogMessage << "Det = " << norm2( Determinant(U) )/vol << std::endl; + std::cout << GridLogMessage << std::endl; + std::cout << GridLogMessage << "-------" << std::endl; + std::cout << GridLogMessage << "Checking Tepid Configuration" << std::endl; + std::cout << GridLogMessage << std::endl; + Sp::TepidConfiguration(pRNG,Umu); + Sp::OmegaInvariance(Umu); + std::cout << GridLogMessage << std::endl; U = PeekIndex(Umu,1); - std::cout << GridLogMessage << "Checking unitarity " << std::endl; - aux = U*adj(U) - Cidentity; - assert( norm2(aux) < 1e-6 ); - std::cout << GridLogMessage << "ok" << std::endl; std::cout << GridLogMessage << std::endl; - std::cout << GridLogMessage << "Checking determinant " << std::endl; std::cout << GridLogMessage << "Det = " << norm2( Determinant(U) )/vol << std::endl; std::cout << GridLogMessage << std::endl; - std::cout << GridLogMessage << "Checking the matrix is in Sp2n " << std::endl; - std::cout << GridLogMessage << std::endl; + + std::cout << GridLogMessage << "Checking the structure is " << std::endl; + std::cout << GridLogMessage << "U = ( W X ) " << std::endl; + std::cout << GridLogMessage << " ( -X^* W^* ) " << std::endl; + std::cout << GridLogMessage << std::endl; + + Complex i(0., 1.); + + + Sp::HotConfiguration(pRNG,Umu); + U = PeekIndex(Umu,0); for (int c1 = 0; c1 < nsp; c1++) //check on W { for (int c2 = 0; c2 < nsp; c2++) @@ -71,17 +96,14 @@ int main (int argc, char **argv) auto Ww = conjugate( Wstar ); auto amizero = sum(W - Ww); auto amizeroo = TensorRemove(amizero); - std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl; + //std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl; assert( amizeroo.real() < 10e-6 ); amizeroo *= i; assert( amizeroo.real() < 10e-6 ); - std::cout << GridLogMessage << "ok " << std::endl; + //std::cout << GridLogMessage << "ok " << std::endl; } } - std::cout < ppseeds({11,12,13,14,15}); - std::vector ssseeds({16,17,18,19,20}); - GridParallelRNG ppRNG(&Grid); pRNG.SeedFixedIntegers(ppseeds); - GridSerialRNG ssRNG; sRNG.SeedFixedIntegers(ssseeds); - - std::cout << GridLogMessage << std::endl; - std::cout << GridLogMessage << "Testing gauge implementation " << std::endl; - - SymplPeriodicGimplR::HotConfiguration(ppRNG, Umu); - - U = PeekIndex(Umu,1); - + Sp::HotConfiguration(pRNG,Umu); + U = PeekIndex(Umu,0); for (int c1 = 0; c1 < nsp; c1++) //check on W { for (int c2 = 0; c2 < nsp; c2++) @@ -121,17 +134,14 @@ int main (int argc, char **argv) auto Ww = conjugate( Wstar ); auto amizero = sum(W - Ww); auto amizeroo = TensorRemove(amizero); - std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl; + //std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl; assert( amizeroo.real() < 10e-6 ); amizeroo *= i; assert( amizeroo.real() < 10e-6 ); - std::cout << GridLogMessage << "ok " << std::endl; + //std::cout << GridLogMessage << "ok " << std::endl; } } - std::cout <::HotConfiguration(pRNG,Umu); U = PeekIndex(Umu,2); - // it is unitary aux = U*adj(U) - identity; std::cout <::OmegaInvariance(U); std::cout <::OmegaInvariance(U); std::cout <::OmegaInvariance(U); + // checks on elements - // an Sp2n matrix is a block matrix of the type - // - // - // U = ( W X ) - // ( -X^* W^* ) - // - // will check this - // - std::cout << GridLogMessage << "Checking that U[ i , j ] - conj U[ i+N/2, j+N/2 ] = 0 for i = 0, ... , N/2 " << std::endl; + + std::cout << GridLogMessage << "Checking the structure is " << std::endl; + std::cout << GridLogMessage << "U = ( W X ) " << std::endl; + std::cout << GridLogMessage << " ( -X^* W^* ) " << std::endl; std::cout <