don't zero created fields

Author: jcosborn

src/field/fieldET.nim

@@ -137,9 +137,18 @@ proc new*[V:static[int],T](x:var FieldObj[V,T]; l:Layout[V]) =
   x.s.new(l.nSitesOuter)
   #fence()
   x.elemSize = sizeOf(T)
+proc newU*[V:static[int],T](x:var FieldObj[V,T]; l:Layout[V]) =
+  # remember to change newFieldArray if the following changes
+  x.l = l
+  x.s.newU(l.nSitesOuter)
+  #fence()
+  x.elemSize = sizeOf(T)
 proc new*[V:static[int],T](x:var Field[V,T]; l:Layout[V]) =
   x.new()
   new(x[], l)
+proc newU*[V:static[int],T](x:var Field[V,T]; l:Layout[V]) =
+  x.new()
+  newU(x[], l)
 proc new*[V:static[int],T](x:var FieldObj[V,T]; y:Field) = x.new(y.l)
 proc new*[V:static[int],T](x:var Field[V,T]; y:Field) = x.new(y.l)
 proc newField*[V:static[int],T](l:Layout[V]; t:typedesc[T]):Field[V,T] =
@@ -148,6 +157,10 @@ proc newOneOf*(x: Field): auto =
   var r: type(x)
   r.new(x.l)
   r
+proc newOneOfU*(x: Field): auto =
+  var r: type(x)
+  r.newU(x.l)
+  r
 template l*(x: FieldUnop): untyped = x.f1.l
 proc newOneOf*(x: FieldUnop): auto =
   var r: evalType(x)

src/gauge/gaugeUtils.nim

@@ -220,7 +220,7 @@ proc staples*[T,F,B](ss,uu,vv:openArray[T];ff:openArray[F];bb:openArray[B]) =
 
 proc plaq*[T](uu: openArray[T]): auto =
   mixin mul, load1, createShiftBufs
-  tic()
+  tic("plaq")
   template getIp(mu,nu: int): int = ((mu*(mu-1)) div 2) + nu
   let u = cast[ptr cArray[T]](unsafeAddr(uu[0]))
   let lo = u[0].l
@@ -233,17 +233,17 @@ proc plaq*[T](uu: openArray[T]): auto =
   for i in 0..<nd-1: sf[nd-1][i] = sf[i][i]
   let np = (nd*(nd-1)) div 2
   var pl = newSeq[float64](np)
-  toc("plaq setup")
+  toc("setup")
   threads:
-    tic()
+    tic("plaq threads")
     #var plt = newSeq[float64](np)
     var plt: array[6,float64]
     var umunu,unumu: type(load1(u[0][0]))
     for mu in 0..<nd:
       for nu in 0..<nd:
         if mu != nu:
           startSB(sf[mu][nu], u[mu][ix])
-    toc("plaq start shifts")
+    toc("start shifts")
     for ir in u[0]:
       for mu in 1..<nd:
         for nu in 0..<mu:
@@ -253,13 +253,13 @@ proc plaq*[T](uu: openArray[T]): auto =
             let ip = getIp(mu,nu)
             let dt = redot(umunu,unumu)
             plt[ip] += simdSum(dt)
-    toc("plaq local")
+    toc("local")
     var needBoundary = false
     for mu in 0..<nd:
       for nu in 0..<nd:
         if mu != nu:
           boundaryWaitSB(sf[mu][nu]): needBoundary = true
-    toc("plaq wait")
+    toc("wait")
     if needBoundary:
       boundarySyncSB()
       for ir in u[0]:
@@ -279,15 +279,15 @@ proc plaq*[T](uu: openArray[T]): auto =
               let ip = getIp(mu,nu)
               let dt = redot(umunu,unumu)
               plt[ip] += simdSum(dt)
-    toc("plaq boundary")
+    toc("boundary")
     threadSum(plt)
     if threadNum == 0:
       for i in 0..<pl.len:
         pl[i] = plt[i]/(lo.physVol.float*float(np*nc))
       rankSum(pl)
-    toc("plaq sum")
+    toc("sum")
   result = pl
-  toc("plaq end", flops=lo.nSites.float*float(np*(2*8*nc*nc*nc-1)))
+  toc("end", flops=lo.nSites.float*float(np*(2*8*nc*nc*nc-1)))
 
 template plaq*(g: Gauge): auto = plaq(g.u)
 
@@ -324,7 +324,7 @@ proc staples*[T,A,F,B](staples,uu,vv:openArray[T]; aa:openArray[A];
 
 proc plaq2*[T](gg:openArray[T]):auto =
   mixin adj
-  tic()
+  tic("plaq2")
   let g = cast[ptr cArray[T]](unsafeAddr(gg[0]))
   let lo = g[0].l
   let nd = lo.nDim
@@ -335,57 +335,61 @@ proc plaq2*[T](gg:openArray[T]):auto =
   var s1 = lo.ColorMatrix()
   #var t1 = lo.ColorMatrix()
   var tr:type(trace(m))
-  toc("plaq2 setup")
+  toc("setup")
   threads:
-    #tic()
+    tic("plaq2 threads")
     m := 0
-    #toc("plaq2 zero")
+    toc("zero")
     for mu in 1..<nd:
       for nu in 0..<mu:
-        #tic()
+        tic("plaq2 loop")
+        #toc("before shift1")
         shift(s0, mu,1, g[nu])
-        #toc("plaq2 shift1")
+        toc("shift1")
         shift(s1, nu,1, g[mu])
-        #toc("plaq2 shift2")
+        toc("shift2")
         #echo "s0: ", trace(s0)
         #echo "s1: ", trace(s1)
         m += (g[mu]*s0) * (g[nu]*s1).adj
         #m += (g[mu]*s0) * (g[nu]*s1)
         #echo mu, " ", nu, " ", trace(m)/nc
-        #toc("plaq2 mul")
-    #toc("plaq2 work")
+        toc("mul")
+    toc("work")
     tr = trace(m)
-    #toc("plaq2 trace")
-  toc("plaq2 threads")
+    toc("trace")
+  toc("end")
   result = tr/(lo.physVol.float*0.5*float(nd*(nd-1)*nc))
 
 proc plaq3*[T](g: seq[T]): auto =
   mixin adj, newTransporters
-  tic()
+  tic("plaq3")
   let lo = g[0].l
   let nd = lo.nDim
   let nc = g[0][0].ncols
   let t = newTransporters(g, g[0], 1)
   var m = lo.ColorMatrix()
   var tr: type(trace(m))
-  toc("plaq3 setup")
+  toc("setup")
   threads:
-    tic()
+    tic("plaq3 threads")
     m := 0
-    toc("plaq3 zero")
+    toc("zero")
     for mu in 1..<nd:
       for nu in 0..<mu:
-        tic()
+        tic("plaq3 loop")
         #m += (t[mu]^*g[nu]) * (t[nu]^*g[mu]).adj
-        discard t[mu]^*g[nu]
-        discard t[nu]^*g[mu]
+        discard t[mu]^*!g[nu]
+        toc("transport1")
+        discard t[nu]^*!g[mu]
+        threadBarrier()
+        toc("transport2")
         m += t[mu].field * t[nu].field.adj
         #echo mu, " ", nu, " ", trace(m)/nc
-        toc("plaq3 mul")
-    toc("plaq3 work")
+        toc("mul")
+    toc("work")
     tr = trace(m)
-    toc("plaq3 trace")
-  toc("plaq3 threads")
+    toc("trace")
+  toc("end")
   result = tr/(lo.physVol.float*0.5*float(nd*(nd-1)*nc))
 
 proc echoPlaq*(g: auto) =

src/gauge/symanzik1loopAction.nim

@@ -43,20 +43,13 @@ proc symanzik1loopAction*[T](coeffs: GaugeActionCoeffs;
     nflops = newThreadSingle(0)
     U = newSeq[GF](nd)
     Uf = newSeq2d[Shft](nd,nd)
-  #var mu: cint = 0
-  #while mu < nd:
   for mu in 0..<nd:
     U[mu] = links[mu]
-    #var nu: cint = 0
-    #while nu < nd:
     for nu in 0..<nd:
-      #if nu == mu: continue
-      #Uf[mu,nu] = links[0].newOneOf
-      #QDP_M_eq_sM(Uf[mu][nu], U[mu], neighbor[nu], QDP_forward, sub)
       if nu != mu:
+        #Uf[mu,nu] = links[0].newOneOf
+        #QDP_M_eq_sM(Uf[mu][nu], U[mu], neighbor[nu], QDP_forward, sub)
         Uf[mu,nu] = newShifter(U[mu], nu, 1)
-      #inc(nu)
-    #inc(mu)
   threads:
     for mu in 0..<nd:
       for nu in 0..<nd:
@@ -73,11 +66,7 @@ proc symanzik1loopAction*[T](coeffs: GaugeActionCoeffs;
     bstpl0 = newSeq2d[GF](nd,nd)
     bstpl = newSeq2d[Shft](nd,nd)
     tc = lo.Complex
-  #var mu: cint = 1
-  #while mu < nd:
   for mu in 1..<nd:
-    #var nu: cint = 0
-    #while nu < mu:
     for nu in 0..<mu:
       if pgm != 0.0:
         UUf[mu,nu] = links[0].newOneOf
@@ -242,68 +231,30 @@ proc symanzik1loopAction*[T](coeffs: GaugeActionCoeffs;
       nflops += 2 * EQMTM + 3 * PEQMTM + 4 * nc * nc
 
     if pgm != 0.0:
-      ##  FIXME: only works for nd=4
+      #  FIXME: only works for nd=4
       if nd != 4:
         qexError("symanzik1loopAction with parallelogram only works for nDim == 4")
       combinefb(pgms, 0, 1, 2)
       combineb(pgms, 0, 2, 1)
       combineb(pgms, 1, 2, 0)
-      ##  rest
-      ## combinefb(pgmt,0,3,1);
-      ## combinefb(pgmt,0,3,2);
+      #  rest
+      # combinefb(pgmt,0,3,1);
+      # combinefb(pgmt,0,3,2);
       combinefb2(pgmt, 0, 3, 1, 2)
       combineb(pgmt, 0, 1, 3)
       combineb(pgmt, 0, 2, 3)
       combinefb(pgmt, 1, 2, 3)
-      ## combineb(pgmt,1,3,0);
-      ## combineb(pgmt,1,3,2);
+      # combineb(pgmt,1,3,0);
+      # combineb(pgmt,1,3,2);
       combineb2(pgmt, 1, 3, 0, 2)
-      ## combineb(pgmt,2,3,0);
-      ## combineb(pgmt,2,3,1);
+      # combineb(pgmt,2,3,0);
+      # combineb(pgmt,2,3,1);
       combineb2(pgmt, 2, 3, 0, 1)
 
     threadSingle:
       acts = plaq * plaqs + rect * rects + pgm * pgms + adpl * adpls
       actt = plaq * plaqt + rect * rectt + pgm * pgmt + adpl * adplt
 
-#[
-  var mu: cint = 0
-  while mu < nd:
-    var nu: cint = 0
-    while nu < nd:
-      if nu == mu:
-        inc(nu)
-        inc(mu)
-        continue
-      QDP_destroy_M(Uf[mu,nu])
-      inc(nu)
-    inc(mu)
-  if pgm:
-    var mu: cint = 0
-    while mu < nd:
-      var nu: cint = 0
-      while nu < nd:
-        if nu == mu:
-          inc(nu)
-          inc(mu)
-          continue
-        QDP_destroy_M(UUf[mu,nu])
-        QDP_destroy_M(fstpl[mu,nu])
-        QDP_destroy_M(bstpl0[mu,nu])
-        QDP_destroy_M(bstpl[mu,nu])
-        inc(nu)
-      inc(mu)
-  else:
-    QDP_destroy_M(UUf[1,0])
-    QDP_destroy_M(UUf[0,1])
-    if rect:
-      QDP_destroy_M(fstpl[1,0])
-      QDP_destroy_M(fstpl[0,1])
-      QDP_destroy_M(bstpl0[1,0])
-      QDP_destroy_M(bstpl0[0,1])
-      QDP_destroy_M(bstpl[1,0])
-      QDP_destroy_M(bstpl[0,1])
- ]#
   let act0 = lo.physVol*(coeffs.plaq + 2*coeffs.rect + coeffs.adjplaq)
   let act1 = lo.physVol*(4*coeffs.pgm)
   result.space = 0.5*(nd-1)*(nd-2)*act0 + act1 - acts