Cleaned up hisq_force & merged with tests that James added to it

src/examples/hisq_force.nim

@@ -1,5 +1,5 @@
 # C.T. Peterson: force test inspired from conversation with Peter Boyle
-# See Grid implementation here:
+# See Grid implementation here: 
 #   -https://github.com/paboyle/Grid/blob/develop/tests/forces/Test_bdy.cc
 import qex
 import gauge/[hisqsmear]
@@ -34,7 +34,7 @@ spf.r2req = frsq
 spf.maxits = 10000
 spf.verbosity = 1
 
-# -- Generic
+# -- Generic 
 
 proc smearRephase(g: auto, sg,sgl: auto): auto {.discardable.} =
   tic()
@@ -66,7 +66,7 @@ proc action(): float =
 proc smearedOneAndThreeLinkForce(f: auto, smearedForce: proc, p: auto, g:auto) =
   # reverse accumulation of the derivative
   # 1. Dslash
-  var
+  var 
     f1 = f.newOneOf()
     f3 = f.newOneOf()
     ff = f.newOneOf()
@@ -76,6 +76,9 @@ proc smearedOneAndThreeLinkForce(f: auto, smearedForce: proc, p: auto, g:auto) =
     discard t[mu] ^* p
     t3[mu] = newShifter(p,mu,3)
     discard t3[mu] ^* p
+    #for i in 0..<4:
+    #  if (i == 0): discard t3[mu] ^* p
+    #  else: discard t3[mu] ^* t3[mu].field
   const n = p[0].len
   threads:
     for mu in 0..<f.len:
@@ -84,15 +87,19 @@ proc smearedOneAndThreeLinkForce(f: auto, smearedForce: proc, p: auto, g:auto) =
           forO b, 0, n-1:
             f1[mu][i][a,b] := p[i][a] * t[mu].field[i][b].adj
             f3[mu][i][a,b] := p[i][a] * t3[mu].field[i][b].adj
+            #f1[mu][i][a,b] := t[mu].field[i][a] * p[i][b].adj
+            #f3[mu][i][a,b] := t3[mu].field[i][a] * p[i][b].adj
 
   # 2. correcting phase
   threads:
-    f1.setBC; f3.setBC;
+    f1.setBC 
+    f3.setBC
     threadBarrier()
-    f1.stagPhase; f3.stagPhase;
+    f1.stagPhase 
+    f3.stagPhase
     threadBarrier()
     for mu in 0..<f.len:
-      for i in f[mu].odd:
+      for i in f[mu].odd: 
         f1[mu][i] *= -1
         f3[mu][i] *= -1
 
@@ -104,7 +111,8 @@ proc smearedOneAndThreeLinkForce(f: auto, smearedForce: proc, p: auto, g:auto) =
     for mu in 0..<f.len:
       for i in f[mu]:
         var s {.noinit.}: typeof(f[0][0])
-        s := ff[mu][i]*g[mu][i].adj
+        #s := g[mu][i]*ff[mu][i]
+        s := ff[mu][i] * g[mu][i].adj
         f[mu][i].projectTAH(s)
 
 proc fforce(f: auto) =
@@ -168,13 +176,11 @@ threads:
 var dS: float
 threads:
   var dSt = 0.0
-  for mu in 0..<p.len:
+  for mu in 0..<p.len: 
     dSt = dSt - reTrMul(p[mu],f[mu])
   threadMaster: dS = dSt
 
 # Compare differences
 let dH = s2+p2-s1-p1
 let (dSdt1,dSdt2) = (dS*eps,s2-s1)
-echo "dt*dS/dt, dS, difference = ", dSdt1,", ", dSdt2, ", ", dSdt1-dSdt2
-
-qexFinalize()
+echo "dt*dS/dt, dS, difference = ", dSdt1,", ", dSdt2, ", ", dSdt1-dSdt2

src/gauge/hisqsmear.nim

@@ -1,18 +1,61 @@
 import qex
-import physics/[hisqLinks]
 import gauge
 import gauge/[fat7l,fat7lderiv]
 
 export hisqLinks
 
+proc uadj[T](u: T): T =
+  var t = newOneOf(u)
+  threads:
+    for mu in 0..<t.len: t[mu] := adj(u[mu])
+  result = t
+
+proc asqtadDeriv[T](
+    deriv: auto, 
+    gauge,mid: T,  
+    coef: Fat7lCoefs,
+    llgauge,llmid: T, 
+    naik: float,
+    perf: var PerfInfo
+  ) =
+  var (f,fll) = (newOneOf(mid),newOneOf(mid))
+  fat7lderiv(f,gauge,mid,coef,fll,llgauge,llmid,naik,perf)
+  threads:
+    for mu in 0..<mid.len:
+      for s in deriv[mu]: deriv[mu][s] := f[mu][s] + fll[mu][s]
+
+proc fat7Deriv[T](
+    deriv: auto,
+    gauge,mid: T,
+    coef: Fat7lCoefs,
+    perf: var PerfInfo
+  ) = 
+  #var t = newOneOf(mid)
+  #threads:
+  #  for mu in 0..<t.len: t[mu] := adj(mid[mu])
+  #deriv.fat7lDeriv(gauge,t,coef,perf)
+  deriv.fat7lDeriv(gauge,mid,coef,perf)
+
+proc projectU[T](v: auto; u: T) =
+  threads:
+    for mu in 0..<u.len: 
+      for s in u[mu]: v[mu][s].projectU(u[mu][s])
+
+proc projectUDeriv[T](dvdu: auto; v,u: T; chain: T) =
+  threads:
+    for mu in 0..<chain.len:
+      for s in chain[mu]: 
+        #dvdu[mu][s].projectUderiv(v[mu][s],u[mu][s],adj(chain[mu][s]))
+        dvdu[mu][s].projectUderiv(v[mu][s],u[mu][s],chain[mu][s])
+
 proc newHISQ*(lepage: float = 0.0; naik: float = 1.0): HisqCoefs =
   result = HisqCoefs(naik: -naik/24.0)
   result.fat7first.setHisqFat7(lepage,0.0)
   result.fat7second.setHisqFat7(2.0-lepage,naik)
 
 proc smearGetForce*[T](
-    self: HisqCoefs;
-    u: T;
+    self: HisqCoefs; 
+    u: T; 
     su,sul: T;
     displayPerformance: bool = false
   ): proc(dsdu: var T; dsdsu,dsdsul: T) =
@@ -26,26 +69,23 @@ proc smearGetForce*[T](
     v = newOneOf(u)
     w = newOneOf(u)
     info: PerfInfo
-
+  
   # Smear
   v.makeImpLinks(u,fat7l1,info) # First fat7
-  threads: # Unitary projection
-    for mu in 0..<w.len:
-      for s in w[mu]: w[mu][s].projectU(v[mu][s])
+  w.projectU(v) # Unitary projection
   makeImpLinks(su,w,fat7l2,sul,w,naik,info) # Second fat7
+  #makeImpLinks(su,u,fat7l2,sul,u,naik,info) # Second fat7
 
   # Chain rule - retains a reference to u,su,sul
   proc smearedForce(dsdu: var T; dsdsu,dsdsul: T) =
-    var
-      dsdx_dxdw = newOneOf(dsdu)
-      dsdx_dxdw_dwdv = newOneOf(dsdu)
-    dsdx_dxdw.fat7lDeriv(w,dsdsu,fat7l2,w,dsdsul,naik,info) # Second fat7
-    threads: # Unitary projection
-      for mu in 0..<dsdx_dxdw_dwdv.len:
-        for s in dsdx_dxdw_dwdv[mu]:
-          dsdx_dxdw_dwdv[mu][s].projectUderiv(w[mu][s],v[mu][s],dsdx_dxdw[mu][s])
-    dsdu.fat7lDeriv(u,dsdx_dxdw_dwdv,fat7l1,info) # First fat7
-
+    var t = newOneOf(dsdu)
+    #dsdu.asqtadDeriv(u,dsdsu,fat7l2,u,dsdsul,naik,info) # Second fat7
+    t.asqtadDeriv(w,dsdsu,fat7l2,w,dsdsul,naik,info) # Second fat7
+    t.projectUDeriv(w,v,t) # Unitary projection
+    dsdu.fat7Deriv(u,t,fat7l1,info) # First fat7
+    if displayPerformance: echo $(info)
+
+  # Display performance (if requested) and return
   if displayPerformance: echo $(info)
   return smearedForce
 
@@ -152,4 +192,4 @@ when isMainModule:
   checkL("all", 1)
 
   echoProf()
-  qexFinalize()
+  qexFinalize()