From 4924b3209e9bb5c358699ebfd8c3d929eaec87c2 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 23 Jan 2024 14:43:58 -0700 Subject: [PATCH] projectU3 yields a unitary matrix --- .gitignore | 4 ++++ Grid/qcd/smearing/HISQSmearing.h | 31 ++++++++++++++++++------------- tests/smearing/Test_fatLinks.cc | 8 +------- 3 files changed, 23 insertions(+), 20 deletions(-) diff --git a/.gitignore b/.gitignore index 40156f9d..94e866e2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,7 @@ +# Doxygen stuff +html/* +latex/* + # Compiled Object files # ######################### *.slo diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 17959495..c8255acc 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -330,7 +330,7 @@ public: // IN--u_mu void projectU3(GF& u_proj, GF& u_mu) const { - LF V(u_mu.Grid()), Q(u_mu.Grid()), sqrtQinv(u_mu.Grid()), id_3(u_mu.Grid()); + LF V(u_mu.Grid()), Q(u_mu.Grid()), sqrtQinv(u_mu.Grid()), id_3(u_mu.Grid()), diff(u_mu.Grid()); CF c0(u_mu.Grid()), c1(u_mu.Grid()), c2(u_mu.Grid()), g0(u_mu.Grid()), g1(u_mu.Grid()), g2(u_mu.Grid()), S(u_mu.Grid()), R(u_mu.Grid()), theta(u_mu.Grid()), u(u_mu.Grid()), v(u_mu.Grid()), w(u_mu.Grid()), den(u_mu.Grid()), f0(u_mu.Grid()), f1(u_mu.Grid()), @@ -338,18 +338,22 @@ public: // Follow MILC 10.1103/PhysRevD.82.074501, eqs (B2-B3) and (C1-C8) for (int mu = 0; mu < Nd; mu++) { - V = PeekIndex(u_mu, mu); - Q = adj(V)*V; - c0 = real(trace(Q)); - c1 = (1/2.)*real(trace(Q*Q)); - c2 = (1/3.)*real(trace(Q*Q*Q)); - S = (1/3.)*c1-(1/18.)*c0*c0; - R = (1/2.)*c2-(1/3. )*c0*c1+(1/27.)*c0*c0*c0; - theta = acos(R*pow(S,-1.5)); - g0 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta-2*M_PI/3.); - g1 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta ); - g2 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta+2*M_PI/3.); -// if (fabs(Q.determinant()/(g0*g1*g2)-1.0) > 1e-5) {} + V = PeekIndex(u_mu, mu); + Q = adj(V)*V; + c0 = real(trace(Q)); + c1 = (1/2.)*real(trace(Q*Q)); + c2 = (1/3.)*real(trace(Q*Q*Q)); + S = (1/3.)*c1-(1/18.)*c0*c0; + if (norm2(S)<1e-28) { + g0 = (1/3.)*c0; g1 = g0; g2 = g1; + } else { + R = (1/2.)*c2-(1/3. )*c0*c1+(1/27.)*c0*c0*c0; + theta = acos(R*pow(S,-1.5)); + g0 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta-2*M_PI/3.); + g1 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta ); + g2 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta+2*M_PI/3.); + } +// if (fabs(Q.determinant()/(g0*g1*g2)-1.0) > 1e-5) { SVD } u = sqrt(g0) + sqrt(g1) + sqrt(g2); v = sqrt(g0*g1) + sqrt(g0*g2) + sqrt(g1*g2); w = sqrt(g0*g1*g2); @@ -360,6 +364,7 @@ public: id_3 = 1.; sqrtQinv = f0*id_3 + f1*Q + f2*Q*Q; + PokeIndex(u_proj, V*sqrtQinv, mu); } }; diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index e2dc5d6d..742cb205 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -51,12 +51,6 @@ struct ConfParameters: Serializable { } }; -// -// one method: input --> fat -// another : input --> long (naik) -// another : input --> unitarize -// - void testSmear(GridCartesian& GRID, LatticeGaugeFieldD Umu, LatticeGaugeFieldD Usmr, LatticeGaugeFieldD Unaik, LatticeGaugeFieldD Ucontrol, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) { @@ -79,7 +73,7 @@ void testSmear(GridCartesian& GRID, LatticeGaugeFieldD Umu, LatticeGaugeFieldD U } else { Grid_error(" |Umu-Unaik|/|Umu| = ",absDiff); } - hisq_fat.projectU3(Uproj,Usmr); + hisq_fat.projectU3(Uproj,Ucontrol); // NerscIO::writeConfiguration(Unaik,"nersc.l8t4b3360.naik"); } }