graph/gauge: fix gaugeActionDeriv2 with proper zero init

src/experimental/graph/gauge.nim

@@ -40,7 +40,9 @@ method valCopy*(z: Ggauge, x: Ggauge) =
     for mu in 0..<u.len:
       u[mu] := v[mu]
 
-method `$`*(x: Ggauge): string = "Gauge"
+method `$`*(x: Ggauge): string =
+  let v = x.gval[0][0][0,0]
+  "Gauge (" & $v.re[0] & ", " & $v.im[0] & ")"
 
 #
 # basic ops
@@ -486,6 +488,14 @@ type Gactcoeff* {.final.} = ref object of Gvalue
 
 proc getactcoeff*(x: Gvalue): GaugeActionCoeffs = Gactcoeff(x).cval
 
+proc `getactcoeff=`*(x: Gvalue, c: GaugeActionCoeffs) =
+  let gc = Gactcoeff(x)
+  gc.cval = c
+
+proc update*(x: Gvalue, c: GaugeActionCoeffs) =
+  x.getactcoeff = c
+  x.updated
+
 converter toGvalue*(x: GaugeActionCoeffs): Gvalue =
   result = Gactcoeff(cval: x)
   result.updated
@@ -655,6 +665,9 @@ proc gaugeActionDeriv2f(v: Gvalue) =
   let g = Ggauge(v.inputs[2])
   let z = Ggauge(v)
   if gc.adjplaq == 0:
+    threads:
+      for mu in 0..<z.gval.len:
+        z.gval[mu] := 0.0
     gc.gaugeDerivDeriv2(g.gval, b.gval, z.gval)
   elif gc.rect == 0 and gc.pgm == 0:
     raiseValueError("unimplemented")

src/experimental/graph/tggauge.nim

@@ -41,7 +41,7 @@ proc ndiff(zt: Gvalue, t: Gscalar): (float, float) =
   (dzdt, e)
 
 template check(ii: tuple[filename:string, line:int, column:int], ast: string, dzdt, e, gdota: float) =
-  if not almostEqual(gdota, dzdt, unitsInLastPlace = 1024*1024):
+  if not almostEqual(gdota, dzdt, unitsInLastPlace = 512*1024):
     checkpoint(ii.filename & ":" & $ii.line & ":" & $ii.column & ": Check failed: " & ast)
     checkpoint("  ndiff: " & $dzdt & " +/- " & $e)
     checkpoint("  grad: " & $gdota)
@@ -289,12 +289,28 @@ suite "gauge action":
     let c = actWilson(beta)
     proc act(x: Gvalue): Gvalue = gaugeAction(c, x)
     proc force(x: Gvalue): Gvalue = gaugeForce(c, x)
-    ckforce(act, force, gg, 4.0*gm)
+    ckforce(act, force, gg, 10.0*gm)
 
   test "wilson force gradient":
     let beta = 5.4
     let c = actWilson(beta)
     proc force(x: Gvalue): Gvalue = gaugeForce(c, x)
     ckgradm(force, gg, gu, gm)
 
+  test "wilson force gradient recomp":
+    let beta = 5.4
+    let c = actWilson(beta)
+    let a = gaugeAction(c, gg)
+    let f2 = gaugeForce(c, gg).norm2
+    let df2 = grad(f2, gg).norm2
+    let rs1 = [a.eval.getfloat, f2.eval.getfloat, df2.eval.getfloat]
+    c.updated
+    gg.updated
+    let rs2 = [a.eval.getfloat, f2.eval.getfloat, df2.eval.getfloat]
+    c.updated
+    gg.updated
+    let rs3 = [a.eval.getfloat, f2.eval.getfloat, df2.eval.getfloat]
+    check rs1 == rs2
+    check rs1 == rs3
+
 qexFinalize()

src/experimental/graph/tgstout.nim

@@ -6,7 +6,7 @@ qexInit()
 letParam:
   lat = @[12,12,12,24]
   dt = 0.1
-  eps = 0.01
+  eps = 0.004
   nstep = 3
   beta = 5.4
   seed:uint = 1234567891
@@ -50,7 +50,7 @@ proc act(t: float): float =
   gc.gaugeAction1 u
 
 var ndt, err: float
-ndiff(ndt, err, act, dt, eps, ordMax=4)
+ndiff(ndt, err, act, dt, eps, ordMax=3)
 echo "numdiff smear dS/dt: ",ndt," +/- ",err
 
 proc stout(g, t: Gvalue, n: int): Gvalue =
@@ -91,8 +91,8 @@ echo "rel dS: ",rds
 echo "rel graph dS/dt: ",rgdt
 echo "rel ndiff dS/dt: ",rndt
 
-doassert rds < 1e-10
-doassert rgdt < 1e-10
-doassert rndt < 1e-10
+doassert rds < 1e-11
+doassert rgdt < 1e-11
+doassert rndt < 1e-11
 
 qexFinalize()