Merge branch 'master' into 'master'

Master

See merge request alramos/latticegpu.jl!3

docs/make.jl

@@ -13,6 +13,8 @@ makedocs(sitename="LatticeGPU", modules=[LatticeGPU], doctest=true,
              "Yang-Mills" => "ym.md",
              "Gradient flow" => "flow.md",
              "Schrödinger Functional" => "sf.md",
+             "Dirac" => "dirac.md",
+             "Solvers" => "solvers.md",
              "Input Output" => "io.md"
              ], 
          repo = "https://igit.ific.uv.es/alramos/latticegpu.jl")

docs/src/dirac.md

@@ -0,0 +1,112 @@
+
+# Dirac operator
+
+The module `Dirac` has the necessary stuctures and function 
+to simulate non-dynamical 4-dimensional Wilson fermions.
+
+There are two main data structures in this module, the structure [`DiracParam`](@ref)
+
+```@docs
+DiracParam
+```
+
+and the workspace [`DiracWorkspace`](@ref)
+
+```@docs
+DiracWorkspace
+```
+
+The workspace stores four fermion fields, namely `.sr`, `.sp`, `.sAp` and `.st`, used
+for different purposes. If the representation is either `SU2fund` of `SU3fund`, an extra
+field with values in `U2alg`/`U3alg` is created to store the clover, used for the improvent.
+
+## Functions
+
+The functions [`Dw!`](@ref), [`g5Dw!`](@ref) and [`DwdagDw!`](@ref) are all related to the 
+Wilson-Dirac operator.
+
+The action of the Dirac operator `Dw!` is the following:
+
+```math
+D_w\psi (\vec{x} = x_1,x_2,x_3,x_4) = (4 + m_0 + i \mu \gamma_5)\psi(\vec{x}) - 
+```
+```math
+    - \frac{1}{2}\sum_{\mu = 1}^4 \theta (\mu) (1-\gamma_\mu) U_\mu(\vec{x}) \psi(\vec{x} + \hat{\mu}) + \theta^* (\mu) (1 + \gamma_\mu) U^{-1}_\mu(\vec{x} - \hat{\mu}) \psi(\vec{x} - \hat{\mu})
+```
+
+where $$m_0$$ and $$\theta$$ are respectively the values `.m0` and `.th` of [`DiracParam`](@ref).
+Note that $$|\theta(\mu)|=1$$ is not built into the code, so it should be imposed explicitly. 
+
+Additionally, if |`dpar.csw`| > 1.0E-10, the clover term is assumed to be stored in `ymws.csw`, which
+can be done via the [`Csw!`](@ref) function. In this case we have the Sheikholeslami–Wohlert (SW) term
+in `Dw!`:
+
+```math
+    \delta D_w^{sw} = \frac{i}{2}c_{sw} \sum_{\pi = 1}^6 F^{cl}_\pi \sigma_\pi \psi(\vec{x})
+```
+where the $$\sigma$$ matrices are those described in the `Spinors` module and the index $$\pi$$ runs
+as specified in `lp.plidx`.
+
+If the boudary conditions, defined in `lp`, are either `BC_SF_ORBI,D` or `BC_SF_AFWB`, the 
+improvement term 
+
+```math
+    \delta D_w^{SF} = (c_t -1) (\delta_{x_4,a} \psi(\vec{x}) + \delta_{x_4,T-a} \psi(\vec{x}))
+```
+is added. Since the time-slice $$t=T$$ is not stored, this accounts to modifying the second
+and last time-slice. 
+
+Note that the Dirac operator for SF boundary conditions assumes that the value of the field 
+in the first time-slice is zero. To enforce this, we have the function
+
+```@docs
+SF_bndfix!
+```
+
+Note that this is not enforced in the Dirac operators, so if the field `so` does not satisfy SF
+boundary conditions, it will not (in general) satisfy them after applying [`Dw!`](@ref) 
+or [`g5Dw!`](@ref). This function is called for the function [`DwdagDw!`](@ref), so in this case
+`so` will always be a proper SF field after calling this function.
+
+The function [`Csw!`](@ref) is used to store the clover in `dws.csw`. It is computed 
+according to the expression
+
+```math
+F_{\mu,\nu} = \frac{1}{8} (Q_{\mu \nu} - Q_{\nu \mu})
+```
+
+where 
+```math
+Q_{\mu\nu} = U_\mu(\vec{x})U_{\nu}(x+\mu)U_{\mu}^{-1}(\vec{x}+\nu)U_{\nu}(\vec{x}) +
+U_{\nu}^{-1}(\vec{x}-\nu) U_\mu (\vec{x}-\nu) U_{\nu}(\vec{x} +\mu - \nu) U^{-1}_{\mu}(\vec{x}) +
+```
+```math
++ U^{-1}_{\mu}(x-\mu)U_\nu^{-1}(\vec{x} - \mu - \nu)U_\mu(\vec{x} - \mu - \nu)U_\nu^{-1}(x-\nu) +
+```
+```math
++U_{\nu}(\vec{x})U_{\mu}^{-1}(\vec{x} + \nu - \mu)U^{-1}_{\nu}(\vec{x} - \mu)U_\mu(\vec{x}-\mu)
+
+```
+
+The correspondence between the tensor field and the GPU-Array is the following:
+```math
+F[b,1,r] \to F_{41}(b,r) ,\quad F[b,2,r] \to F_{42}(b,r) ,\quad F[b,3,r] \to F_{43}(b,r) 
+```
+```math
+F[b,4,r] \to F_{31}(b,r) ,\quad F[b,5,r] \to F_{32}(b,r) ,\quad F[b,6,r] \to F_{21}(b,r) 
+```
+where $$(b,r)$$ labels the lattice points as explained in the module `Space`
+
+The function [`pfrandomize!`](@ref), userfull for stochastic sources, is also present. It
+randomizes a  fermion field either in all the space or in a specifit time-slice.
+
+The generic interface of these functions reads
+
+```@docs
+Dw!
+g5Dw!
+DwdagDw!
+Csw!
+pfrandomize!
+mtwmdpar
+```

docs/src/solvers.md

@@ -0,0 +1,86 @@
+
+# Solvers
+
+The module `Solvers` contains the functions to invert the Dirac
+operator as well as functions to obtain specific propagators.
+
+
+## CG.jl
+
+The function [`CG!`](@ref) implements the Conjugate gradient 
+algorith for the operator A
+
+```@docs
+CG!
+```
+
+where the tolerance is normalized with respect to $$|$$``si``$$|^2$$.
+The operator A must have the same input structure as all the Dirac
+operators. If the maximum number of iterations `maxiter` is reached,
+the function will throw an error. The estimation for $$A^{-1}x$$ is
+stored in ``si``, and the number of iterations is returned.
+
+
+Note that all the fermion field in ``dws`` are used
+inside the function and will be modified. In particular, the final residue
+is given by $$|$$``dws.sr``$$|^2$$.
+
+
+
+## Propagators.jl
+
+In this file, we define a couple of useful functions to obtain certain
+propagators.
+
+```@docs
+propagator!
+```
+
+Internally, this function solves the equation 
+
+```math
+    D_w^\dagger D_w \psi = \gamma_5 D_w \gamma_5 \eta
+```
+
+where $$\eta$$ is either a point-source with the specified color and spin
+or a random source in a time-slice and stores the value in ``pro``. 
+To solve this equation, the [`CG!`](@ref) function is used.
+
+
+For the case of SF boundary conditions, we have the boundary-to-bulk 
+propagator, implemented by the function [`bndpropagator!`](@ref)
+
+```@docs
+bndpropagator!
+```
+
+This propagator is defined by the equation:
+
+```math
+D_W S(x) = \frac{c_t}{\sqrt{V}} \delta_{x_0,1} U_0^\dagger(0,\vec{x}) P_+
+```
+
+The analog for the other boundary is implemented in the function [`Tbndpropagator!`](@ref)
+
+```@docs
+Tbndpropagator!
+```
+
+defined by the equation:
+
+```math
+D_W R(x) = \frac{c_t}{\sqrt{V}} \delta_{x_0,T-1} U_0(T-1,\vec{x}) P_-
+```
+
+Where $$P_\pm = (1 \pm \gamma_0)/2$$. The boundary-to-boundary 
+propagator
+
+```math
+    - \frac{c_t}{\sqrt{V}} \sum_\vec{x} U_0 ^\dagger (T-1,\vec{x}) P_+ S(T-1,\vec{x})
+```
+
+is computed by the function [`bndtobnd`](@ref)
+
+```@docs
+bndtobnd
+```

docs/src/spinors.md

@@ -0,0 +1,86 @@
+# Spinors
+
+The module Spinors defines the necessary functions for the structure `Spinor{NS,G}`,
+which is a NS-tuple with values in G.
+
+The functions `norm`, `norm2`, `dot`, `*`, `/`, `/`, `+`, `-`, `imm` and `mimm`, 
+if defined for G, are extended to Spinor{NS,G} for general NS.
+
+For the 4d case where NS = 4 there are some specific functions to implement different
+operations with the gamma matrices. The convention for these matrices is
+
+
+```math
+\gamma _4 = \left( 
+    \begin{array}{cccc}
+        0 & 0 & -1 & 0\\
+        0 & 0 & 0 & -1\\
+        -1 & 0 & 0 & 0\\
+        0 & -1 & 0 & 0\\
+    \end{array}
+    \right)
+    \quad 
+    \gamma_1 = \left( 
+    \begin{array}{cccc}
+        0 & 0 & 0 & -i\\
+        0 & 0 & -i & 0\\
+        0 & i & 0 & 0\\
+        i & 0 & 0 & 0\\
+    \end{array}
+    \right)
+```
+```math 
+    \gamma _2 = \left( 
+    \begin{array}{cccc}
+        0 & 0 & 0 & -1\\
+        0 & 0 & 1 & 0\\
+        0 & 1 & 0 & 0\\
+        -1 & 0 & 0 & 0\\
+    \end{array}
+    \right)
+    \quad 
+    \gamma_3 = \left( 
+    \begin{array}{cccc}
+        0 & 0 & -i & 0\\
+        0 & 0 & 0 & i\\
+        i & 0 & 0 & 0\\
+        0 & -i & 0 & 0\\
+    \end{array}
+    \right)
+```
+
+
+The function [`dmul`](@ref) implements the multiplication over the $$\gamma$$ matrices
+
+```@docs
+dmul
+```
+
+The function [`pmul`](@ref) implements the $$ (1 \pm \gamma_N) $$ proyectors. The functions 
+[`gpmul`](@ref) and [`gdagpmul`](@ref) do the same and then multiply each element by `g`and 
+g^-1 repectively.
+
+```@docs
+pmul
+gpmul
+gdagpmul
+```
+
+## Some examples
+
+Here we just display some examples for these functions. We display it with `ComplexF64`
+instead of `SU3fund` or `SU2fund` for simplicity.
+
+
+```@setup exs
+import Pkg # hide
+Pkg.activate("/home/alberto/code/julia/LatticeGPU/") # hide
+using LatticeGPU # hide
+```
+```@repl exs
+spin = Spinor{4,Complex{Float64}}((1.0,im*0.5,2.3,0.0))
+println(spin)
+println(dmul(Gamma{4},spin))
+println(pmul(Pgamma{2,-1},spin))
+
+```

src/Dirac/Dirac.jl

@@ -19,20 +19,31 @@ using ..Fields
 using ..YM
 using ..Spinors
 
+"""
+    struct DiracParam{T,R}
 
+Stores the parameters of the Dirac operator. It can be generated via the constructor `function DiracParam{T}(::Type{R},m0,csw,th,tm,ct)`. The first argument can be ommited and is taken to be `SU3fund`.
+The parameters are:
+
+- `m0::T` : Mass of the fermion
+- `csw::T` : Improvement coefficient for the Csw term
+- `th{Ntuple{4,Complex{T}}}` : Phase for the fermions included in the boundary conditions, reabsorbed in the Dirac operator.
+- `tm` : Twisted mass parameter
+- `ct` : Boundary improvement term, only used for Schrödinger Funtional boundary conditions. 
+"""
 struct DiracParam{T,R}
     m0::T
     csw::T
     th::NTuple{4,Complex{T}}
+    tm::T
     ct::T
 
-
+    function DiracParam{T}(::Type{R},m0,csw,th,tm,ct) where {T,R}
-    function DiracParam{T}(::Type{R},m0,csw,th,ct) where {T,R} 
+        return new{T,R}(m0,csw,th,tm,ct)
-        return new{T,R}(m0,csw,th,ct)
     end
 
-    function DiracParam{T}(m0,csw,th,ct) where {T} 
+    function DiracParam{T}(m0,csw,th,tm,ct) where {T}
-        return new{T,SU3fund}(m0,csw,th,ct)
+        return new{T,SU3fund}(m0,csw,th,tm,ct)
     end
 end
 function Base.show(io::IO, dpar::DiracParam{T,R}) where {T,R}
@@ -40,11 +51,24 @@ function Base.show(io::IO, dpar::DiracParam{T,R}) where {T,R}
     println(io, "Wilson fermions in the: ", R, " representation")
     println(io, " - Bare mass: ", dpar.m0," // Kappa = ",0.5/(dpar.m0+4))
     println(io, " - Csw :                ", dpar.csw)
-    println(io, " - c_t:                 ", dpar.ct)
     println(io, " - Theta:               ", dpar.th)
+    println(io, " - Twisted mass:        ", dpar.tm)
+    println(io, " - c_t:                 ", dpar.ct)
     return nothing
 end
 
+
+"""
+    struct DiracWorkspace{T}    
+
+Workspace needed to work with fermion fields. It contains four scalar fermion fields and, for the SU2fund and SU3fund, a U(N) field to store the clover term.
+
+It can be created with the constructor `DiracWorkspace(::Type{G}, ::Type{T}, lp::SpaceParm{4,6,B,D})`. For example:
+
+dws = DiracWorkspace(SU2fund,Float64,lp);
+dws = DiracWorkspace(SU3fund,Float64,lp);
+
+"""
 struct DiracWorkspace{T}
     sr
     sp
@@ -146,18 +170,26 @@ function krnl_csw!(csw::AbstractArray{T}, U, Ubnd, ipl, lp::SpaceParm{4,M,B,D})
     return nothing
 end
 
+
+"""
+    function Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D})
+
+Computes the Dirac operator (with the Wilson term) `\`\``D_w``\`\` with gauge field U and parameters `dpar` of the field `si` and stores it in `so`. 
+If `dpar.csw` is different from zero, the clover term should be stored in `dws.csw` via the Csw! function and is automatically included in the operator.
+
+"""
 function Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D}) where {B,D}
        
     if abs(dpar.csw) > 1.0E-10
         @timeit "Dw" begin
             CUDA.@sync begin
-               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, lp)
+               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, lp)
             end
         end
     else
         @timeit "Dw" begin
             CUDA.@sync begin
-               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.th, lp)
+               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.tm, dpar.th, lp)
             end
         end
     end
@@ -165,7 +197,7 @@ function Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6
     return nothing
 end
 
-function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, lp::SpaceParm{4,6,B,D}) where {B,D}
+function krnl_Dwimpr!(so, U, si, Fcsw, m0, tm, th, csw, lp::SpaceParm{4,6,B,D}) where {B,D}
 
     b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
 
@@ -180,7 +212,7 @@ function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, lp::SpaceParm{4,6,B,D}) wher
 
     @inbounds begin 
         
-        so[b,r] = (4+m0)*si[b,r]  + 0.5*csw*im*( Fcsw[b,1,r]*dmul(Gamma{10},si[b,r]) + Fcsw[b,2,r]*dmul(Gamma{11},si[b,r]) + Fcsw[b,3,r]*dmul(Gamma{12},si[b,r])
+        so[b,r] = (4+m0)*si[b,r]+ im*tm*dmul(Gamma{5},si[b,r]) + 0.5*csw*im*( Fcsw[b,1,r]*dmul(Gamma{10},si[b,r]) + Fcsw[b,2,r]*dmul(Gamma{11},si[b,r]) + Fcsw[b,3,r]*dmul(Gamma{12},si[b,r])
                                                                              +Fcsw[b,4,r]*dmul(Gamma{15},si[b,r]) + Fcsw[b,5,r]*dmul(Gamma{14},si[b,r]) + Fcsw[b,6,r]*dmul(Gamma{13},si[b,r]))
 
         so[b,r] -= 0.5*(th[1]*gpmul(Pgamma{1,-1},U[b,1,r],si[bu1,ru1]) +conj(th[1])*gdagpmul(Pgamma{1,+1},U[bd1,1,rd1],si[bd1,rd1]) +
@@ -193,7 +225,7 @@ function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, lp::SpaceParm{4,6,B,D}) wher
     return nothing
 end
 
-function krnl_Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
+function krnl_Dw!(so, U, si, m0, tm, th, lp::SpaceParm{4,6,B,D}) where {B,D}
 
     b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
 
@@ -208,7 +240,7 @@ function krnl_Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
 
     @inbounds begin 
         
-        so[b,r] = (4+m0)*si[b,r]
+        so[b,r] = (4+m0)*si[b,r] + im*tm*dmul(Gamma{5},si[b,r])
 
         so[b,r] -= 0.5*(th[1]*gpmul(Pgamma{1,-1},U[b,1,r],si[bu1,ru1]) +conj(th[1])*gdagpmul(Pgamma{1,+1},U[bd1,1,rd1],si[bd1,rd1]) +
                         th[2]*gpmul(Pgamma{2,-1},U[b,2,r],si[bu2,ru2]) +conj(th[2])*gdagpmul(Pgamma{2,+1},U[bd2,2,rd2],si[bd2,rd2]) +
@@ -225,13 +257,13 @@ function Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{SpacePa
     if abs(dpar.csw) > 1.0E-10
         @timeit "Dw" begin
             CUDA.@sync begin
-               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, dpar.ct, lp)
+               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, dpar.ct, lp)
             end
         end
     else
         @timeit "Dw" begin
             CUDA.@sync begin
-               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.th, dpar.ct, lp)
+               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.tm, dpar.th, dpar.ct, lp)
             end
         end
     end
@@ -239,7 +271,7 @@ function Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{SpacePa
     return nothing
 end
 
-function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
+function krnl_Dwimpr!(so, U, si, Fcsw, m0, tm, th, csw, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
 
     # The field si is assumed to be zero at t = 0
 
@@ -258,7 +290,7 @@ function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, ct, lp::Union{SpaceParm{4,6,
 
         @inbounds begin 
             
-            so[b,r] = (4+m0)*si[b,r]  + 0.5*csw*im*( Fcsw[b,1,r]*dmul(Gamma{10},si[b,r]) + Fcsw[b,2,r]*dmul(Gamma{11},si[b,r]) + Fcsw[b,3,r]*dmul(Gamma{12},si[b,r])
+            so[b,r] = (4+m0)*si[b,r]  + im*tm*dmul(Gamma{5},si[b,r]) + 0.5*csw*im*( Fcsw[b,1,r]*dmul(Gamma{10},si[b,r]) + Fcsw[b,2,r]*dmul(Gamma{11},si[b,r]) + Fcsw[b,3,r]*dmul(Gamma{12},si[b,r])
                                                                                    +Fcsw[b,4,r]*dmul(Gamma{15},si[b,r]) + Fcsw[b,5,r]*dmul(Gamma{14},si[b,r]) + Fcsw[b,6,r]*dmul(Gamma{13},si[b,r]))
 
             
@@ -276,7 +308,7 @@ function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, ct, lp::Union{SpaceParm{4,6,
     return nothing
 end
 
-function krnl_Dw!(so, U, si, m0, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
+function krnl_Dw!(so, U, si, m0, tm, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
 
     # The field si is assumed to be zero at t = 0
 
@@ -295,7 +327,7 @@ function krnl_Dw!(so, U, si, m0, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},S
 
         @inbounds begin 
             
-            so[b,r] = (4+m0)*si[b,r]
+            so[b,r] = (4+m0)*si[b,r] + im*tm*dmul(Gamma{5},si[b,r])
             so[b,r] -= 0.5*(th[1]*gpmul(Pgamma{1,-1},U[b,1,r],si[bu1,ru1]) +conj(th[1])*gdagpmul(Pgamma{1,+1},U[bd1,1,rd1],si[bd1,rd1]) +
                             th[2]*gpmul(Pgamma{2,-1},U[b,2,r],si[bu2,ru2]) +conj(th[2])*gdagpmul(Pgamma{2,+1},U[bd2,2,rd2],si[bd2,rd2]) +
                             th[3]*gpmul(Pgamma{3,-1},U[b,3,r],si[bu3,ru3]) +conj(th[3])*gdagpmul(Pgamma{3,+1},U[bd3,3,rd3],si[bd3,rd3]) +
@@ -310,18 +342,24 @@ function krnl_Dw!(so, U, si, m0, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},S
     return nothing
 end
 
+"""
+    function g5Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D})    
+
+Computes \`\` \\gamma_5 \`\` times the Dirac operator (with the Wilson term) with gauge field U and parameters `dpar` of the field `si` and stores it in `so`. 
+If `dpar.csw` is different from zero, the clover term should be stored in `dws.csw` via the Csw! function and is automatically included in the operator.
+"""
 function g5Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D}) where {B,D}
        
     if abs(dpar.csw) > 1.0E-10
         @timeit "g5Dw" begin
             CUDA.@sync begin
-               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, lp)
+               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, lp)
             end
         end
     else
         @timeit "g5Dw" begin
             CUDA.@sync begin
-               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.th, lp)
+               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.tm, dpar.th, lp)
             end
         end
     end
@@ -329,7 +367,7 @@ function g5Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4
     return nothing
 end
 
-function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, th, csw, lp::SpaceParm{4,6,B,D}) where {B,D}
+function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, tm, th, csw, lp::SpaceParm{4,6,B,D}) where {B,D}
 
     b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
 
@@ -352,13 +390,13 @@ function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, th, csw, lp::SpaceParm{4,6,B,D}) wh
                         th[3]*gpmul(Pgamma{3,-1},U[b,3,r],si[bu3,ru3]) +conj(th[3])*gdagpmul(Pgamma{3,+1},U[bd3,3,rd3],si[bd3,rd3]) +
                         th[4]*gpmul(Pgamma{4,-1},U[b,4,r],si[bu4,ru4]) +conj(th[4])*gdagpmul(Pgamma{4,+1},U[bd4,4,rd4],si[bd4,rd4]) )
 
-        so[b,r] = dmul(Gamma{5}, so[b,r])
+        so[b,r] = dmul(Gamma{5}, so[b,r])+ im*tm*si[b,r]
     end
 
     return nothing
 end
 
-function krnl_g5Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
+function krnl_g5Dw!(so, U, si, m0, tm, th, lp::SpaceParm{4,6,B,D}) where {B,D}
 
     b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
 
@@ -380,7 +418,7 @@ function krnl_g5Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
                         th[3]*gpmul(Pgamma{3,-1},U[b,3,r],si[bu3,ru3]) +conj(th[3])*gdagpmul(Pgamma{3,+1},U[bd3,3,rd3],si[bd3,rd3]) +
                         th[4]*gpmul(Pgamma{4,-1},U[b,4,r],si[bu4,ru4]) +conj(th[4])*gdagpmul(Pgamma{4,+1},U[bd4,4,rd4],si[bd4,rd4]) )
 
-        so[b,r] = dmul(Gamma{5}, so[b,r])
+        so[b,r] = dmul(Gamma{5}, so[b,r]) + im*tm*si[b,r]
     end
 
     return nothing
@@ -391,13 +429,13 @@ function g5Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{Space
     if abs(dpar.csw) > 1.0E-10
         @timeit "g5Dw" begin
             CUDA.@sync begin
-               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, dpar.ct, lp)
+               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, dpar.ct, lp)
             end
         end
     else
         @timeit "g5Dw" begin
             CUDA.@sync begin
-               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.th, dpar.ct, lp)
+               CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.tm, dpar.th, dpar.ct, lp)
             end
         end
     end
@@ -405,7 +443,7 @@ function g5Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{Space
     return nothing
 end
 
-function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, th, csw, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
+function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, tm, th, csw, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
 
     # The field si is assumed to be zero at t = 0
 
@@ -439,12 +477,12 @@ function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, th, csw, ct, lp::Union{SpaceParm{4,
         end
     end
 
-    so[b,r] = dmul(Gamma{5}, so[b,r])
+    so[b,r] = dmul(Gamma{5}, so[b,r])+ im*tm*si[b,r]
 
     return nothing
 end
 
-function krnl_g5Dw!(so, U, si, m0, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
+function krnl_g5Dw!(so, U, si, m0, tm, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
 
     # The field si is assumed to be zero at t = 0
 
@@ -475,11 +513,17 @@ function krnl_g5Dw!(so, U, si, m0, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D}
         end
     end
 
-    so[b,r] = dmul(Gamma{5}, so[b,r])
+    so[b,r] = dmul(Gamma{5}, so[b,r]) + im*tm*si[b,r]
 
     return nothing
 end
 
+"""
+    function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D})
+
+Applies the operator \`\` \\gamma_5 D_w \`\` twice to `si` and stores the result in `so`. This is equivalent to appling the operator \`\` D_w^\\dagger D_w \`\`
+The Dirac operator is the same as in the functions `Dw!` and `g5Dw!`
+"""
 function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
 
     if abs(dpar.csw) > 1.0E-10
@@ -487,30 +531,32 @@ function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{Sp
             
             @timeit "g5Dw" begin
                 CUDA.@sync begin
-                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(dws.st, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, dpar.ct, lp)
+                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(dws.st, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, dpar.ct, lp)
                 end
             end
-
+            SF_bndfix!(dws.st,lp)
             @timeit "g5Dw" begin
                 CUDA.@sync begin
-                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, dws.st, dws.csw, dpar.m0, dpar.th, dpar.csw, dpar.ct, lp)
+                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, dws.st, dws.csw, dpar.m0, -dpar.tm, dpar.th, dpar.csw, dpar.ct, lp)
                 end
             end
+            SF_bndfix!(so,lp)
         end
     else
         @timeit "DwdagDw" begin
             
             @timeit "g5Dw" begin
                 CUDA.@sync begin
-                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(dws.st, U, si, dpar.m0, dpar.th, dpar.ct, lp)
+                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(dws.st, U, si, dpar.m0, dpar.tm, dpar.th, dpar.ct, lp)
                 end
             end
-
+            SF_bndfix!(dws.st,lp)
             @timeit "g5Dw" begin
                 CUDA.@sync begin
-                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, dws.st, dpar.m0, dpar.th, dpar.ct, lp)
+                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, dws.st, dpar.m0, -dpar.tm, dpar.th, dpar.ct, lp)
                 end
             end
+            SF_bndfix!(so,lp)
         end
     end
 
@@ -524,13 +570,13 @@ function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpacePar
             
             @timeit "g5Dw" begin
                 CUDA.@sync begin
-                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(dws.st, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, lp)
+                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(dws.st, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, lp)
                 end
             end
 
             @timeit "g5Dw" begin
                 CUDA.@sync begin
-                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, dws.st, dws.csw, dpar.m0, dpar.th, dpar.csw, lp)
+                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, dws.st, dws.csw, dpar.m0, -dpar.tm, dpar.th, dpar.csw, lp)
                 end
             end
         end
@@ -539,13 +585,13 @@ function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpacePar
             
             @timeit "g5Dw" begin
                 CUDA.@sync begin
-                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(dws.st, U, si, dpar.m0, dpar.th, lp)
+                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(dws.st, U, si, dpar.m0, dpar.tm, dpar.th, lp)
                 end
             end
 
             @timeit "g5Dw" begin
                 CUDA.@sync begin
-                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, dws.st, dpar.m0, dpar.th, lp)
+                    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, dws.st, dpar.m0, -dpar.tm, dpar.th, lp)
                 end
             end
         end
@@ -554,12 +600,27 @@ function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpacePar
     return nothing
 end
 
+"""
+    function mtwmdpar(dpar::DiracParam)
 
+Returns `dpar` with oposite value of the twisted mass.
+"""
+function mtwmdpar(dpar::DiracParam{P,R}) where {P,R}
+    return DiracParam{P}(R,dpar.m0,dpar.csw,dpar.th,-dpar.tm,dpar.ct)
+end
+
+
+"""
+    SF_bndfix!(sp, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}})
+
+Sets all the values of `sp` in the  first time slice to zero.
+"""
 function SF_bndfix!(sp, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
+    @timeit "SF boundary fix" begin
         CUDA.@sync begin
             CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfbndfix!(sp, lp)
         end
-    
+    end
     return nothing
 end
 
@@ -643,7 +704,7 @@ function krnl_assign_pf_su2!(f::AbstractArray, p , lp::SpaceParm, t::Int64)
     return nothing
 end
 
-export Dw!, g5Dw!, DwdagDw!, SF_bndfix!, Csw!, pfrandomize!
+export Dw!, g5Dw!, DwdagDw!, SF_bndfix!, Csw!, pfrandomize!, mtwmdpar
 
 
 include("DiracIO.jl")

src/Groups/AlgebraU2.jl

@@ -1,12 +1,23 @@
 
 
 
+"""
+    struct U2alg{T} <: Algebra
+
+Elements of the `U(2)` Algebra. The type `T <: AbstractFloat` can be used to define single or double precision elements.
+"""
 struct U2alg{T} <: Algebra
     u11::T
     u22::T
     u12::Complex{T}
 end
 
+
+"""
+    antsym(a::SU2{T}) where T <: AbstractFloat
+
+Returns the antisymmetrization of the SU2 element `a`, that is `\`\` `a - a^{\\dagger}` `\`. This method returns al element of `U2alg{T}`.
+"""
 function antsym(a::SU2{T}) where T <: AbstractFloat
     return U2alg{T}(2.0*imag(a.t1),-2.0*imag(a.t1),2.0*a.t2)
 end

src/Groups/AlgebraU3.jl

@@ -1,6 +1,10 @@
 
 
+"""
+    struct U3alg{T} <: Algebra
 
+Elements of the `U(3)` Algebra. The type `T <: AbstractFloat` can be used to define single or double precision elements.
+"""
 struct U3alg{T} <: Algebra
     u11::T
     u22::T
@@ -10,6 +14,11 @@ struct U3alg{T} <: Algebra
     u23::Complex{T}
 end
 
+"""
+    antsym(a::SU3{T}) where T <: AbstractFloat
+
+Returns the antisymmetrization of the SU3 element `a`, that is `\`\` `a - a^{\\dagger}` `\`. This method returns al element of `U3alg{T}`.
+"""
 function antsym(a::SU3{T}) where T <: AbstractFloat
     t1 = 2.0*imag(a.u11)
     t2 = 2.0*imag(a.u22)

src/LatticeGPU.jl

@@ -46,19 +46,19 @@ export Eoft_clover, Eoft_plaq, Qtop
 export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3
 export flw, flw_adapt
 export sfcoupling, bndfield, setbndfield
-export import_lex64, import_cern64, import_bsfqcd, save_cnfg, read_cnfg
+export import_lex64, import_cern64, import_bsfqcd, save_cnfg, read_cnfg, read_gp
 
 include("Spinors/Spinors.jl")
 
 using .Spinors    
-export Spinor, Pgamma
+export Spinor, Pgamma, Gamma
 export imm, mimm
 export pmul, gpmul, gdagpmul, dmul
 
 include("Dirac/Dirac.jl")
 using .Dirac
 export DiracWorkspace, DiracParam
-export Dw!, g5Dw!, DwdagDw!, SF_bndfix!, Csw!, pfrandomize!
+export Dw!, g5Dw!, DwdagDw!, SF_bndfix!, Csw!, pfrandomize!, mtwmdpar
 export read_prop, save_prop, read_dpar
 
 include("Solvers/Solvers.jl")

src/Solvers/CG.jl

@@ -9,11 +9,6 @@
 ### created: Tue Nov 30 11:10:57 2021
 ###                               
 
-"""
-    function CG!
-
-Solves the linear equation `Ax = si`
-"""
 function krnl_dot!(sum,fone,ftwo)
     b=Int64(CUDA.threadIdx().x)
     r=Int64(CUDA.blockIdx().x)
@@ -23,7 +18,7 @@ function krnl_dot!(sum,fone,ftwo)
 return nothing
 end
 
-function field_dot(fone::AbstractArray,ftwo::AbstractArray,sumf,lp) where {T}
+function field_dot(fone::AbstractArray,ftwo::AbstractArray,sumf,lp)
         
     CUDA.@sync begin
             CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_dot!(sumf,fone,ftwo)
@@ -32,6 +27,12 @@ function field_dot(fone::AbstractArray,ftwo::AbstractArray,sumf,lp) where {T}
     return sum(sumf)
 end
 
+
+"""
+    function CG!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T}, maxiter::Int64 = 10, tol=1.0)
+
+Solves the linear equation `Ax = si`
+"""
 function CG!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T}, maxiter::Int64 = 10, tol=1.0) where {T}
 
     dws.sr .= si

src/Solvers/Propagators.jl

@@ -25,7 +25,7 @@ function propagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Space
         CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp)
     end
        
-    g5Dw!(pro,U,dws.sp,dpar,dws,lp)
+    g5Dw!(pro,U,dws.sp,mtwmdpar(dpar),dws,lp)
       
     CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
     return nothing
@@ -46,11 +46,7 @@ function propagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Space
         CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp)
     end
        
-    g5Dw!(pro,U,dws.sp,dpar,dws,lp)
+    g5Dw!(pro,U,dws.sp,mtwmdpar(dpar),dws,lp)
-
-    CUDA.@sync begin
-        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp)
-    end
       
     CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
     return nothing
@@ -60,7 +56,7 @@ end
 
     function bndpropagator!(pro,U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64)
         
-Saves the propagator in from the t=0 boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's'. The factor c_t is included while the factor 1/sqrt(V) is not.
+Saves the propagator from the t=0 boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's' in 'pro'. The factor c_t is included while the factor 1/sqrt(V) is not.
 For the propagator from T to the bulk, use the function Tbndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64)
 
 """
@@ -85,6 +81,7 @@ function bndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Sp
         return nothing
     end
 
+    SF_bndfix!(pro,lp)
     fill!(dws.sp,zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
     
     CUDA.@sync begin
@@ -95,15 +92,15 @@ function bndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Sp
         CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp)
     end
        
-    g5Dw!(pro,U,dpar.ct*dws.sp,dpar,dws,lp)
+    g5Dw!(pro,U,dpar.ct*dws.sp,mtwmdpar(dpar),dws,lp)
     
     CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
-    return pro
+    return nothing
 end
 
 """
 
-    function Tbndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64)
+    function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64)
         
 Returns the propagator from the t=T boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's'. The factor c_t is included while the factor 1/sqrt(V) is not.
 For the propagator from t=0 to the bulk, use the function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64)
@@ -129,6 +126,7 @@ function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::S
         return nothing
     end
 
+    SF_bndfix!(pro,lp)
     fill!(dws.sp,zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
     
     CUDA.@sync begin
@@ -139,10 +137,10 @@ function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::S
             CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp)
     end
        
-    g5Dw!(pro,U,dpar.ct*dws.sp,dpar,dws,lp)
+    g5Dw!(pro,U,dpar.ct*dws.sp,mtwmdpar(dpar),dws,lp)
     
     CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
-    return pro
+    return nothing
 end
 
 

src/Spinors/Spinors.jl

@@ -169,7 +169,7 @@ end
 
 
 """
-    gpmul(pgamma{N,S}, g::G, a::Spinor) G <: Group
+    gpmul(Pgamma{N,S}, g::G, a::Spinor) G <: Group
 
 Returns ``g(1+s\\gamma_N)a``
 """
@@ -226,7 +226,7 @@ end
 end
 
 """
-    gdagpmul(pgamma{N,S}, g::G, a::Spinor) G <: Group
+    gdagpmul(Pgamma{N,S}, g::G, a::Spinor) G <: Group
 
 Returns ``g^+ (1+s\\gamma_N)a``
 """
@@ -284,33 +284,33 @@ end
 
 
 # dummy structs for dispatch:
-# Basis of \\Gamma_n
+# Basis of \\gamma_n
 struct Gamma{N}
 end
 
 """
-    dmul(n::Int64, a::Spinor)
+    dmul(Gamma{n}, a::Spinor)
 
-Returns ``\\Gamma_n a``
+Returns ``\\gamma_n a``
 
-indexing for Dirac basis ``\\Gamma_n``:
+indexing for Dirac basis ``\\gamma_n``:
 
- 1  gamma1
+ 1  gamma1;
- 2  gamma2
+ 2  gamma2;
- 3  gamma3
+ 3  gamma3;
- 4  gamma0
+ 4  gamma0;
- 5  gamma5
+ 5  gamma5;
- 6  gamma1 gamma5
+ 6  gamma1 gamma5;
- 7  gamma2 gamma5
+ 7  gamma2 gamma5;
- 8  gamma3 gamma5
+ 8  gamma3 gamma5;
- 9  gamma0 gamma5
+ 9  gamma0 gamma5;
-10  sigma01
+10  sigma01;
-11  sigma02
+11  sigma02;
-12  sigma03
+12  sigma03;
-13  sigma21
+13  sigma21;
-14  sigma32
+14  sigma32;
-15  sigma31
+15  sigma31;
-16  identity
+16  identity;
 
 """
 @inline dmul(::Type{Gamma{1}}, a::Spinor{NS,G}) where {NS,G} = Spinor{NS,G}((mimm(a.s[4]), mimm(a.s[3]), imm(a.s[2]), imm(a.s[1])))

src/YM/YM.jl

@@ -179,6 +179,6 @@ include("YMsf.jl")
 export sfcoupling, bndfield, setbndfield
 
 include("YMio.jl")
-export import_lex64, import_cern64, import_bsfqcd, save_cnfg, read_cnfg
+export import_lex64, import_cern64, import_bsfqcd, save_cnfg, read_cnfg, read_gp
 
 end

src/YM/YMio.jl

@@ -297,3 +297,50 @@ function import_cern64(fname, ibc, lp::SpaceParm; log=true)
 
     return CuArray(Ucpu)
 end
+
+
+
+"""
+    read_gp(fname::String)
+
+Reads Gauge parameters from file `fname` using the native (BDIO) format. Returns GaugeParm and SpaceParm.
+"""
+function read_gp(fname::String)
+
+    UID_HDR = 14
+    fb = BDIO_open(fname, "r")
+    while BDIO_get_uinfo(fb) != UID_HDR
+        BDIO_seek!(fb)
+    end
+    ihdr = Vector{Int32}(undef, 2)
+    BDIO_read(fb, ihdr)
+    if (ihdr[1] != convert(Int32, 1653996111)) && (ihdr[2] != convert(Int32, 2))
+        error("Wrong file format [header]")
+    end
+
+    run = BDIO.BDIO_read_str(fb)
+
+    while BDIO_get_uinfo(fb) != 1
+        BDIO_seek!(fb)
+    end
+
+    ifoo = Vector{Int32}(undef, 4)
+    BDIO_read(fb, ifoo)
+    ndim = convert(Int64, ifoo[1])
+    npls = convert(Int64, round(ndim*(ndim-1)/2))
+    ibc  = convert(Int64, ifoo[2])
+    nf   = ifoo[4]
+    
+    ifoo = Vector{Int32}(undef, ndim+convert(Int32, npls))
+    BDIO_read(fb, ifoo)
+    iL  = ntuple(i -> convert(Int64, ifoo[i]),ndim)
+    ntw = ntuple(i -> convert(Int64, ifoo[i+ndim]), npls)
+
+    dfoo = Vector{Float64}(undef, 4)
+    BDIO_read(fb, dfoo)
+
+    lp = SpaceParm{ndim}(iL, (4,4,4,4), ibc, ntw)
+    gp = GaugeParm{Float64}(SU3{Float64}, dfoo[1], dfoo[2])
+    
+    return gp, lp
+end

test/dirac/test_fp_fa.jl

@@ -8,20 +8,20 @@ using TimerOutputs
 
 function fP_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1.0e-16)
 
-@timeit "fP inversion (x12)" begin
+    @timeit "fP inversion (x12)" begin
 
-lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0));
+    lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0));
-exptheta = exp.(im.*theta./lp.iL);
+    exptheta = exp.(im.*theta./lp.iL);
 
-dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,1.0);
+    dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,0.0,1.0);
-dws = DiracWorkspace(SU3fund,Float64,lp);
+    dws = DiracWorkspace(SU3fund,Float64,lp);
 
-U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64}));
+    U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64}));
-psi = scalar_field(Spinor{4,SU3fund{Float64}},lp);
+    psi = scalar_field(Spinor{4,SU3fund{Float64}},lp);
 
-res = zeros(lp.iL[4])
+    res = zeros(lp.iL[4])
 
-for s in 1:4 for c in 1:3
+    for s in 1:4 for c in 1:3
         bndpropagator!(psi,U,dpar,dws,lp,1000,prec,c,s);
 
         for t in 1:lp.iL[4]
@@ -33,11 +33,11 @@ for s in 1:4 for c in 1:3
 
         end
 
-end end
+    end end
 
-end
+    end
 
-@timeit "fP analitical solution" begin
+    @timeit "fP analitical solution" begin
 
         #THEORETICAL SOLUTION: hep-lat/9606016 eq (2.33)
 
@@ -53,14 +53,39 @@ end
             res_th[i] = (2*3*sinh(omega)/(Rpp^2)) * ( (Mpp + sinh(omega))*exp(-2*omega*(i-1)) - (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1))) )
         end
 
-end
+    end
     return sum(abs.(res-res_th))
 
 end
 
 function fA_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1.0e-16)
 
-@timeit "fA inversion (x12)" begin
+    @timeit "fA inversion (x12)" begin
+
+        lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0));
+        exptheta = exp.(im.*theta./lp.iL);
+
+        dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,0.0,1.0);
+        dws = DiracWorkspace(SU3fund,Float64,lp);
+
+        U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64}));
+        psi = scalar_field(Spinor{4,SU3fund{Float64}},lp);
+
+        res = im*zeros(lp.iL[4])
+
+        for s in 1:4 for c in 1:3
+            bndpropagator!(psi,U,dpar,dws,lp,1000,prec,c,s);
+
+            for t in 1:lp.iL[4]
+            #for i in 1:lp.iL[1]    for j in 1:lp.iL[2]        for k in 1:lp.iL[3]
+                        i=abs(rand(Int))%lp.iL[1] +1;j=abs(rand(Int))%lp.iL[2] +1;k=abs(rand(Int))%lp.iL[3] +1;
+                        CUDA.@allowscalar (res[t] += -dot(psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...],dmul(Gamma{4},psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...]))/2)
+            #end end end
+            #res[t] = res[t]/(lp.iL[1]*lp.iL[2]*lp.iL[3])
+
+            end
+
+        end end
 
     lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0));
     exptheta = exp.(im.*theta./lp.iL);
@@ -103,8 +128,23 @@ end
         res_th[i] = (6/(Rpp^2)) * ( 2*(Mpp - sinh(omega))*(Mpp + sinh(omega))*exp(-2*omega*lp.iL[4]) 
         - Mpp*((Mpp + sinh(omega))*exp(-2*omega*(i-1)) + (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1)))))
     end
+        #THEORETICAL SOLUTION: hep-lat/9606016 eq (2.32)
 
-end
+    @timeit "fA analitical solution" begin
+        res_th = zeros(lp.iL[4])
+
+        pp3 = ntuple(i -> theta[i]/lp.iL[i],3)
+        omega = 2 * asinh(0.5* sqrt(( sum((sin.(pp3)).^2) + (m + 2*(sum((sin.(pp3./2)).^2) ))^2) / (1+m+2*(sum((sin.(pp3./2)).^2) )) ) )
+        pp = (-im*omega,pp3...)
+        Mpp = m + 2* sum((sin.(pp./2)).^2)
+        Rpp = Mpp*(1-exp(-2*omega*lp.iL[4])) + sinh(omega) * (1+exp(-2*omega*lp.iL[4]))
+
+        for i in 2:lp.iL[4]
+            res_th[i] = (6/(Rpp^2)) * ( 2*(Mpp - sinh(omega))*(Mpp + sinh(omega))*exp(-2*omega*lp.iL[4])
+            - Mpp*((Mpp + sinh(omega))*exp(-2*omega*(i-1)) + (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1)))))
+        end
+
+    end
 
     return sum(abs.(res-res_th))
 

test/dirac/test_solver_plw.jl

@@ -7,7 +7,7 @@ using LatticeGPU, CUDA, TimerOutputs
 function Dwpw_test(;p=0,s=1,c=1)
 lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), 0, (0,0,0,0,0,0))
 gp = GaugeParm{Float64}(SU3{Float64}, 6.0, 1.0)
-dpar = DiracParam{Float64}(SU3fund,1.3,0.0,(1.0,1.0,1.0,1.0),0.0)
+dpar = DiracParam{Float64}(SU3fund,1.3,0.0,(1.0,1.0,1.0,1.0),0.0,0.0)
 dws = DiracWorkspace(SU3fund,Float64,lp);
 
 p==0 ? p = Int.(round.(lp.iL.*rand(4),RoundUp)) : nothing
@@ -96,15 +96,15 @@ end
 
 
 begin
-dif = 0.0
+diff = 0.0
 for i in 1:3 for j in 1:4
-    dif += Dwpw_test(c=i,s=j)
+    global diff += Dwpw_test(c=i,s=j)
 end end
 
-if dif < 1.0e-15
+if diff < 1.0e-15
-    print("Dwpl test passed with average error ", dif/12,"!\n")    
+    print("Dwpl test passed with average error ", diff/12,"!\n")
 else
-    error("Dwpl test failed with difference: ",dif,"\n")
+    error("Dwpl test failed with difference: ",diff,"\n")
 end
 
 

test/dirac/test_solver_rand.jl

@@ -9,7 +9,7 @@ using CUDA, LatticeGPU, TimerOutputs
     lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), 0, (0,0,0,0,0,0))
     gp = GaugeParm{Float64}(SU3{Float64}, 6.0, 1.0)
     ymws = YMworkspace(SU3, Float64, lp)
-    dpar = DiracParam{Float64}(SU3fund,2.3,0.0,(1.0,1.0,1.0,1.0),0.0)
+    dpar = DiracParam{Float64}(SU3fund,2.3,0.0,(1.0,1.0,1.0,1.0),0.0,0.0)
     dws = DiracWorkspace(SU3fund,Float64,lp);
 
     randomize!(ymws.mom, lp, ymws)
@@ -32,7 +32,8 @@ end
 CUDA.@sync begin
         CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(rpsi)
 end
-g5Dw!(prop,U,rpsi,dpar,dws,lp)
+
+g5Dw!(prop,U,rpsi,mtwmdpar(dpar),dws,lp)
 CG!(prop,U,DwdagDw!,dpar,lp,dws,10000,1.0e-14)
 
 Dw!(dws.sp,U,prop,dpar,dws,lp)

test/runtests.jl

@@ -1,3 +1,6 @@
 
-#include("SAD/test_sad.jl")
+include("SAD/test_sad.jl")
 include("flow/test_adapt.jl")
+include("dirac/test_fp_fa.jl")
+include("dirac/test_solver_plw.jl")
+include("dirac/test_solver_rand.jl")