From 08188caab1f5825d8d0ffd4b96a1e8aca0e5c250 Mon Sep 17 00:00:00 2001
From: ctpeterson <curtistaylorpetersonwork@gmail.com>
Date: Thu, 19 Dec 2024 14:44:29 -0500
Subject: [PATCH] Added minor modifications to fat7lderiv to make it compatible
 with cleaned up hisqsmear.nim

---
 src/gauge/fat7lderiv.nim | 240 ++++++++++++++++++++++++++++-----------
 1 file changed, 174 insertions(+), 66 deletions(-)

diff --git a/src/gauge/fat7lderiv.nim b/src/gauge/fat7lderiv.nim
index 70a0a9c..6ecee7d 100644
--- a/src/gauge/fat7lderiv.nim
+++ b/src/gauge/fat7lderiv.nim
@@ -4,6 +4,121 @@ import base/profile, layout
 let STAPLE_FLOPS = 3*198+216+18
 let SIDE_FORCE_FLOPS = 7*198+3*216+3*18
 
+#[
+proc staple(auto out, auto in0, auto link0, int mu, int nu) =
+#define link     ftmp0[0]
+#define linkmu   ftmp[0][mu]
+#define in       ftmp0[1]
+#define innu     ftmp[1][nu]
+#define linkin   mtmp[0]
+#define back     btmp0[0]
+#define backnu   btmp[0][nu]
+#define linkinnu mtmp[1]
+
+  QDP_M_eq_M(link, link0, QDP_all);
+  QDP_M_eq_sM(linkmu, link, QDP_neighbor[mu], QDP_forward, QDP_all);
+  QDP_M_eq_M(in, in0, QDP_all);
+  QDP_M_eq_sM(innu, in, QDP_neighbor[nu], QDP_forward, QDP_all);
+  QDP_M_eq_Ma_times_M(linkin, link, in, QDP_all);
+  QDP_M_eq_M_times_M(back, linkin, linkmu, QDP_all);
+  QDP_M_eq_sM(backnu, back, QDP_neighbor[nu], QDP_backward, QDP_all);
+  QDP_M_eq_M_times_M(linkinnu, link, innu, QDP_all);
+  QDP_discard_M(innu);
+  QDP_M_peq_M_times_Ma(out, linkinnu, linkmu, QDP_all);
+  QDP_discard_M(linkmu);
+  QDP_M_peq_M(out, backnu, QDP_all);
+  QDP_discard_M(backnu);
+#define STAPLE_FLOPS (3*198+216+18)
+
+#undef link
+#undef linkmu
+#undef in
+#undef innu
+#undef linkin
+#undef back
+#undef backnu
+#undef linkinnu
+}
+]#
+
+#[
+static void
+side_force(QDP_ColorMatrix *force, QDP_ColorMatrix *bot0, QDP_ColorMatrix *side0,
+           QDP_ColorMatrix *top0, int mu, int nu, QDP_ColorMatrix *stpl)
+{
+#define side     ftmp0[0]
+#define sidemu   ftmp[0][mu]
+#define top      ftmp0[1]
+#define topnu    ftmp[1][nu]
+#define bot      ftmp0[2]
+#define botnu    ftmp[2][nu]
+#define sidebot  mtmp[0]
+#define sidetop  mtmp[1]
+#define topnusidebot  btmp0[0]
+#define fbmu          btmp[0][mu]
+#define botnusidetop  btmp0[1]
+#define fmbmu         btmp[1][mu]
+#define sidebotsidemu btmp0[2]
+#define stm           btmp[2][nu]
+#define botnusidemu   mtmp[2]
+#define botsidemu     mtmp[3]
+
+  // force += bot * sidemu * topnu+
+  // force -= bot-mu+ * side-mu * topnu-mu
+  // -= top <-> bot
+  // stpl += side * botnu * sidemu+
+  // stpl += side-nu+ * bot-nu * sidemu-nu
+  QDP_M_eq_M(side, side0, QDP_all);
+  QDP_M_eq_sM(sidemu, side, QDP_neighbor[mu], QDP_forward, QDP_all);
+  QDP_M_eq_M(top, top0, QDP_all);
+  QDP_M_eq_sM(topnu, top, QDP_neighbor[nu], QDP_forward, QDP_all);
+  QDP_M_eq_M(bot, bot0, QDP_all);
+  QDP_M_eq_sM(botnu, bot, QDP_neighbor[nu], QDP_forward, QDP_all);
+  QDP_M_eq_Ma_times_M(sidebot, side, bot, QDP_all);
+  QDP_M_eq_Ma_times_M(sidetop, side, top, QDP_all);
+  QDP_M_eq_Ma_times_M(topnusidebot, topnu, sidebot, QDP_all);
+  QDP_M_eq_sM(fbmu, topnusidebot, QDP_neighbor[mu], QDP_backward, QDP_all);
+  QDP_M_eq_Ma_times_M(botnusidetop, botnu, sidetop, QDP_all);
+  QDP_M_eq_sM(fmbmu, botnusidetop, QDP_neighbor[mu], QDP_backward, QDP_all);
+  QDP_M_eq_M_times_M(sidebotsidemu, sidebot, sidemu, QDP_all);
+  QDP_M_eq_sM(stm, sidebotsidemu, QDP_neighbor[nu], QDP_backward, QDP_all);
+  QDP_M_eq_M_times_Ma(botnusidemu, botnu, sidemu, QDP_all);
+  QDP_discard_M(botnu);
+  QDP_M_peq_M_times_M(stpl, side, botnusidemu, QDP_all);
+ //QDP_M_meq_M_times_Ma(force, top, botnusidemu, QDP_all);
+  QDP_M_peq_M_times_Ma(force, top, botnusidemu, QDP_all);
+  QDP_M_eq_M_times_M(botsidemu, bot, sidemu, QDP_all);
+  QDP_discard_M(sidemu);
+  QDP_M_peq_M_times_Ma(force, botsidemu, topnu, QDP_all);
+  QDP_discard_M(topnu);
+  //QDP_M_meq_Ma(force, fbmu, QDP_all);
+  QDP_M_peq_Ma(force, fbmu, QDP_all);
+  QDP_discard_M(fbmu);
+  QDP_M_peq_Ma(force, fmbmu, QDP_all);
+  QDP_discard_M(fmbmu);
+  QDP_M_peq_M(stpl, stm, QDP_all);
+  QDP_discard_M(stm);
+#define SIDE_FORCE_FLOPS (7*198+3*216+3*18)
+
+#undef side
+#undef sidemu
+#undef top
+#undef topnu
+#undef bot
+#undef botnu
+#undef sidebot
+#undef sidetop
+#undef topnusidebot
+#undef fbmu
+#undef botnusidetop
+#undef fmbmu
+#undef sidebotsidemu
+#undef stm
+#undef botnusidemu
+#undef botsidemu
+}
+]#
+
 proc fat7lDeriv*(deriv: auto, gauge: auto, mid: auto, coef: Fat7lCoefs,
                 llderiv: auto, llgauge: auto, llmid: auto, naik: float,
                 perf: var PerfInfo) =
@@ -18,7 +133,7 @@ proc fat7lDeriv*(deriv: auto, gauge: auto, mid: auto, coef: Fat7lCoefs,
   let have3 = (coef3!=0.0) or have5
   type lcm = type(gauge[0])
   let lo = gauge[0].l
-  proc newlcm: lcm = result.new(lo)
+  proc newlcm: lcm = result.new(gauge[0].l)
   var
     stpl3: array[4,lcm]
     mid5: array[4,lcm]
@@ -27,8 +142,8 @@ proc fat7lDeriv*(deriv: auto, gauge: auto, mid: auto, coef: Fat7lCoefs,
     s1a,s1b: array[4,Shifter[lcm,type(gauge[0][0])]]
     tm1,tm2: lcm
     sm1: array[4,Shifter[lcm,type(gauge[0][0])]]
-  if have3 or (naik!=0):
-    if have5 or (naik!=0):
+  if have3:
+    if have5:
       for mu in 0..3:
         stpl3[mu] = newlcm()
         mid5[mu] = newlcm()
@@ -178,7 +293,7 @@ proc fat7lDeriv*(deriv: auto, gauge: auto, mid: auto, coef: Fat7lCoefs,
     if naik!=0.0:
       for mu in 0..3:
         # force += mid * Umumu' * Umu'
-        # force += U-mu' * mid-mu * Umu'
+        # force -= U-mu' * mid-mu * Umu'
         # force += U-mu' * U-mu-mu' * mid-mu-mu
         #QDP_M_eq_M(Uf, U, QDP_all);
         #QDP_M_eq_sM(Umu, Uf, QDP_neighbor[mu], QDP_forward, QDP_all);
@@ -214,23 +329,6 @@ proc fat7lDeriv*(deriv: auto, gauge: auto, mid: auto, coef: Fat7lCoefs,
   perf.flops += nflops * gauge[0].l.localGeom.prod
   perf.secs += getElapsedTime()
 
-proc fat7lDeriv*(
-    deriv: auto,
-    gauge: auto,
-    mid: auto,
-    coef: Fat7lCoefs,
-    lgauge: auto,
-    lmid: auto,
-    naik: float,
-    perf: var PerfInfo
-  ) =
-  var (fx,fxl) = (newOneOf(mid),newOneOf(lmid))
-  fat7lDeriv(fx,gauge,mid,coef,fxl,lgauge,lmid,naik,perf)
-  threads:
-    for mu in 0..<deriv.len:
-      for s in deriv[mu]:
-        deriv[mu][s] += fx[mu][s] + fxl[mu][s]
-
 proc fat7lDeriv*(deriv: auto, gauge: auto, mid: auto, coef: Fat7lCoefs,
                 perf: var PerfInfo) =
   fat7lDeriv(deriv, gauge, mid, coef, deriv, gauge, mid, 0.0, perf)
@@ -267,8 +365,7 @@ when isMainModule:
   coef.lepage = 0.0
 
   #g.unit
-  #g.random r
-  g.gaussian r
+  g.random r
   #dg.randomTAH r
   dg.gaussian r
   for mu in 0..<g.len:
@@ -296,8 +393,6 @@ when isMainModule:
     gc.gaugeDeriv2(g, fd)
     check(a2-a, 1)
 
-  checkG()
-
   proc checkS(name: string, tol: float) =
     echo "Checking ", name
     for mu in 0..<fd.len:
@@ -313,9 +408,11 @@ when isMainModule:
     let a = gc.gaugeAction2(fl)
     let a2 = gc.gaugeAction2(fl2)
     gc.gaugeDeriv2(fl, ch)
-    fat7lDeriv(fd, g, ch, coef, info)
+    fat7lderiv(fd, g, ch, coef, info)
     check(a2-a, tol)
 
+  checkG()
+
   checkS("oneLink", 1)
   coef.oneLink = 0.0
   coef.threeStaple = 1.0
@@ -325,7 +422,7 @@ when isMainModule:
   checkS("fiveStaple", 2)
   coef.fiveStaple = 0.0
   coef.sevenStaple = 1.0
-  checkS("sevenStaple", 5)
+  checkS("sevenStaple", 3)
   coef.sevenStaple = 0.0
   coef.lepage = 1.0
   checkS("Lepage", 25)
@@ -336,48 +433,59 @@ when isMainModule:
   coef.lepage = 0.0
   checkS("not Lepage", 60)
   coef.lepage = 1.0
-  checkS("all fat", 90)
+  var naik = 1.0
 
-  proc checkL(name: string, tol: float) =
-    echo "Checking ", name
-    for mu in 0..<fd.len:
-      fd[mu] := 0
-      ch[mu] := 0
-      ld[mu] := 0
-    info.clear
-    resetTimers()
-    makeImpLinks(fl, g, coef, ll, g, naik, info)
-    echo "  ", info
-    makeImpLinks(fl2, g2, coef, ll2, g2, naik, info)
-    echo "  ", info
-    let a = gc.gaugeAction2(fl) + gc.gaugeAction2(ll)
-    let a2 = gc.gaugeAction2(fl2) + gc.gaugeAction2(ll2)
-    gc.gaugeDeriv2(fl, ch)
-    gc.gaugeDeriv2(ll, ld)
-    fat7lDeriv(fd, g, ch, coef, g, ld, naik, info)
-    #let bias = fl.norm2 / (fl.len*lo.physVol)
-    #let lbias = ll.norm2 / (ll.len*lo.physVol)
-    #let a = 0.5*(fl.norm2subtract(bias) + ll.norm2subtract(lbias))
-    #let a2 = 0.5*(fl2.norm2subtract(bias) + ll2.norm2subtract(lbias))
-    #echo a, " ", a2
-    #fat7lDeriv(fd, g, fl, coef, g, ll, naik, info)
-    check(a2-a, tol)
+  var fl = lo.newGauge()
+  var fl2 = lo.newGauge()
+  var ll = lo.newGauge()
+  var ll2 = lo.newGauge()
+  var dfl = lo.newGauge()
+  var g2 = lo.newGauge()
+  var dg = lo.newGauge()
+  var fd = lo.newGauge()
+  var ld = lo.newGauge()
+  for mu in 0..<dg.len:
+    dg[mu] := 0.00001 * g[mu]
+    g2[mu] := g[mu] + dg[mu]
+    fd[mu] := 0
 
-  coef.oneLink = 0.0
-  coef.threeStaple = 0.0
-  coef.fiveStaple = 0.0
-  coef.sevenStaple = 0.0
-  coef.lepage = 0.0
-  naik = 0.5
-  checkL("long link", 1)
+  echo g.plaq
+  echo g2.plaq
+  makeImpLinks(fl, g, coef, info)
+  info.clear
+  resetTimers()
+  makeImpLinks(fl, g, coef, info)
+  echo info
+  makeImpLinks(fl2, g2, coef, info)
+  echo info
+  echo fl.plaq
+  echo fl2.plaq
 
-  coef.oneLink = 1.0
-  coef.threeStaple = 1.0
-  coef.fiveStaple = 1.0
-  coef.sevenStaple = 1.0
-  coef.lepage = 1.0
-  naik = 100.0
-  checkL("all", 6)
+  fat7lderiv(fd, g, dg, coef, info)
+  echo info
+  for mu in 0..3:
+    dfl[mu] := fl2[mu] - fl[mu]
+    #echo dfl[mu].norm2
+    #echo fd[mu].norm2
+    dfl[mu] -= fd[mu]
+    echo dfl[mu].norm2
+
+  for mu in 0..<dg.len:
+    fd[mu] := 0
+    ld[mu] := 0
+  makeImpLinks(fl, g, coef, ll, g, naik, info)
+  makeImpLinks(fl2, g2, coef, ll2, g2, naik, info)
+  fat7lderiv(fd, g, dg, coef, ld, g, dg, naik, info)
+  echo info
+  for mu in 0..3:
+    dfl[mu] := fl2[mu] - fl[mu]
+    #echo dfl[mu].norm2
+    #echo fd[mu].norm2
+    dfl[mu] -= fd[mu]
+    echo dfl[mu].norm2
+    dfl[mu] := ll2[mu] - ll[mu]
+    dfl[mu] -= ld[mu]
+    echo dfl[mu].norm2
 
   echoProf()
   qexFinalize()