fix OMP pragma

update test tolerances
make staggered solver work with zero mass
update Grid wrapper

src/base/omp.nim

@@ -1,6 +1,7 @@
 import os
 
 when defined(noOpenmp):
+  static: echo "OpenMP disabled"
   template omp_set_num_threads*(x: cint) = discard
   template omp_get_num_threads*(): cint = 1
   template omp_get_max_threads*(): cint = 1
@@ -10,6 +11,7 @@ when defined(noOpenmp):
     block:
       body
 else:
+  static: echo "Using OpenMP"
   when existsEnv("OMPFLAG"):
     const ompFlag = getEnv("OMPFLAG")
   else:
@@ -25,7 +27,7 @@ else:
   template ompPragma(p:string):untyped =
     #forceOmpOn()
     #{. emit:["#pragma omp ", p] .}
-    {. emit:["_Pragma(\"omp ", p, "\");"] .}
+    {. emit:["_Pragma(\"omp ", p, "\")"] .}
   template ompBlock*(p:string; body:untyped):untyped =
     #{. emit:"#pragma omp " & p .}
     #{. emit:"{ /* Inserted by ompBlock " & p & " */".}

src/gauge/gaugeUtils.nim

@@ -920,6 +920,7 @@ proc lrmul(res:auto, l:auto, r:auto, la,ra:bool) =
 
 proc gaugeProd*(g:auto, ptree:OrdPathTree, origin=true):auto =
   tic("gaugeProd")
+  GC_fullCollect() # need to free space since we might allocate many fields
   type
     F = typeof(g[0])
     S = typeof(g[0][0])

src/gauge/wflow.nim

@@ -18,7 +18,7 @@ Zi = eps Z(Wi)
 
 ]#
 
-template gaugeFlow*(g: array|seq, steps: int, eps: float, measure: untyped): untyped =
+template gaugeFlow*(g: array|seq, steps: int, eps: float, measure: untyped) =
   ## Wilson flow.
   ## The input gauge field will be modified.
   ## `wflowT` is injected for `measure`.
@@ -65,7 +65,7 @@ template gaugeFlow*(g: array|seq, steps: int, eps: float, measure: untyped): unt
         break
     toc("end")
   flowProc()
-template gaugeFlow*(g: array|seq, eps: float, measure: untyped): untyped =
+template gaugeFlow*(g: array|seq, eps: float, measure: untyped) =
   gaugeFlow(g, 0, eps):
     measure
 

src/grid/Grid.nim

@@ -13,3 +13,6 @@ else:  # put stubs here
 
   proc gridSolveEE*(s:Staggered; r,t:Field; m:SomeNumber; sp: var SolverParams) =
     qexError "Grid not compiled in (define gridDir)"
+
+  proc gridSolveOO*(s:Staggered; r,t:Field; m:SomeNumber; sp: var SolverParams) =
+    qexError "Grid not compiled in (define gridDir)"

src/grid/GridDefs.nim

@@ -1,5 +1,5 @@
-import ospaths
-import io/qio
+import os
+#import io/qio
 
 const gridDir {.strdefine.} = getHomeDir() & "/lqcd/install/grid"
 const gridPassC = "-I" & gridDir / "include"
@@ -30,9 +30,9 @@ type
   GridWilsonLoops*[T] {.importcpp:"Grid::WilsonLoops",gh.} = object
   GridFermion*[T] {.importcpp:"'0::FermionField",gh,byref.} = object
   GridNaiveStaggeredFermionR* {.
-    importcpp:"Grid::NaiveStaggeredFermionR",gh.} = object
+    importcpp:"Grid::NaiveStaggeredFermionD",gh.} = object
   GridImprovedStaggeredFermionR* {.
-    importcpp:"Grid::ImprovedStaggeredFermionR",gh.} = object
+    importcpp:"Grid::ImprovedStaggeredFermionD",gh.} = object
 
 type
   GridParity* = enum
@@ -51,7 +51,7 @@ proc newStdVector*[R,T](a,b: ptr T): stdvector[R] {.
   constructor,importcpp:"'0(#,#)".}
 
 template `+`*[T](a: ptr T, n: int): untyped =
-  cast[ptr T](cast[ByteAddress](a) + n*sizeof(T))
+  cast[ptr T](cast[int](a) + n*sizeof(T))
 
 proc newStdVector*[T](a: openArray[T]): stdvector[T] =
   let a0 = unsafeaddr a[0]
@@ -108,7 +108,7 @@ template vector_obj*[T:GridFermion](x: T): untyped =
 template scalar_obj*[T:GridFermion](x: T): untyped =
   GridScalarObject[T]
 
-proc vectorizeFromLexOrdArray*[T](x: stdVector[T], r: any) {.
+proc vectorizeFromLexOrdArray*[T](x: stdVector[T], r: auto) {.
   importc,gh.}
 proc vectorizeFromLexOrdArray*[T](x: stdVector[T], r: GridLatticeGaugeField) {.
   importcpp:"vectorizeFromLexOrdArray(#,#)",gh.}
@@ -130,11 +130,11 @@ proc checkerboard*(x: ptr GridFermion): int =
   var r = 0
   #{.emit:"auto t = x->Checkerboard();".}
   #{.emit:"if(t==Grid::Odd){r=1;}".}
-  {.emit:"r = x->Checkerboard();".}
+  {.emit:[r," = ",x,"->Checkerboard();"].}
   r
 
 proc checkerboard*(x: ptr GridFermion, y: GridParity) =
-  {.emit:"x->Checkerboard() = y;".}
+  {.emit:[x,"->Checkerboard() = ",y,";"].}
 template checkerboard*(x: var GridFermion, y: GridParity) =
   checkerboard(addr x, y)
 template even*(x: var GridFermion) = x.checkerboard(gpEven)

src/grid/GridUtils.nim

@@ -80,7 +80,7 @@ proc assignStag*(r0: var GridFermion, x: Field) =
       #echo i
       for j in 0..<nd:
         var l = lo.coords[j][i].cint - c0[j]
-        {.emit:"c[j] = l;".}
+        {.emit:["c[j] = ",l,";"].}
       for ic in cint(0)..2:
         var tr,ti: float
         tr := x{i}[ic].re
@@ -120,9 +120,9 @@ proc assignStag*(r0: var Field, x0: var GridFermion) =
     for i in subset.singleSites:
       for j in 0..<nd:
         var l = lo.coords[j][i].cint - c0[j]
-        {.emit:"c[j] = l;".}
+        {.emit:["c[j] = ",l,";"].}
       {.emit:["autoView(dst, ",x[],", CpuRead);"].}
-      {.emit:"peekLocalSite(t, dst, c);".}
+      {.emit:["peekLocalSite(",t,", dst, c);"].}
       for ic in 0..2:
         var tr,ti: float
         {.emit:"tr = t._internal._internal._internal[ic].real();".}

src/grid/gridStag.nim

@@ -2,6 +2,7 @@ import qex
 import physics/stagSolve
 import grid/gridImpl
 
+#[
 template getGrid(g: Field): untyped =
   let lo = g.l
   let latt_size = newCoordinate(lo.physGeom)
@@ -10,8 +11,10 @@ template getGrid(g: Field): untyped =
   let simd_layout = GridDefaultSimd(lo.nDim, Nsimd(GridVComplex));
   let mpi_layout = newCoordinate(lo.rankGeom)
   let grid = newGridCartesian(latt_size,simd_layout,mpi_layout)
+]#
 
-proc gridSolveEE*(s:Staggered; r,t:Field; m:SomeNumber; sp: var SolverParams) =
+proc gridSolveXX*(s:Staggered; r,t:Field; m:SomeNumber; sp: var SolverParams;
+                  parEven = true) =
   let lo = r.l
   let latt_size = newCoordinate(lo.physGeom)
   #let simd_layout = newCoordinate(lo.innerGeom)
@@ -27,46 +30,54 @@ proc gridSolveEE*(s:Staggered; r,t:Field; m:SomeNumber; sp: var SolverParams) =
 
   if s.g.len == 4: # plain staggered
     type ferm = GridNaiveStaggeredFermionR
-    var gsrc0 = rbgrid.fermion(ferm)
-    var gsoln0 = rbgrid.fermion(ferm)
-    gsrc0.even
-    gsoln0.even
-    gsrc0 := t
+    var gsrc = rbgrid.fermion(ferm)
+    var gsoln = rbgrid.fermion(ferm)
+    if parEven:
+      gsrc.even
+      gsoln.even
+    else:
+      gsrc.odd
+      gsoln.odd
+    gsrc := t
     {.emit:"using namespace Grid;".}
-    {.emit:"gsoln0 = Zero();".}
+    {.emit:[gsoln," = Zero();"].}
     s.g.stagPhase([0,1,3,7])
     gfl := s.g[0..3]
     s.g.stagPhase([0,1,3,7])
-    {.emit:"using Stag = NaiveStaggeredFermionR;".}
+    {.emit:["using Stag = ",GridNaiveStaggeredFermionR,";"].}
     {.emit:"using FermionField = Stag::FermionField;".}
-    {.emit:"Stag Ds(grid,rbgrid,2.*mass,2.,1.);".}
+    {.emit:["Stag Ds(grid,rbgrid,2.*",mass,",2.,1.);"].}
     {.emit:"Ds.ImportGauge(gfl);".}
     {.emit:"SchurStaggeredOperator<Stag,FermionField> HermOp(Ds);".}
-    {.emit:"ConjugateGradient<FermionField> CG(res, maxit, false);".}
-    {.emit:"CG(HermOp, gsrc0, gsoln0);".}
-    {.emit:"sp.iterations = CG.IterationsToComplete;".}
+    {.emit:["ConjugateGradient<FermionField> CG(",res,", ",maxit,", false);"].}
+    {.emit:["CG(HermOp, ",gsrc,", ",gsoln,");"].}
+    {.emit:[sp,".iterations = CG.IterationsToComplete;"].}
     var rr = r
-    rr := gsoln0
+    rr := gsoln
   elif s.g.len == 8: # Naik staggered
     type ferm = GridImprovedStaggeredFermionR
     var gsrc = rbgrid.fermion(ferm)
     var gsoln = rbgrid.fermion(ferm)
-    gsrc.even
-    gsoln.even
+    if parEven:
+      gsrc.even
+      gsoln.even
+    else:
+      gsrc.odd
+      gsoln.odd
     gsrc := t
     {.emit:"using namespace Grid;".}
-    {.emit:"gsoln = Zero();".}
+    {.emit:[gsoln," = Zero();"].}
     var gll = grid.gauge()
     gfl := @[s.g[0],s.g[2],s.g[4],s.g[6]]
     gll := @[s.g[1],s.g[3],s.g[5],s.g[7]]
-    {.emit:"using ImpStag = ImprovedStaggeredFermionR;".}
+    {.emit:["using ImpStag = ",GridImprovedStaggeredFermionR,";"].}
     {.emit:"using FermionField = ImpStag::FermionField;".}
     {.emit:"ImpStag Ds(grid,rbgrid,2.*mass,2.,2.,1.);".}
     {.emit:"Ds.ImportGaugeSimple(gll,gfl);".}
     {.emit:"SchurStaggeredOperator<ImpStag,FermionField> HermOp(Ds);".}
     {.emit:"ConjugateGradient<FermionField> CG(res, maxit, false);".}
-    {.emit:"CG(HermOp, gsrc, gsoln);".}
-    {.emit:"sp.iterations = CG.IterationsToComplete;".}
+    {.emit:["CG(HermOp, ",gsrc,", ",gsoln,");"].}
+    {.emit:[sp,".iterations = CG.IterationsToComplete;"].}
     var rr = r
     rr := gsoln
   else:
@@ -80,3 +91,9 @@ proc gridSolveEE*(s:Staggered; r,t:Field; m:SomeNumber; sp: var SolverParams) =
     #soln2 := gsrc
     soln2 := gsoln
   ]#
+
+proc gridSolveEE*(s:Staggered; r,t:Field; m:SomeNumber; sp: var SolverParams) =
+  gridSolveXX(s, r, t, m, sp, parEven=true)
+
+proc gridSolveOO*(s:Staggered; r,t:Field; m:SomeNumber; sp: var SolverParams) =
+  gridSolveXX(s, r, t, m, sp, parEven=false)

src/maths/types.nim

@@ -428,10 +428,16 @@ proc `[]`*(x: Indexed): auto {.alwaysInline,noInit.} =
   #let tIndexedBracket0 = xx
   #let tIndexedBracket0 = getPtr xx; template x:untyped {.gensym.} = tIndexedBracket0[]
   obj(x)[idx(x)]
-template `[]`*(xx:Indexed; i:SomeInteger):untyped =
-  #let tIndexedBracket1 = xx
-  let tIndexedBracket1 = getPtr xx; template x:untyped {.gensym.} = tIndexedBracket1[]
-  optIndexed(obj(x)[i], idx(x))
+#template `[]`*(xx:Indexed; i:SomeInteger):untyped =
+#  #let tIndexedBracket1 = xx
+#  let tIndexedBracket1 = getPtr xx; template x:untyped {.gensym.} = tIndexedBracket1[]
+#  optIndexed(obj(x)[i], idx(x))
+proc `[]`*(x: Indexed; i: SomeInteger, r: var auto) {.alwaysInline.} =
+  r = optIndexed(obj(x)[i], idx(x))
+template `[]`*(x: Indexed; i: SomeInteger): auto =
+  var r {.noInit.}: typeof(optIndexed(obj(x)[i], idx(x)))
+  x[i, r]
+  r
 #template `[]`*(xx:Indexed; i,j:SomeInteger):untyped =
 #  #let tIndexedBracket2 = xx
 #  let tIndexedBracket2 = getPtr xx; template x:auto {.gensym.} = tIndexedBracket2[]

src/physics/stagD.nim

@@ -413,7 +413,7 @@ proc stagDb*(sd:StaggeredD; r:Field; g:openArray[Field2];
   #  rir := m*getVec(x[ir], ic)
 
 # r = m2 - Deo * Doe
-proc stagD2ee*(sde,sdo:StaggeredD; r:Field; g:openArray[Field2];
+proc stagD2xx*(sdx,sdy:StaggeredD; r:Field; g:openArray[Field2];
                x:Field; m2:SomeNumber) =
   tic()
   var t{.global.}:evalType(x)
@@ -424,24 +424,32 @@ proc stagD2ee*(sde,sdo:StaggeredD; r:Field; g:openArray[Field2];
     threadBarrier()
   #threadBarrier()
   #stagD(sdo, t, g, x, 0.0)
-  toc("stagD2ee init")
+  toc("stagD2xx init")
   block:
-    stagDP(sdo, t, g, x, 0):
+    stagDP(sdy, t, g, x, 0):
       rir := 0
-  toc("stagD2ee DP")
+  toc("stagD2xx DP")
   threadBarrier()
-  toc("stagD2ee barrier")
+  toc("stagD2xx barrier")
   #stagD(sde, r, g, t, 0.0)
   block:
-    stagDM(sde, r, g, t, 6):
+    stagDM(sdx, r, g, t, 6):
       rir := (4.0*m2)*x[ir]
-  toc("stagD2ee DM")
+  toc("stagD2xx DM")
   #threadBarrier()
   #r[sde.sub] := m2*x - r
   #for ir in r[sde.subset]:
   #  msubVSVV(r[ir], m2, x[ir], r[ir])
   #r[sde.sub] := 0.25*r
 
+proc stagD2ee*(sde,sdo:StaggeredD; r:Field; g:openArray[Field2];
+               x:Field; m2:SomeNumber) =
+  stagD2xx(sde, sdo, r, g, x, m2)
+
+proc stagD2oo*(sde,sdo:StaggeredD; r:Field; g:openArray[Field2];
+               x:Field; m2:SomeNumber) =
+  stagD2xx(sdo, sde, r, g, x, m2)
+
 # r = m2 + Deo * Doe
 # modified D: Df + Db
 proc stagD2eeFB*(sde,sdo:StaggeredD; r:Field; g:openArray[Field2];