small additions

src/base/stdUtils.nim

@@ -184,6 +184,12 @@ proc `+=`*[T](r: var openArray[T], x: openArray[T]) {.inline.} =
   for i in 0..<n:
     r[i] += x[i]
 
+proc `-`*[T](x: seq[T]): seq[T] {.inline.} =
+  let n = x.len
+  result.newSeq(n)
+  for i in 0..<n:
+    result[i] = -x[i]
+
 #proc sum*[T](x: openArray[T]): T =
 #  for i in 0..<x.len: result += x[i]
 

src/eigens/linalgFuncs.nim

@@ -101,7 +101,7 @@ proc newZmat*(nr,nc: int): zmat =
 template mat_nrows*(x: dmat): int = x.nrows
 template mat_ncols*(x: dmat): int = x.ncols
 template `[]`*(x: dmat|zmat, i,j: int): untyped = x.dat[i+j*x.nrows]
-template `[]=`*(x: dmat|zmat, i,j: int, y: untyped): untyped =
+template `[]=`*(x: dmat|zmat, i,j: int, y: untyped) =
   x.dat[i+j*x.nrows] = y
 template dmat_get*(x: dmat, i,j: int): float64 = x[i,j]
 #template dmat_set(x: dmat, i,j: int, y: float64) = x[i,j] = y

src/field/fieldET.nim

@@ -60,6 +60,8 @@ type
   notSomeField2* = concept x
     x isnot SomeField2
 
+template elemType*(x:Field):typedesc = evalType(x[0])
+template elemType*[V:static[int],T](x:typedesc[Field[V,T]]):typedesc = evalType(T)
 template numberType*(x:Field):untyped = numberType(x[0])
 
 #template fieldUnop*(o: static[FieldOps], x: SomeField): auto =
@@ -218,7 +220,8 @@ template has*[F:Field](x: typedesc[F], y: typedesc): bool =
       has(getT F, y)
     else: false
 
-template `[]`*[F:Field](x:typedesc[F]; i:int):typedesc = F.T
+#template `[]`*[F:Field](x:typedesc[F]; i:int):typedesc = F.T
+template index*[F:Field](x:typedesc[F]; i:typedesc[int]):typedesc = F.T
 template `[]`*(x:Field; i:int):untyped = x.s[i]
 #template `[]=`*(x:Field; i:int; y:typed) =
 proc `[]=`*(x:Field; i:int; y:auto) =

src/maths/complexProxy.nim

@@ -295,13 +295,13 @@ template conjImagU(xr,xi: untyped): untyped =
   -conj(xi)
 unaryOverloads(`~`, conj, conjRealU, conjImagU)
 
-#template adjRealU(xr,xi: untyped): untyped =
-#  mixin adj
-#  adj(xr)
-#template adjImagU(xr,xi: untyped): untyped =
-#  mixin adj
-#  -adj(xi)
-#unaryOverloads(`@`, adj, adjRealU, adjImagU)
+template adjRealU(xr,xi: untyped): untyped =
+  mixin adj
+  adj(xr)
+template adjImagU(xr,xi: untyped): untyped =
+  mixin adj
+  -adj(xi)
+unaryOverloads(`@`, adj, adjRealU, adjImagU)
 
 template transposeRealU(xr,xi: untyped): untyped =
   mixin transpose
@@ -362,7 +362,7 @@ template `/`*(x: ImagProxy): untyped = inv(x)
 #template invComplexU(x: untyped): untyped = x.adj * x.norm2.inv
 #template inv*(x: ComplexProxy): untyped = safecall(invComplexU, x)
 proc inv*(x: ComplexProxy): auto {.alwaysInline,noInit.} =
-  mixin inv
+  mixin adj, inv
   x.adj * x.norm2.inv
 template `/`*(x: ComplexProxy): untyped = inv(x)
 

src/physics/stagD.nim

@@ -51,6 +51,24 @@ template initStagD3T*(l:var Layout; T:typedesc; ss:string):untyped =
 proc initStagD3*(x:Field; sub:string):auto =
   result = initStagD3T(x.l, evalType(x[0]), sub)
 
+proc subdirs*(s: Staggered, dirs: seq[int]): Staggered =
+  result.se.sf.newSeq(0)
+  result.se.sb.newSeq(0)
+  result.se.sub = s.se.sub
+  result.se.subset = s.se.subset
+  result.so.sb.newSeq(0)
+  result.so.sf.newSeq(0)
+  result.so.sub = s.so.sub
+  result.so.subset = s.so.subset
+  result.g.newSeq(0)
+  for i in 0..<dirs.len:
+    let d = dirs[i]
+    result.se.sf.add s.se.sf[d]
+    result.se.sb.add s.se.sb[d]
+    result.so.sf.add s.so.sf[d]
+    result.so.sb.add s.so.sb[d]
+    result.g.add s.g[d]
+
 proc rephase*(s: Staggered) =
   s.g.setBC
   threadBarrier()

src/physics/stagSolve.nim

@@ -11,6 +11,8 @@ import maths
 import quda/qudaWrapper
 import grid/Grid
 
+var precon = false
+
 proc solveEO*(s: Staggered; r,x: Field; m: SomeNumber; sp0: var SolverParams) =
   var sp = sp0
   sp.subset.layoutSubset(r.l, sp.subsetName)
@@ -55,12 +57,16 @@ proc solveXX*(s: Staggered; r,x: Field; m: SomeNumber; sp0: var SolverParams;
         stagD2oo(s.se, s.so, a, s.g, b, m*m)
         toc("stagD2oo")
       #threadBarrier()
-    var oa = (apply: op)
     var cg = newCgState(r, x)
     if parEven:
       sp.subset.layoutSubset(r.l, "even")
     else:
       sp.subset.layoutSubset(r.l, "odd")
+    #if precon:
+      #var oap = (apply: op, applyPrecon: oppre)
+      #cg.solve(oap, sp)
+    #else:
+    var oa = (apply: op)
     cg.solve(oa, sp)
     toc("cg.solve")
     sp.calls = 1
@@ -332,8 +338,15 @@ when isMainModule:
   var r = lo.ColorVector()
   var rs = newRNGField(RngMilc6, lo, intParam("seed", 987654321).uint64)
 
+  if fn == "":
+    var warm = floatParam("warm", 0.15)
+    threads:
+      #g.random rs
+      g.warm warm, rs
+  let plaq = g.plaq
+  echo "plaq: ", plaq.sum, " ", plaq
+
   threads:
-    g.random rs
     g.setBC
     g.stagPhase
     #v1 := 0
@@ -362,6 +375,14 @@ when isMainModule:
     threads:
       r := v1 - v2
       echo "err2: ", r.norm2
+    threads:
+      v1 := 0
+    resetTimers()
+    precon = true
+    s.solve(v1, v3, m, sp)
+    threads:
+      r := v1 - v2
+      echo "err2: ", r.norm2
 
   block:
     v1 := 0
@@ -371,7 +392,7 @@ when isMainModule:
     echo "even point"
     test()
     echo sp.getStats()
-
+#[
   block:
     v1 := 0
     let p = lo.rankIndex([0,0,0,1])
@@ -386,7 +407,7 @@ when isMainModule:
     echo "random"
     test()
     echo sp.getStats()
-
+]#
   if intParam("timers", 0)!=0:
     echoTimers()
 

src/solvers/cg.nim

@@ -6,8 +6,8 @@ export solverBase
 
 type
   CgState*[T] = object
-    r,Ap,b: T
-    p,x,z: T
+    r*,Ap,b*: T
+    p,x*,z: T
     b2,r2,r2old,r2stop,rz,rzold: float
     iterations: int
     precon: bool