Initial copy of LGPU

LICENSE

+
+  "THE BEER-WARE LICENSE":
+Alberto Ramos wrote this file. As long as you retain this 
+notice you can do whatever you want with this stuff. If we meet some 
+day, and you think this stuff is worth it, you can buy me a beer in 
+return. <alberto.ramos@cern.ch>

Project.toml

+name = "LatticeGPU"
+uuid = "958c3683-801a-4582-9cfa-2d6e2ae1763b"
+authors = ["Alberto Ramos <alberto.ramos@ific.uv.es>"]
+version = "0.1.0"
+
+[deps]
+BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"

README.md

+# LatticeGPU.jl: A framework for lattice computations on the GPU
+
+This is a library/code to simulate the pure `SU(N)` gauge theory on
+GPUs. 
+
+

docs/make.jl

+using Documenter
+
+import Pkg
+Pkg.activate("../")
+using LatticeGPU
+
+makedocs(sitename="LatticeGPU", modules=[LatticeGPU], doctest=true,
+         pages = [
+             "LatticeGPU.jl" => "index.md", 
+             "Space-time" => "space.md",
+             "Groups and algebras" => "groups.md"
+             ], 
+         repo = "https://igit.ific.uv.es/alramos/latticegpu.jl")

docs/src/groups.md

+
+# Groups and Algebras
+
+The module `Groups` contain generic data types to deal with group and
+algebra elements. Group elements $$g\in SU(N)$$ are represented in
+some compact notation. For the case $$N=2$$ we use two complex numbers
+(Caley-Dickson representation, i.e. $$g=(z_1,z_2)$$ with
+$$|z_1|^2 + |z_2|^2=1$$). For the case $$N=3$$ we only store two
+rows of the $$N\times N$$ unitary matrix.
+
+Group operations translate straightforwardly in matrix operations,
+except for $$SU(2)$$, where we use
+```math
+\forall g=(z_1, z_2)\in SU(2)\, \Longrightarrow\, g^{-1} =
+(\overline z_1, -z_2)
+```
+and 
+```math
+\forall g=(z_1, z_2), g'=(w_1,w_2)\in SU(2)\, \Longrightarrow \, g g'= (z_1w_1 - \overline w_2 z_2, w_2z_1 + z_2\overline w_1)\,.
+```
+
+Group elements can be "casted" into arrays. The abstract type
+[`GMatrix`](@ref) represent generic $$ N\\times N$$ complex arrays.
+
+
+Algebra elements $$X \in \mathfrak{su}(N)$$ are traceless
+anti-hermitian matrices represented trough $$N^2-1$$ real numbers
+$$X^a$$. We have
+```math
+  X = \sum_{a=1}^{N^2-1} X^a T^a
+```
+where $$T^a$$ are a basis of the $$\mathfrak{su}(N)$$ algebra with
+the conventions
+```math
+  (T^a)^+ = -T^a \,\qquad
+  {\rm Tr}(T^aT^b) = -\frac{1}{2}\delta_{ab}\,.
+```
+
+If we denote by $$\{\mathbf{e}_i\}$$ the usual ortonormal basis of
+$$\mathbb R^N$$, the $$N^2-1$$ matrices $$T^a$$ can always be written in the
+following way 
+1. symmetric matrices ($$N(N-1)/2$$) 
+   ```math
+   (T^a)_{ij} = \frac{\imath}{2}(\mathbf{e}_i\otimes\mathbf{e}_j +
+   \mathbf{e}_j\otimes \mathbf{e}_i)\quad (1\le i < j \le N)
+   ```
+1. anti-symmetric matrices ($$N(N-1)/2$$) 
+   ```math
+   (T^a)_{ij} = \frac{1}{2}(\mathbf{e}_i\otimes\mathbf{e}_j -
+   \mathbf{e}_j\otimes \mathbf{e}_i)\quad (1\le i < j \le N)
+   ```
+1. diagonal matrices ($$N-1$$). With $$l=1,...,N-1$$ and $$a=l+N(N-1)$$
+   ```math
+   (T^a)_{ij} = \frac{\imath}{2}\sqrt{\frac{2}{l(l+1)}}
+   \left(\sum_{j=1}^l\mathbf{e}_j\otimes\mathbf{e}_j -
+     l\mathbf{e}_{l+1}\otimes \mathbf{e}_{l+1}\right)
+   \quad (l=1,\dots,N-1;\, a=l+N(N-1))
+   ```
+
+For example in the case of $$\mathfrak{su}(2)$$, and denoting the Pauli
+matrices by $$\sigma^a$$ we have
+```math
+  T^a = \frac{\imath}{2}\sigma^a \quad(a=1,2,3)\,,
+```
+while for $$\mathfrak{su}(3)$$ we have $$T^a=\imath \lambda^a/2$$ (with
+$$\lambda^a$$ the Gell-Mann matrices) up to a permutation.
+
+
+## Some examples
+
+```@setup exs
+import Pkg # hide
+Pkg.activate("/home/alberto/code/julia/LatticeGPU/") # hide
+using LatticeGPU # hide
+```
+
+Here we just show some example codes playing with groups and
+elements. The objective is to get an idea on how group operations 
+We can generate some random group elements.
+```@repl exs
+# Generate random groups elements, 
+# check they are actually from the grup
+g = rand(SU2{Float64})
+println("Are we in a group?: ", isgroup(g))
+g = rand(SU3{Float64})
+println("Are we in a group?: ", isgroup(g))
+```
+
+We can also do some simple operations 
+```@repl exs
+X = rand(SU3alg{Float64})
+g1 = exp(X);
+g2 = exp(X, -2.0); # e^{-2X}
+g = g1*g1*g2;
+println("Is this one? ", dev_one(g))
+```
+
+## Types 
+
+### Groups
+
+The generic interface reads
+
+```@docs
+Group
+```
+
+Concrete supported implementations are
+```@docs
+SU2
+SU3
+```
+
+### Algebras
+
+The generic interface reads
+```@docs
+Algebra
+```
+
+With the following concrete implementations
+```@docs
+SU2alg
+SU3alg
+```
+
+### GMatrix
+
+The generic interface reads
+```@docs
+GMatrix
+```
+
+With the following concrete implementations
+```@docs
+M2x2
+M3x3
+```
+
+
+## Generic `Group` methods
+
+```@docs
+inverse
+dag
+tr
+dev_one
+unitarize
+isgroup
+projalg
+```
+
+## Generic `Algebra` methods
+
+```@docs
+exp
+expm
+alg2mat
+```
+
+

docs/src/index.md

+
+# LatticeGPU.jl documentation
+
+
+```@contents
+```
+
+

docs/src/space.md

+
+# Lattice structure
+
+D-dimensional lattice points are labeled by two ordered integer
+numbers: the point-in-block index ($$b$$ in the figure below) and the
+block index ($$r$$ in the fire below). the routines [`up`](@ref) and
+[`dw`](@ref) allow you to displace to the neighboring points of the
+lattice. 
+![blocks](blocks.png "D dimensional lattice points are labeled by its
+index (a `Tuple` of integers). Given a point (by its index), there are
+routines to jump to neighboring points.")
+
+Directions are numbered fro 1 to D, and then Euclidean time is always
+assumed to be the last coordinate. Planes are also numbered from 1 to
+$$N(N-1)/2$$. The structure [`SpaceParm`](@ref)
+contains the information of the 
+lattice structure and has to be passed as argument to most routines. 
+
+```@docs
+SpaceParm
+up
+dw
+updw
+point_coord
+point_time
+point_index
+point_color
+```

src/Dirac/Dirac.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos and Carlos Pena wrote this file. As long as you retain this  
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy us a beer in 
+### return. <alberto.ramos@cern.ch> <carlos.pena@uam.es>
+###
+### file:    Dirac.jl
+### created: Thu Nov 18 17:20:24 2021
+###                               
+
+
+module Dirac
+
+using CUDA, TimerOutputs
+using ..Space
+using ..Groups
+using ..Fields
+using ..YM
+using ..Spinors
+
+struct DiracWorkspace{T}
+    sr
+    sp
+    sAp
+    st
+    
+    function DiracWorkspace(::Type{G}, ::Type{T}, lp::SpaceParm{4,6,B,D}) where {G,T <: AbstractFloat, B,D}
+
+        sr  = scalar_field(Spinor{4,G}, lp)
+        sp  = scalar_field(Spinor{4,G}, lp)
+        sAp = scalar_field(Spinor{4,G}, lp)
+        st  = scalar_field(Spinor{4,G}, lp)
+        return new{T}(sr,sp,sAp,st)
+    end
+end
+export DiracWorkspace
+
+function Dw!(so, U, si, m0, lp::SpaceParm)
+
+    @timeit "Dw" begin
+        CUDA.@sync begin
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, m0, [1,1,1,1], lp)
+        end
+    end
+    
+    return nothing
+end
+
+function DwdagDw!(so, U, si, m0, st, lp::SpaceParm)
+
+    @timeit "DwdagDw" begin
+        CUDA.@sync begin
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(st, U, si, m0, [1,1,1,1], lp)
+        end
+        CUDA.@sync begin
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, st, m0, [1,1,1,1], lp)
+        end
+    end
+    
+    return nothing
+end
+
+
+function krnl_Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
+
+    b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
+
+    # For SF:
+    #  - cttilde affects mass term at x0 = a, T-a
+    #  - Spinor can be periodic as long as 0 at x_0=0
+    @inbounds begin 
+        so[b,r] = (4+m0)*si[b,r]
+        for id in 1:4
+            bu, ru = up((b,r), id, lp)
+            bd, rd = dw((b,r), id, lp)
+            
+            so[b,r] -= ( th[id]*gpmul(Pgamma{id,-1},U[b,id,r],si[bu,ru]) +
+                conj(th[id])*gdagpmul(Pgamma{id,+1},U[bd,id,rd],si[bd,rd]) )/2
+        end
+    end
+
+    return nothing
+end
+
+function krnl_g5Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
+
+    b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
+
+    @inbounds begin 
+        so[b,r] = (4+m0)*si[b,r]
+        for id in 1:4
+            bu, ru = up((b,r), id, lp)
+            bd, rd = dw((b,r), id, lp)
+            
+            so[b,r] -= ( th[id]*gpmul(Pgamma{id,-1},U[b,id,r],si[bu,ru]) +
+                conj(th[id])*gdagpmul(Pgamma{id,+1},U[bd,id,rd],si[bd,rd]) )/2
+        end
+        so[b,r] = dmul(Gamma{5}, so[b,r])
+    end
+        
+    return nothing
+end
+
+export Dw!, DwdagDw!
+
+end

src/Fields/Fields.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    Fields.jl
+### created: Wed Oct  6 17:37:03 2021
+###                               
+
+module Fields
+
+using CUDA
+using ..Space
+
+
+vector_field(::Type{T}, lp::SpaceParm)     where {T} = CuArray{T, 3}(undef, lp.bsz, lp.ndim, lp.rsz)
+scalar_field(::Type{T}, lp::SpaceParm)     where {T} = CuArray{T, 2}(undef, lp.bsz, lp.rsz)
+nscalar_field(::Type{T}, n, lp::SpaceParm) where {T} = CuArray{T, 3}(undef, lp.bsz, n, lp.rsz)
+
+scalar_field_point(::Type{T}, lp::SpaceParm{N,M,D}) where {T,N,M,D} = CuArray{T, N}(undef, lp.iL...)
+
+export vector_field, scalar_field, nscalar_field, scalar_field_point
+
+end
+

src/Groups/AlgebraSU2.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    AlgebraSU2.jl
+### created: Sun Oct  3 09:24:25 2021
+###                               
+
+SU2alg(x::T)                       where T <: AbstractFloat = SU2alg{T}(x,0.0,0.0)
+SU2alg(v::Vector{T})               where T <: AbstractFloat = SU2alg{T}(v[1],v[2],v[3])
+
+"""
+    projalg([z::Complex,] g::T) where T <: Union{Group, GMatrix}
+
+Projects the group element/matrix `g` (or `zg`) to the algebra. This amounts to the following operation:
+
+\`\` X = {\\rm projalg}(g) = -X^a\\frac{1}{2} {\\rm tr}\\{T^a(g - g^{\\dagger})\\} \`\`
+
+where the substraction of group elements has to be understood in the matrix sense. This method returns an [`Algebra`](@ref) type. The generators of the algebra are the elements \`\` T^a\`\`. 
+"""
+projalg(g::SU2{T})                 where T <: AbstractFloat = SU2alg{T}(imag(g.t2), real(g.t2), imag(g.t1))
+projalg(z::Complex{T}, g::SU2{T})  where T <: AbstractFloat = SU2alg{T}(imag(z*g.t2), real(z*g.t2), imag(z*g.t1))
+dot(a::SU2alg{T}, b::SU2alg{T})    where T <: AbstractFloat = a.t1*b.t1 + a.t2*b.t2 + a.t3*b.t3
+norm(a::SU2alg{T})                 where T <: AbstractFloat = sqrt(a.t1^2 + a.t2^2 + a.t3^2)
+norm2(a::SU2alg{T})                where T <: AbstractFloat = a.t1^2 + a.t2^2 + a.t3^2
+
+Base.:+(a::SU2alg{T})              where T <: AbstractFloat = SU2alg{T}(a.t1,a.t2,a.t3)
+Base.:-(a::SU2alg{T})              where T <: AbstractFloat = SU2alg{T}(-a.t1,-a.t2,-a.t3)
+Base.:+(a::SU2alg{T},b::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(a.t1+b.t1,a.t2+b.t2,a.t3+b.t3)
+Base.:-(a::SU2alg{T},b::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(a.t1-b.t1,a.t2-b.t2,a.t3-b.t3)
+
+Base.:*(a::SU2alg{T},b::Number)    where T <: AbstractFloat = SU2alg{T}(a.t1*b,a.t2*b,a.t3*b)
+Base.:*(b::Number,a::SU2alg{T})    where T <: AbstractFloat = SU2alg{T}(a.t1*b,a.t2*b,a.t3*b)
+Base.:/(a::SU2alg{T},b::Number)    where T <: AbstractFloat = SU2alg{T}(a.t1/b,a.t2/b,a.t3/b)
+
+"""
+    alg2mat(a::T) where T <: Algebra
+
+Returns the [`Algebra`](@ref) element as a [`GMatrix`](@ref) type.
+"""
+function alg2mat(a::SU2alg{T}) where T <: AbstractFloat
+
+    u11::Complex{T} = complex(0.0, a.t3)/2
+    u22::Complex{T} = conj(u11)
+    u12::Complex{T} = complex(a.t2,a.t1)/2
+    u21::Complex{T} = -conj(u12)
+    
+    return M2x2{T}(u11,u12,u21,u22)
+end
+
+Base.:*(a::SU2alg,b::SU2) = alg2mat(a)*b
+Base.:*(a::SU2,b::SU2alg) = a*alg2mat(b)
+Base.:/(a::SU2alg,b::SU2) = alg2mat(a)/b
+Base.:\(a::SU2,b::SU2alg) = a\alg2mat(b)
+
+"""
+    exp(a::T, t::Number=1) where {T <: Algebra}
+
+Computes `exp(ta)`
+"""
+function Base.exp(a::SU2alg{T}) where T <: AbstractFloat
+    
+    rm = sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0
+    if (abs(rm) < 0.05)
+        rms = rm^2/2.0
+        ca = 1.0 - rms    *(1.0 - (rms/6.0 )*(1.0 - rms/15.0))
+        sa = 0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0))
+    else
+        ca = CUDA.cos(rm)
+	sa = CUDA.sin(rm)/(2.0*rm)
+    end
+
+    t1 = complex(ca,sa*a.t3)
+    t2 = complex(sa*a.t2,sa*a.t1)
+    return SU2{T}(t1,t2)
+end
+
+function Base.exp(a::SU2alg{T}, t::T) where T <: AbstractFloat
+    
+    rm = t*sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0
+    if (abs(rm) < 0.05)
+        rms = rm^2/2.0
+        ca = 1.0 - rms    *(1.0 - (rms/6.0 )*(1.0 - rms/15.0))
+        sa = t*(0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0)))
+    else
+        ca = CUDA.cos(rm)
+	sa = t*CUDA.sin(rm)/(2.0*rm)
+    end
+
+    t1 = complex(ca,sa*a.t3)
+    t2 = complex(sa*a.t2,sa*a.t1)
+    return SU2{T}(t1,t2)
+end
+
+
+"""
+    expm(g::G, a::A, t=1) where {G <: Group, A <: Algebra}
+
+Computes `g*exp(ta)`
+
+"""
+function expm(g::SU2{T}, a::SU2alg{T}) where T <: AbstractFloat
+    
+    rm = sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0
+    if (abs(rm) < 0.05)
+        rms = rm^2/2.0
+        ca = 1.0 - rms    *(1.0 - (rms/6.0 )*(1.0 - rms/15.0))
+        sa = 0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0))
+    else
+        ca = CUDA.cos(rm)
+	sa = CUDA.sin(rm)/(2.0*rm)
+    end
+
+    t1 = complex(ca,sa*a.t3)*g.t1-complex(sa*a.t2,sa*a.t1)*conj(g.t2)
+    t2 = complex(ca,sa*a.t3)*g.t2+complex(sa*a.t2,sa*a.t1)*conj(g.t1)
+    return SU2{T}(t1,t2)
+end
+
+function expm(g::SU2{T}, a::SU2alg{T}, t::T) where T <: AbstractFloat
+    
+    rm = t*sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0
+    if (abs(rm) < 0.05)
+        rms = rm^2/2.0
+        ca = 1.0 - rms    *(1.0 - (rms/6.0 )*(1.0 - rms/15.0))
+        sa = t*(0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0)))
+    else
+        ca = CUDA.cos(rm)
+	sa = t*CUDA.sin(rm)/(2.0*rm)
+    end
+
+    t1 = complex(ca,sa*a.t3)*g.t1-complex(sa*a.t2,sa*a.t1)*conj(g.t2)
+    t2 = complex(ca,sa*a.t3)*g.t2+complex(sa*a.t2,sa*a.t1)*conj(g.t1)
+    return SU2{T}(t1,t2)
+               
+end

src/Groups/AlgebraSU3.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    AlgebraSU3.jl
+### created: Sun Oct  3 09:13:07 2021
+###                               
+
+function projalg(a::SU3{T}) where T <: AbstractFloat
+
+    sr3ov2::T = 0.866025403784438646763723170752
+
+    ditr = ( imag(a.u11) + imag(a.u22) + 2.0*imag(a.u11*a.u22 - a.u12*a.u21) )/3.0
+    m12 = (a.u12 - conj(a.u21))/2.0
+    m13 = (a.u13 - (a.u12*a.u23 - a.u13*a.u22) )/2.0
+    m23 = (a.u23 - (a.u13*a.u21 - a.u11*a.u23) )/2.0
+
+    return SU3alg{T}(imag( m12 ), imag( m13 ), imag( m23 ),
+                     real( m12 ), real( m13 ), real( m23 ),
+                     (imag(a.u11)-imag(a.u22))/2.0,
+                     sr3ov2*(ditr))
+end
+
+function projalg(z::Complex{T}, a::SU3{T}) where T <: AbstractFloat
+
+    sr3ov2::T = 0.866025403784438646763723170752
+
+    zu11 = z*a.u11
+    zu12 = z*a.u12
+    zu13 = z*a.u13
+    zu21 = z*a.u21
+    zu22 = z*a.u22
+    zu23 = z*a.u23
+
+    ditr = ( imag(zu11) + imag(zu22) - 2.0*imag(z*conj(a.u11*a.u22 - a.u12*a.u21)) )/3.0
+    m12 = (zu12 - conj(zu21))/2.0
+    m13 = (zu13 - conj(z)*(a.u12*a.u23 - a.u13*a.u22) )/2.0
+    m23 = (zu23 - conj(z)*(a.u13*a.u21 - a.u11*a.u23) )/2.0
+
+    return SU3alg{T}(imag( m12 ), imag( m13 ), imag( m23 ),
+                     real( m12 ), real( m13 ), real( m23 ),
+                     (imag(zu11)-imag(zu22))/2.0,
+                     sr3ov2*(ditr))
+end
+
+dot(a::SU3alg{T},b::SU3alg{T})     where T <: AbstractFloat = a.t1*b.t1 + a.t2*b.t2 + a.t3*b.t3 + a.t4*b.t4 + a.t5*b.t5 + a.t6*b.t6 + a.t7*b.t7 + a.t8*b.t8
+norm2(a::SU3alg{T})                where T <: AbstractFloat = a.t1^2 + a.t2^2 + a.t3^2 + a.t4^2 + a.t5^2 + a.t6^2 + a.t7^2 + a.t8^2
+norm(a::SU3alg{T})                 where T <: AbstractFloat = sqrt(a.t1^2 + a.t2^2 + a.t3^2 + a.t4^2 + a.t5^2 + a.t6^2 + a.t7^2 + a.t8^2)
+Base.zero(::Type{SU3alg{T}})       where T <: AbstractFloat = SU3alg{T}(zero(T),zero(T),zero(T),zero(T),zero(T),zero(T),zero(T),zero(T))
+
+Base.:+(a::SU3alg{T})              where T <: AbstractFloat = SU3alg{T}(a.t1,a.t2,a.t3,a.t4,a.t5,a.t6,a.t7,a.t8)
+Base.:-(a::SU3alg{T})              where T <: AbstractFloat = SU3alg{T}(-a.t1,-a.t2,-a.t3,-a.t4,-a.t5,-a.t6,-a.t7,-a.t8)
+Base.:+(a::SU3alg{T},b::SU3alg{T}) where T <: AbstractFloat = SU3alg{T}(a.t1+b.t1,a.t2+b.t2,a.t3+b.t3,a.t4+b.t4,a.t5+b.t5,a.t6+b.t6,a.t7+b.t7,a.t8+b.t8)
+Base.:-(a::SU3alg{T},b::SU3alg{T}) where T <: AbstractFloat = SU3alg{T}(a.t1-b.t1,a.t2-b.t2,a.t3-b.t3,a.t4-b.t4,a.t5-b.t5,a.t6-b.t6,a.t7-b.t7,a.t8-b.t8)
+Base.:*(a::SU3alg{T},b::Number)    where T <: AbstractFloat = SU3alg{T}(b*a.t1,b*a.t2,b*a.t3,b*a.t4,b*a.t5,b*a.t6,b*a.t7,b*a.t8)
+Base.:*(b::Number,a::SU3alg{T})    where T <: AbstractFloat = SU3alg{T}(b*a.t1,b*a.t2,b*a.t3,b*a.t4,b*a.t5,b*a.t6,b*a.t7,b*a.t8)
+Base.:/(a::SU3alg{T},b::Number)    where T <: AbstractFloat = SU3alg{T}(a.t1/b,a.t2/b,a.t3/b,a.t4/b,a.t5/b,a.t6/b,a.t7/b,a.t8/b)
+
+Base.:*(a::SU3alg{T},b::SU3alg{T}) where T = convert(M3x3{T}, a)*convert(M3x3{T}, a)
+
+
+function alg2mat(a::SU3alg{T}) where T <: AbstractFloat
+
+    two::T = 2.0
+    rct::T = 3.46410161513775458
+
+    x8p::T = a.t8/rct
+    x7p::T = a.t7/two
+    u11::Complex{T} = complex(0.0, x7p + x8p)
+    u22::Complex{T} = complex(0.0,-x7p + x8p)
+    u33::Complex{T} = complex(0.0,-2.0*x8p)
+    u12::Complex{T} = complex(a.t4,a.t1)/two
+    u13::Complex{T} = complex(a.t5,a.t2)/two
+    u23::Complex{T} = complex(a.t6,a.t3)/two
+    u21::Complex{T} = -conj(u12)
+    u31::Complex{T} = -conj(u13)
+    u32::Complex{T} = -conj(u23)
+    
+    return M3x3{T}(u11,u12,u13, u21,u22,u23, u31,u32,u33)
+end
+
+Base.:*(a::SU3alg,b::SU3) = alg2mat(a)*b
+Base.:*(a::SU3,b::SU3alg) = a*alg2mat(b)
+Base.:/(a::SU3alg,b::SU3) = alg2mat(a)/b
+Base.:\(a::SU3,b::SU3alg) = a\alg2mat(b)
+
+@inline function exp_iter(dch::Complex{T}, tch::T) where T <: AbstractFloat
+
+    c::NTuple{22, T} = ( 1.957294106339126128e-20, 4.110317623312164853e-19,
+                         8.220635246624329711e-18, 1.561920696858622643e-16,
+                         2.811457254345520766e-15, 4.779477332387385293e-14,
+                         7.647163731819816473e-13, 1.147074559772972473e-11,
+                         1.605904383682161451e-10, 2.087675698786809894e-09,
+                         2.505210838544171879e-08, 2.755731922398589067e-07,
+                         2.755731922398589065e-06, 2.480158730158730158e-05,
+                         1.984126984126984127e-04, 1.388888888888888888e-03,
+                         8.333333333333333333e-03, 4.166666666666666666e-02,
+                         1.666666666666666666e-01, 0.5, 1.0, 1.0 )
+    
+    q0 = complex(c[1])
+    q1 = complex(0.0)
+    q2 = complex(0.0)
+    @inbounds for i in 2:length(c)
+        qt0 = q0
+        qt1 = q1
+        q0  = complex(c[i]) + dch*q2
+        q1  = qt0 - tch*q2
+        q2  = qt1 
+    end 
+    
+    return q0, q1, q2
+end
+
+    
+function expm(g::SU3{T}, a::SU3alg{T}, t::Number) where T <: AbstractFloat
+
+    tpw = t^2
+    M = alg2mat(a)
+    Msq = M*M
+    dch::Complex{T} = tpw*t*(M.u11*M.u22*M.u33 + M.u13*M.u21*M.u32 +
+                             M.u31*M.u12*M.u23 - M.u11*M.u23*M.u32 -
+                             M.u12*M.u21*M.u33 - M.u13*M.u22*M.u31)
+    tch::T = -tpw*(real(Msq.u11)+real(Msq.u22)+real(Msq.u33))/2.0
+
+    q0, q1, q2 = exp_iter(dch, tch)
+    
+    q1 = t*q1
+    q2 = tpw*q2
+
+    g2 = SU3{T}(q1*M.u11 + q2*Msq.u11+q0, q1*M.u12 + q2*Msq.u12, q1*M.u13 + q2*Msq.u13,
+                q1*M.u21 + q2*Msq.u21, q1*M.u22 + q2*Msq.u22+q0, q1*M.u23 + q2*Msq.u23)
+
+    return g2*g
+end
+
+function expm(g::SU3{T}, a::SU3alg{T}) where T <: AbstractFloat
+
+    M = alg2mat(a)
+    Msq = M*M
+    dch::Complex{T} = M.u11*M.u22*M.u33 + M.u13*M.u21*M.u32 +
+        M.u31*M.u12*M.u23 - M.u11*M.u23*M.u32 -
+        M.u12*M.u21*M.u33 - M.u13*M.u22*M.u31
+    tch::T = -(real(Msq.u11)+real(Msq.u22)+real(Msq.u33))/2.0
+    
+    q0, q1, q2 = exp_iter(dch, tch)
+    
+    g2 = SU3{T}(q1*M.u11 + q2*Msq.u11+q0, q1*M.u12 + q2*Msq.u12, q1*M.u13 + q2*Msq.u13,
+                q1*M.u21 + q2*Msq.u21, q1*M.u22 + q2*Msq.u22+q0, q1*M.u23 + q2*Msq.u23)
+    
+    return g2*g
+end
+
+function Base.exp(a::SU3alg{T}) where T <: AbstractFloat
+
+    M = alg2mat(a)
+    Msq = M*M
+    dch::Complex{T} = M.u11*M.u22*M.u33 + M.u13*M.u21*M.u32 +
+        M.u31*M.u12*M.u23 - M.u11*M.u23*M.u32 -
+        M.u12*M.u21*M.u33 - M.u13*M.u22*M.u31
+    tch::T = -(real(Msq.u11)+real(Msq.u22)+real(Msq.u33))/2.0
+
+    q0, q1, q2 = exp_iter(dch, tch)
+    
+    g2 = SU3{T}(q1*M.u11 + q2*Msq.u11+q0, q1*M.u12 + q2*Msq.u12, q1*M.u13 + q2*Msq.u13,
+                q1*M.u21 + q2*Msq.u21, q1*M.u22 + q2*Msq.u22+q0, q1*M.u23 + q2*Msq.u23)
+
+    return g2
+end
+
+function Base.exp(a::SU3alg{T}, t::Number) where T <: AbstractFloat
+
+    tpw = t^2
+    M = alg2mat(a)
+    Msq = M*M
+    dch::Complex{T} = tpw*t*(M.u11*M.u22*M.u33 + M.u13*M.u21*M.u32 +
+                             M.u31*M.u12*M.u23 - M.u11*M.u23*M.u32 -
+                             M.u12*M.u21*M.u33 - M.u13*M.u22*M.u31)
+    tch::T = -tpw*(real(Msq.u11)+real(Msq.u22)+real(Msq.u33))/2.0
+
+    q0, q1, q2 = exp_iter(dch, tch)
+    
+    q1 = t*q1
+    q2 = tpw*q2
+
+    g2 = SU3{T}(q1*M.u11 + q2*Msq.u11+q0, q1*M.u12 + q2*Msq.u12, q1*M.u13 + q2*Msq.u13,
+                q1*M.u21 + q2*Msq.u21, q1*M.u22 + q2*Msq.u22+q0, q1*M.u23 + q2*Msq.u23)
+
+    return g2
+end
+

src/Groups/FundamentalSU3.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    FundamentalSU3.jl
+### created: Tue Nov 16 15:07:00 2021
+###                               
+
+
+SU3fund(a::T, b::T, c::T)          where T <: AbstractFloat = SU3fund{T}(complex(a), complex(b), complex(c))
+
+"""
+    dag(a::SU3fund{T})
+
+Returns the conjugate of a fundamental element. 
+"""
+dag(a::SU3fund{T})                 where T <: AbstractFloat = SU3fund{T}(conj(a.t1), conj(a.t2), conj(a.t3))
+
+"""
+    norm(a::SU3fund{T})
+
+Returns the norm of a fundamental element. Same result as dot(a,a).
+"""
+norm(a::SU3fund{T})                where T <: AbstractFloat = sqrt((abs2(a.t1) + abs2(a.t2) + abs2(a.t1)))
+
+"""
+    norm(a::SU3fund{T})
+
+Returns the norm of a fundamental element. Same result as sqrt(dot(a,a)).
+"""
+norm2(a::SU3fund{T})               where T <: AbstractFloat = (abs2(a.t1) + abs2(a.t2) + abs2(a.t1))
+
+"""
+    dot(a::SU3fund{T},b::SU3fund{T})
+
+Returns the scalar product of two fundamental elements. The convention is for the product to the linear in the second argument, and anti-linear in the first argument. 
+"""
+dot(g1::SU3fund{T},g2::SU3fund{T}) where T <: AbstractFloat = conj(g1.t1)*g2.t1+g1.t2*conj(g2.t2)+g1.t3*conj(g2.t3)
+
+"""
+    *(g::SU3{T},b::SU3fund{T})
+
+Returns ga
+"""
+Base.:*(g::SU3{T},b::SU3fund{T})     where T <: AbstractFloat = SU3fund{T}(g.u11*b.t1 + g.u12*b.t2 + g.u13*b.t3,
+                                                                           g.u21*b.t1 + g.u22*b.t2 + g.u23*b.t3,
+                                                                           conj(g.u12*g.u23 - g.u13*g.u22)*b.t1 +
+                                                                               conj(g.u13*g.u21 - g.u11*g.u23)*b.t2 +
+                                                                               conj(g.u11*g.u22 - g.u12*g.u21)*b.t3)
+
+"""
+    \\(g::SU3{T},b::SU3fund{T})
+
+Returns g^dag b
+"""
+Base.:\(g::SU3{T},b::SU3fund{T})     where T <: AbstractFloat = SU3fund{T}(conj(g.u11)*b.t1 + conj(g.u21)*b.t2 + (g.u12*g.u23 - g.u13*g.u22)*b.t3,
+                                                                           conj(g.u12)*b.t1 + conj(g.u22)*b.t2 + (g.u13*g.u21 - g.u11*g.u23)*b.t3,
+                                                                           conj(g.u13)*b.t1 + conj(g.u23)*b.t2 + (g.u11*g.u22 - g.u12*g.u21)*b.t3)
+
+Base.:+(a::SU3fund{T},b::SU3fund{T}) where T <: AbstractFloat = SU3fund{T}(a.t1+b.t1,a.t2+b.t2,a.t3+b.t3)
+Base.:-(a::SU3fund{T},b::SU3fund{T}) where T <: AbstractFloat = SU3fund{T}(a.t1-b.t1,a.t2-b.t2,a.t3-b.t3)
+Base.:+(a::SU3fund{T})               where T <: AbstractFloat = SU3fund{T}(a.t1,a.t2,a.t3)
+Base.:-(a::SU3fund{T})               where T <: AbstractFloat = SU3fund{T}(-a.t1,-a.t2,-a.t3)
+imm(a::SU3fund{T})                   where T <: AbstractFloat = SU3fund{T}(complex(-imag(a.t1),real(a.t1)),
+                                                                           complex(-imag(a.t2),real(a.t2)),
+                                                                           complex(-imag(a.t3),real(a.t3)))
+mimm(a::SU3fund{T})                  where T <: AbstractFloat = SU3fund{T}(complex(imag(a.t1),-real(a.t1)),
+                                                                           complex(imag(a.t2),-real(a.t2)),
+                                                                           complex(imag(a.t3),-real(a.t3)))
+
+# Operations with numbers
+Base.:*(a::SU3fund{T},b::Number) where T <: AbstractFloat = SU3fund{T}(b*a.t1,b*a.t2,b*a.t3)
+Base.:*(b::Number,a::SU3fund{T}) where T <: AbstractFloat = SU3fund{T}(b*a.t1,b*a.t2,b*a.t3)
+Base.:/(a::SU3fund{T},b::Number) where T <: AbstractFloat = SU3fund{T}(a.t1/b,a.t2/b,a.t3/b)
+

src/Groups/GroupSU2.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    GroupSU2.jl
+### created: Sun Jul 11 17:23:12 2021
+###                               
+
+#
+# SU(2) group elements represented trough Cayley-Dickson
+#       construction
+# https://en.wikipedia.org/wiki/Cayley%E2%80%93Dickson_construction
+using CUDA, Random
+
+SU2(a::T, b::T)          where T <: AbstractFloat = SU2{T}(complex(a), complex(b))
+
+"""
+    inverse(g::T) where T <: Group
+
+Returns the group inverse of `g`.
+"""
+inverse(b::SU2{T})       where T <: AbstractFloat = SU2{T}(conj(b.t1), -b.t2)
+
+"""
+    dag(g::T) where T <: Group
+
+Returns the group inverse of `g`.
+"""
+dag(a::SU2{T})           where T <: AbstractFloat = inverse(a)
+norm(a::SU2{T})          where T <: AbstractFloat = sqrt(abs2(a.t1) + abs2(a.t2))
+norm2(a::SU2{T})         where T <: AbstractFloat = abs2(a.t1) + abs2(a.t2)
+
+"""
+    tr(g::T) where T <: Group
+
+Returns the trace of the groups element `g`.
+"""
+tr(g::SU2{T})            where T <: AbstractFloat = complex(2*real(g.t1), 0.0)
+
+"""
+    dev_one(T) where T <: Group
+
+Returns the distance to the unit group element
+"""
+dev_one(g::SU2{T})       where T <: AbstractFloat = sqrt(( abs2(g.t1 - one(T)) + abs2(g.t2))/2)
+
+"""
+    unitarize(a::T) where {T <: Group}
+
+Return a unitarized element of the group.
+"""
+function unitarize(a::SU2{T}) where T <: AbstractFloat
+    dr = sqrt(abs2(a.t1) + abs2(a.t2))
+    if (dr == 0.0)
+        return SU2{T}(0.0,0.0)
+    end
+    return SU2{T}(a.t1/dr,a.t2/dr)
+end
+
+Base.:*(a::SU2{T},b::SU2{T}) where T <: AbstractFloat = SU2{T}(a.t1*b.t1-a.t2*conj(b.t2),a.t1*b.t2+a.t2*conj(b.t1))
+Base.:/(a::SU2{T},b::SU2{T}) where T <: AbstractFloat = SU2{T}(a.t1*conj(b.t1)+a.t2*conj(b.t2),-a.t1*b.t2+a.t2*b.t1)
+Base.:\(a::SU2{T},b::SU2{T}) where T <: AbstractFloat = SU2{T}(conj(a.t1)*b.t1+a.t2*conj(b.t2),conj(a.t1)*b.t2-a.t2*conj(b.t1))
+
+"""
+    isgroup(g::T) where T <: Group
+
+Returns `true` if `g` is a group element, `false` otherwise.
+"""
+function isgroup(a::SU2{T}) where T <: AbstractFloat
+    tol = 1.0E-10
+    if (abs2(a.t1) + abs2(a.t2) - 1.0 < 1.0E-10)
+        return true
+    else
+        return false
+    end
+end
+

src/Groups/GroupSU3.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    GroupSU3.jl
+### created: Sun Jul 11 17:23:02 2021
+###                               
+
+inverse(a::SU3{T})       where T <: AbstractFloat = SU3{T}(conj(a.u11),conj(a.u21),(a.u12*a.u23 - a.u13*a.u22), conj(a.u12),conj(a.u22),(a.u13*a.u21 - a.u11*a.u23))
+dag(a::SU3{T})           where T <: AbstractFloat = inverse(a)
+tr(a::SU3{T})            where T <: AbstractFloat = a.u11+a.u22+conj(a.u11*a.u22 - a.u12*a.u21)
+dev_one(g::SU3{T})       where T <: AbstractFloat = sqrt(( abs2(g.u11 - one(T)) + abs2(g.u12) + abs2(g.u13) + abs2(g.u21) + abs2(g.u22 - one(T)) + abs2(g.u23) )/6)
+
+function unitarize(g::SU3{T}) where T <: AbstractFloat
+
+    dv = sqrt(abs2(g.u11)+abs2(g.u12)+abs2(g.u13))
+    gu11 = g.u11/dv
+    gu12 = g.u12/dv
+    gu13 = g.u13/dv
+    
+    z    = g.u21*conj(gu11) + g.u22*conj(gu12) + g.u23*conj(gu13)
+    gu21 = g.u21 - z*gu11
+    gu22 = g.u22 - z*gu12
+    gu23 = g.u23 - z*gu13
+    dv = sqrt(abs2(gu21)+abs2(gu22)+abs2(gu23))
+    
+    return SU3{T}(gu11, gu12, gu13, gu21/dv, gu22/dv, gu23/dv)
+end
+    
+function Base.:*(a::SU3{T},b::SU3{T}) where T <: AbstractFloat
+
+    bu31 = conj(b.u12*b.u23 - b.u13*b.u22)
+    bu32 = conj(b.u13*b.u21 - b.u11*b.u23)
+    bu33 = conj(b.u11*b.u22 - b.u12*b.u21)
+
+    return SU3{T}(a.u11*b.u11 + a.u12*b.u21 + a.u13*bu31,
+                  a.u11*b.u12 + a.u12*b.u22 + a.u13*bu32, 
+                  a.u11*b.u13 + a.u12*b.u23 + a.u13*bu33, 
+                  a.u21*b.u11 + a.u22*b.u21 + a.u23*bu31, 
+                  a.u21*b.u12 + a.u22*b.u22 + a.u23*bu32,
+                  a.u21*b.u13 + a.u22*b.u23 + a.u23*bu33)
+end
+
+function Base.:/(a::SU3{T},b::SU3{T}) where T <: AbstractFloat
+
+    bu31 = (b.u12*b.u23 - b.u13*b.u22)
+    bu32 = (b.u13*b.u21 - b.u11*b.u23)
+    bu33 = (b.u11*b.u22 - b.u12*b.u21)
+
+    return SU3{T}(a.u11*conj(b.u11) + a.u12*conj(b.u12) + a.u13*conj(b.u13),
+                  a.u11*conj(b.u21) + a.u12*conj(b.u22) + a.u13*conj(b.u23), 
+                  a.u11*(bu31)      + a.u12*(bu32)      + a.u13*(bu33), 
+                  a.u21*conj(b.u11) + a.u22*conj(b.u12) + a.u23*conj(b.u13), 
+                  a.u21*conj(b.u21) + a.u22*conj(b.u22) + a.u23*conj(b.u23),
+                  a.u21*(bu31)      + a.u22*(bu32)      + a.u23*(bu33))
+end
+
+function Base.:\(a::SU3{T},b::SU3{T}) where T <: AbstractFloat
+
+    au31 = (a.u12*a.u23 - a.u13*a.u22)
+    au32 = (a.u13*a.u21 - a.u11*a.u23)
+    bu31 = conj(b.u12*b.u23 - b.u13*b.u22)
+    bu32 = conj(b.u13*b.u21 - b.u11*b.u23)
+    bu33 = conj(b.u11*b.u22 - b.u12*b.u21)
+
+    return SU3{T}(conj(a.u11)*b.u11 + conj(a.u21)*b.u21 + (au31)*bu31,
+                  conj(a.u11)*b.u12 + conj(a.u21)*b.u22 + (au31)*bu32, 
+                  conj(a.u11)*b.u13 + conj(a.u21)*b.u23 + (au31)*bu33, 
+                  conj(a.u12)*b.u11 + conj(a.u22)*b.u21 + (au32)*bu31, 
+                  conj(a.u12)*b.u12 + conj(a.u22)*b.u22 + (au32)*bu32,
+                  conj(a.u12)*b.u13 + conj(a.u22)*b.u23 + (au32)*bu33)
+end
+
+function isgroup(a::SU3{T}) where T <: AbstractFloat
+
+    tol = 1.0E-10
+    g = a/a
+    if ( (abs(g.u11 - 1.0) < tol) &&
+         (abs(g.u12) < tol) &&
+         (abs(g.u13) < tol) &&
+         (abs(g.u21) < tol) &&
+         (abs(g.u22 - 1.0) < tol) &&
+         (abs(g.u23) < tol) )
+        return true
+    else
+        return false
+    end
+end
+

src/Groups/GroupU1.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    GroupU1.jl
+### created: Tue Sep 21 09:33:44 2021
+###                               
+
+using CUDA, Random
+
+import Base.:*, Base.:+, Base.:-,Base.:/,Base.:\,Base.exp,Base.zero,Base.one
+import Random.rand
+struct U1{T} <: Group
+    t1::T
+    t2::T
+end
+U1(a::T)          where T <: AbstractFloat = U1{T}(a,zero(T))
+inverse(b::U1{T}) where T <: AbstractFloat = U1{T}(b.t1,-b.t2)
+dag(a::U1{T})    where T <: AbstractFloat = inverse(a)
+norm(a::U1{T})    where T <: AbstractFloat = sqrt(a.t1^2+a.t2^2)
+norm2(a::U1{T})   where T <: AbstractFloat = a.t1^2+a.t2^2
+tr(g::U1{T})      where T <: AbstractFloat = complex(a.t1)
+Base.one(::Type{U1{T}}) where T <: AbstractFloat = U1{T}(one(T), zero(T))
+function Random.rand(rng::AbstractRNG, ::Random.SamplerType{U1{T}}) where T <: AbstractFloat
+    r = randn(rng,T)
+    return U1{T}(CUDA.cos(r),CUDA.sin(r))
+end
+
+"""
+    function normalize(a::U1)
+
+Return a normalized element of `SU(2)`
+"""
+function normalize(a::U1{T}) where T <: AbstractFloat
+    dr = norm(a)
+    return U1{T}(a.t1/dr, a.t2/dr)
+end
+
+Base.:*(a::U1{T},b::U1{T}) where T <: AbstractFloat = U1{T}(a.t1*b.t1-a.t2*b.t2, a.t1*b.t2+a.t2*b.t1)
+Base.:/(a::U1{T},b::U1{T}) where T <: AbstractFloat = U1{T}(a.t1*b.t1+a.t2*b.t2, -a.t1*b.t2+a.t2*b.t1)
+Base.:\(a::U1{T},b::U1{T}) where T <: AbstractFloat = U1{T}(a.t1*b.t1+a.t2*b.t2, a.t1*b.t2-a.t2*b.t1)
+
+struct U1alg{T} <: Algebra
+    t::T
+end
+projalg(g::U1{T})             where T <: AbstractFloat = U1alg{T}(g.t2)
+dot(a::U1alg{T}, b::U1alg{T}) where T <: AbstractFloat = a.t*b.t
+norm(a::U1alg{T})             where T <: AbstractFloat = abs(a.t)
+norm2(a::U1alg{T})            where T <: AbstractFloat = a.t^2
+Base.zero(::Type{U1alg{T}})   where T <: AbstractFloat = U1alg{T}(zero(T))
+Random.rand(rng::AbstractRNG, ::Random.SamplerType{U1alg{T}}) where T <: AbstractFloat = U1alg{T}(randn(rng,T))
+
+Base.:+(a::U1alg{T})             where T <: AbstractFloat = U1alg{T}(a.t)
+Base.:-(a::U1alg{T})             where T <: AbstractFloat = U1alg{T}(-a.t)
+Base.:+(a::U1alg{T},b::U1alg{T}) where T <: AbstractFloat = U1alg{T}(a.t+b.t)
+Base.:-(a::U1alg{T},b::U1alg{T}) where T <: AbstractFloat = U1alg{T}(a.t-b.t)
+
+Base.:*(a::U1alg{T},b::Number)   where T <: AbstractFloat = U1alg{T}(a.t*b)
+Base.:*(b::Number,a::U1alg{T})   where T <: AbstractFloat = U1alg{T}(a.t*b)
+Base.:/(a::U1alg{T},b::Number)   where T <: AbstractFloat = U1alg{T}(a.t/b)
+
+isgroup(a::U1{T}) where T <: AbstractFloat = (abs(a.t) -1.0) < 1.0E-10
+
+Base.exp(a::U1alg{T}) where T <: AbstractFloat = U1{T}(CUDA.cos(a.t), CUDA.sin(a.t))
+Base.exp(a::U1alg{T}, t::T) where T <: AbstractFloat  = U1{T}(CUDA.cos(t*a.t), CUDA.sin(t*a.t))
+
+expm(g::U1{T}, a::U1alg{T}) where T <: AbstractFloat = U1{T}(CUDA.cos(a.t), CUDA.sin(a.t))*g
+expm(g::U1{T}, a::U1alg{T}, t::T) where T <: AbstractFloat = U1{T}(CUDA.cos(t*a.t), CUDA.sin(t*a.t))*g
+
+export U1, U1alg, inverse, dag, tr, projalg, expm, exp, norm, norm2, isgroup

src/Groups/Groups.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    Groups.jl
+### created: Sun Jul 11 18:02:16 2021
+###                               
+
+
+module Groups
+
+using CUDA, Random
+import Base.:*, Base.:+, Base.:-,Base.:/,Base.:\,Base.exp,Base.one,Base.zero
+import Random.rand
+
+"""
+    abstract type Group
+
+This abstract type encapsulates group types (i.e. \`\`SU(N)\`\`). The following operations/methods are assumed to be defined for a each group.
+- `*,/,\\`: Usual group operations.
+- [`inverse`](@ref), [`dag`](@ref): The group inverse.
+- [`tr`](@ref)
+- [`dev_one`](@ref)
+- [`unitarize`](@ref)
+- [`isgroup`](@ref)
+- [`projalg`](@ref)
+"""
+abstract type Group end
+
+"""
+    abstract type Algebra
+
+This abstract type encapsulates algebra types  (i.e. \`\` {\\rm su}(N)\`\`). The following operations/methods are assumed to be defined for each algebra:
+- `+,-`: Usual operation of Algebra elements.
+- `*,/`: Multiplication/division by Numbers
+- `*`: Multiplication of [`Algebra`](@ref) and [`Group`](@ref) or [`GMatrix`](@ref) elements. These routines use the matrix form of the group elements, and return a [`GMatrix`](@ref) type.
+- [`exp`](@ref), [`expm`](@ref): Exponential map. Returns a [`Group`](@ref) element.
+"""
+abstract type Algebra end
+
+"""
+    abstract type GM
+
+This abstract type encapsulates \`\` N\\times N \`\` matrix types. The usual matrix operations (`+,-,*,/`) are expected to be defined for this case.
+"""
+abstract type GMatrix end
+
+export Group, Algebra, GMatrix
+
+##
+# SU(2) and 2x2 matrix operations
+##
+include("SU2Types.jl")
+export SU2, SU2alg, M2x2
+
+include("GroupSU2.jl")
+include("M2x2.jl")
+include("AlgebraSU2.jl")
+## END SU(2)
+
+##
+# SU(3) and 3x3 matrix operations
+##
+include("SU3Types.jl")
+export SU3, SU3alg, M3x3, SU3fund
+
+include("GroupSU3.jl")
+include("M3x3.jl")
+include("AlgebraSU3.jl")
+include("FundamentalSU3.jl")
+export imm, mimm
+## END SU(3)
+
+include("GroupU1.jl")
+export U1, U1alg
+
+
+export dot, expm, exp, dag, unitarize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat, dev_one
+
+
+end # module

src/Groups/M2x2.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    M3x3.jl
+### created: Sun Oct  3 09:03:34 2021
+###                               
+
+
+Base.:*(a::M2x2{T},b::M2x2{T}) where T <: AbstractFloat = M2x2{T}(a.u11*b.u11 + a.u12*b.u21,
+                                                                  a.u11*b.u12 + a.u12*b.u22, 
+                                                                  a.u21*b.u11 + a.u22*b.u21, 
+                                                                  a.u21*b.u12 + a.u22*b.u22)
+
+Base.:*(a::SU2{T},b::M2x2{T}) where T <: AbstractFloat = M2x2{T}(a.t1*b.u11+a.t2*b.u21,
+                                                                 a.t1*b.u12+a.t2*b.u22,
+                                                                 -conj(a.t2)*b.u11+conj(a.t1)*b.u21,
+                                                                 -conj(a.t2)*b.u12+conj(a.t1)*b.u22)
+    
+Base.:*(a::M2x2{T},b::SU2{T}) where T <: AbstractFloat = M2x2{T}(a.u11*b.t1-a.u12*conj(b.t2),
+                                                                 a.u11*b.t2+a.u12*conj(b.t1),
+                                                                 a.u21*b.t1-a.u22*conj(b.t2),
+                                                                 -a.u21*b.t2+a.u22*conj(b.t1))
+
+Base.:/(a::M2x2{T},b::SU2{T}) where T <: AbstractFloat = M2x2{T}(a.u11*conj(b.t1)-a.u12*conj(b.t2),
+                                                                 a.u11*b.t2+a.u12*b.t1,
+                                                                 a.u21*conj(b.t1)-a.u22*conj(b.t2),
+                                                                 -a.u21*b.t2+a.u22*b.t1)
+
+Base.:\(a::SU2{T},b::M2x2{T}) where T <: AbstractFloat = M2x2{T}(conj(a.t1)*b.u11+a.t2*b.u21,
+                                                                 conj(a.t1)*b.u12+a.t2*b.u22,
+                                                                 -conj(a.t2)*b.u11+a.t1*b.u21,
+                                                                 -conj(a.t2)*b.u12+a.t1*b.u22)
+
+Base.:*(a::Number,b::M2x2{T}) where T <: AbstractFloat  = M2x2{T}(a*b.u11, a*b.u12,
+                                                                  a*b.u21, a*b.u22)
+
+Base.:*(b::M2x2{T},a::Number) where T <: AbstractFloat  = M2x2{T}(a*b.u11, a*b.u12,
+                                                                  a*b.u21, a*b.u22)
+
+Base.:+(a::M2x2{T},b::M2x2{T}) where T <: AbstractFloat = M2x2{T}(a.u11+b.u11, a.u12+b.u12,
+                                                                  a.u21+b.u21, a.u22+b.u22)
+
+Base.:-(a::M2x2{T},b::M2x2{T}) where T <: AbstractFloat = M2x2{T}(a.u11-b.u11, a.u12-b.u12,
+                                                                  a.u21-b.u21, a.u22-b.u22)
+
+Base.:-(b::M2x2{T}) where T <: AbstractFloat            = M2x2{T}(-b.u11, -b.u12,
+                                                                  -b.u21, -b.u22)
+
+Base.:+(b::M2x2{T}) where T <: AbstractFloat            = M2x2{T}(b.u11, b.u12,
+                                                                  b.u21, b.u22)
+
+function projalg(a::M2x2{T}) where T <: AbstractFloat
+
+    m12 = (a.u12 - conj(a.u21))/2
+
+    return SU2alg{T}(imag( m12 ), real( m12 ), (imag(a.u11) - imag(a.u22))/2)
+end
+

src/Groups/M3x3.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    M3x3.jl
+### created: Sun Oct  3 09:03:34 2021
+###                               
+
+
+Base.:*(a::M3x3{T},b::M3x3{T}) where T <: AbstractFloat = M3x3{T}(a.u11*b.u11 + a.u12*b.u21 + a.u13*b.u31,
+                                                                  a.u11*b.u12 + a.u12*b.u22 + a.u13*b.u32, 
+                                                                  a.u11*b.u13 + a.u12*b.u23 + a.u13*b.u33, 
+                                                                  a.u21*b.u11 + a.u22*b.u21 + a.u23*b.u31, 
+                                                                  a.u21*b.u12 + a.u22*b.u22 + a.u23*b.u32,
+                                                                  a.u21*b.u13 + a.u22*b.u23 + a.u23*b.u33,
+                                                                  a.u31*b.u11 + a.u32*b.u21 + a.u33*b.u31, 
+                                                                  a.u31*b.u12 + a.u32*b.u22 + a.u33*b.u32,
+                                                                  a.u31*b.u13 + a.u32*b.u23 + a.u33*b.u33)
+
+function Base.:*(a::SU3{T},b::M3x3{T}) where T <: AbstractFloat
+    
+    au31 = conj(a.u12*a.u23 - a.u13*a.u22)
+    au32 = conj(a.u13*a.u21 - a.u11*a.u23)
+    au33 = conj(a.u11*a.u22 - a.u12*a.u21)
+
+    return M3x3{T}(a.u11*b.u11 + a.u12*b.u21 + a.u13*b.u31,
+                   a.u11*b.u12 + a.u12*b.u22 + a.u13*b.u32, 
+                   a.u11*b.u13 + a.u12*b.u23 + a.u13*b.u33, 
+                   a.u21*b.u11 + a.u22*b.u21 + a.u23*b.u31, 
+                   a.u21*b.u12 + a.u22*b.u22 + a.u23*b.u32,
+                   a.u21*b.u13 + a.u22*b.u23 + a.u23*b.u33,
+                   au31*b.u11  + au32*b.u21  + au33*b.u31, 
+                   au31*b.u12  + au32*b.u22  + au33*b.u32,
+                   au31*b.u13  + au32*b.u23  + au33*b.u33)
+end
+
+    
+function Base.:*(a::M3x3{T},b::SU3{T}) where T <: AbstractFloat
+
+    bu31 = conj(b.u12*b.u23 - b.u13*b.u22)
+    bu32 = conj(b.u13*b.u21 - b.u11*b.u23)
+    bu33 = conj(b.u11*b.u22 - b.u12*b.u21)
+
+    return M3x3{T}(a.u11*b.u11 + a.u12*b.u21 + a.u13*bu31,
+                   a.u11*b.u12 + a.u12*b.u22 + a.u13*bu32, 
+                   a.u11*b.u13 + a.u12*b.u23 + a.u13*bu33, 
+                   a.u21*b.u11 + a.u22*b.u21 + a.u23*bu31, 
+                   a.u21*b.u12 + a.u22*b.u22 + a.u23*bu32,
+                   a.u21*b.u13 + a.u22*b.u23 + a.u23*bu33,
+                   a.u31*b.u11 + a.u32*b.u21 + a.u33*bu31, 
+                   a.u31*b.u12 + a.u32*b.u22 + a.u33*bu32,
+                   a.u31*b.u13 + a.u32*b.u23 + a.u33*bu33)
+end
+
+function Base.:/(a::M3x3{T},b::SU3{T}) where T <: AbstractFloat
+
+    bu31 = (b.u12*b.u23 - b.u13*b.u22)
+    bu32 = (b.u13*b.u21 - b.u11*b.u23)
+    bu33 = (b.u11*b.u22 - b.u12*b.u21)
+
+    return M3x3{T}(a.u11*conj(b.u11) + a.u12*conj(b.u12) + a.u13*conj(b.u13),
+                   a.u11*conj(b.u21) + a.u12*conj(b.u22) + a.u13*conj(b.u23), 
+                   a.u11*(bu31) + a.u12*(bu32) + a.u13*(bu33), 
+                   a.u21*conj(b.u11) + a.u22*conj(b.u12) + a.u23*conj(b.u13), 
+                   a.u21*conj(b.u21) + a.u22*conj(b.u22) + a.u23*conj(b.u23),
+                   a.u21*(bu31) + a.u22*(bu32) + a.u23*(bu33),
+                   a.u31*conj(b.u11) + a.u32*conj(b.u12) + a.u33*conj(b.u13), 
+                   a.u31*conj(b.u21) + a.u32*conj(b.u22) + a.u33*conj(b.u23),
+                   a.u31*(bu31) + a.u32*(bu32) + a.u33*(bu33))
+end
+Base.:/(a::M3x3{T}, b::Number) where T <: AbstractFloat = M3x3{T}(a.u11/b, a.u12/b, a.u13/b,
+                                                                  a.u21/b, a.u22/b, a.u23/b,
+                                                                  a.u31/b, a.u32/b, a.u33/b)
+
+function Base.:\(a::SU3{T},b::M3x3{T}) where T <: AbstractFloat
+
+    au31 = (a.u12*a.u23 - a.u13*a.u22)
+    au32 = (a.u13*a.u21 - a.u11*a.u23)
+    au33 = (a.u11*a.u22 - a.u12*a.u21)
+
+    return M3x3{T}(conj(a.u11)*b.u11 + conj(a.u21)*b.u21 + (au31)*b.u31,
+                   conj(a.u11)*b.u12 + conj(a.u21)*b.u22 + (au31)*b.u32, 
+                   conj(a.u11)*b.u13 + conj(a.u21)*b.u23 + (au31)*b.u33, 
+                   conj(a.u12)*b.u11 + conj(a.u22)*b.u21 + (au32)*b.u31, 
+                   conj(a.u12)*b.u12 + conj(a.u22)*b.u22 + (au32)*b.u32,
+                   conj(a.u12)*b.u13 + conj(a.u22)*b.u23 + (au32)*b.u33,
+                   conj(a.u13)*b.u11 + conj(a.u23)*b.u21 + (au33)*b.u31, 
+                   conj(a.u13)*b.u12 + conj(a.u23)*b.u22 + (au33)*b.u32,
+                   conj(a.u13)*b.u13 + conj(a.u23)*b.u23 + (au33)*b.u33)
+
+end
+
+Base.:*(a::Number,b::M3x3{T}) where T <: AbstractFloat  = M3x3{T}(a*b.u11, a*b.u12, a*b.u13,
+                                                                  a*b.u21, a*b.u22, a*b.u23,
+                                                                  a*b.u31, a*b.u32, a*b.u33)
+
+Base.:*(b::M3x3{T},a::Number) where T <: AbstractFloat  = M3x3{T}(a*b.u11, a*b.u12, a*b.u13,
+                                                                  a*b.u21, a*b.u22, a*b.u23,
+                                                                  a*b.u31, a*b.u32, a*b.u33)
+
+Base.:+(a::M3x3{T},b::M3x3{T}) where T <: AbstractFloat = M3x3{T}(a.u11+b.u11, a.u12+b.u12, a.u13+b.u13,
+                                                                  a.u21+b.u21, a.u22+b.u22, a.u23+b.u23,
+                                                                  a.u31+b.u31, a.u32+b.u32, a.u33+b.u33)
+
+Base.:-(a::M3x3{T},b::M3x3{T}) where T <: AbstractFloat = M3x3{T}(a.u11-b.u11, a.u12-b.u12, a.u13-b-u13,
+                                                                  a.u21-b.u21, a.u22-b.u22, a.u23-b-u23,
+                                                                  a.u31-b.u31, a.u32-b.u32, a.u33-b-u33)
+
+Base.:-(b::M3x3{T}) where T <: AbstractFloat            = M3x3{T}(-b.u11, -b.u12, -b.u13,
+                                                                  -b.u21, -b.u22, -b.u23,
+                                                                  -b.u31, -b.u32, -b.u33)
+
+Base.:+(b::M3x3{T}) where T <: AbstractFloat            = M3x3{T}(b.u11, b.u12, b.u13,
+                                                                  b.u21, b.u22, b.u23,
+                                                                  b.u31, b.u32, b.u33)
+function projalg(a::M3x3{T}) where T <: AbstractFloat
+
+    sr3ov2::T = 0.866025403784438646763723170752
+
+    ditr = ( imag(a.u11) + imag(a.u22) - 2.0*imag(a.u33) )/3.0
+    m12 = (a.u12 - conj(a.u21))/2.0
+    m13 = (a.u13 - conj(a.u31))/2.0
+    m23 = (a.u23 - conj(a.u32))/2.0
+
+    return SU3alg{T}(imag( m12 ), imag( m13 ), imag( m23 ),
+                     real( m12 ), real( m13 ), real( m23 ),
+                     (imag(a.u11)-imag(a.u22))/2.0,
+                     sr3ov2*(ditr))
+end
+
+dag(a::M3x3{T}) where T = M3x3{T}(conj(a.u11), conj(a.u21), conj(a.u31),
+                                  conj(a.u12), conj(a.u22), conj(a.u32),
+                                  conj(a.u13), conj(a.u23), conj(a.u33))
+tr(a::M3x3{T}) where T = a.u11+a.u22+a.u33

src/Groups/SU2Types.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    SU2Types.jl
+### created: Sun Oct  3 09:22:48 2021
+###                               
+
+
+"""
+    struct SU2{T} <: Group
+
+\`\`SU(2)\`\` elements. The type `T <: AbstractFloat` can be used to define single or double precision elements.
+"""
+struct SU2{T} <: Group
+    t1::Complex{T}
+    t2::Complex{T}
+end
+
+"""
+    struct M2x2{T} <: GMatrix
+
+Generic \`\` 2\\times 2\`\` complex matrix.
+"""
+struct M2x2{T} <: GMatrix
+    u11::Complex{T}
+    u12::Complex{T}
+    u21::Complex{T}
+    u22::Complex{T}
+end
+
+"""
+    struct SU2alg{T} <: Algebra
+
+\`\`{\\rm su}(2)\`\` Algebra elements.
+"""
+struct SU2alg{T} <: Algebra
+    t1::T
+    t2::T
+    t3::T
+end
+
+Base.zero(::Type{SU2alg{T}}) where T <: AbstractFloat = SU2alg{T}(zero(T),zero(T),zero(T))
+Base.zero(::Type{M2x2{T}})   where T <: AbstractFloat = M2x2{T}(zero(T),zero(T),zero(T),zero(T))
+Base.one(::Type{SU2{T}})     where T <: AbstractFloat = SU2{T}(one(T),zero(T))
+Base.one(::Type{M2x2{T}})    where T <: AbstractFloat = M2x2{T}(one(T),zero(T),zero(T),one(T))
+
+Base.convert(::Type{M2x2{T}}, a::SU2alg{T}) where T = alg2mat(a)
+Base.convert(::Type{M2x2{T}}, a::SU2{T}) where T = M2x2{T}(a.t1,a.t2,-conj(a.t2), conj(a.t1))
+
+Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU2alg{T}}) where T <: AbstractFloat = SU2alg{T}(randn(rng,T),randn(rng,T),randn(rng,T))
+Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU2{T}})    where T <: AbstractFloat = exp(SU2alg{T}(randn(rng,T),randn(rng,T),randn(rng,T)))

src/Groups/SU3Types.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    SU3Types.jl
+### created: Sun Oct  3 09:05:23 2021
+###                               
+
+#
+# Use memory efficient representation: Only store
+# first two rows. Third row constructed on the fly.
+#
+# a.u31 = conj(a.u12*a.u23 - a.u13*a.u22)
+# a.u32 = conj(a.u13*a.u21 - a.u11*a.u23)
+# a.u33 = conj(a.u11*a.u22 - a.u12*a.u21)
+#
+
+import Base.convert
+
+"""
+    struct SU3{T} <: Group
+
+\`\`SU(3)\`\` elements. The type `T <: AbstractFloat` can be used to define single or double precision elements.
+"""
+struct SU3{T} <: Group
+    u11::Complex{T}
+    u12::Complex{T}
+    u13::Complex{T}
+    u21::Complex{T}
+    u22::Complex{T}
+    u23::Complex{T}
+end
+Base.one(::Type{SU3{T}}) where T <: AbstractFloat = SU3{T}(one(T),zero(T),zero(T),zero(T),one(T),zero(T))
+
+"""
+    struct M3x3{T} <: Group
+
+3x3 Complex matrix. The type `T <: AbstractFloat` can be used to define single or double precision elements.
+"""
+struct M3x3{T} <: GMatrix
+    u11::Complex{T}
+    u12::Complex{T}
+    u13::Complex{T}
+    u21::Complex{T}
+    u22::Complex{T}
+    u23::Complex{T}
+    u31::Complex{T}
+    u32::Complex{T}
+    u33::Complex{T}
+end
+Base.one(::Type{M3x3{T}}) where T <: AbstractFloat = M3x3{T}(one(T),zero(T),zero(T),zero(T),one(T),zero(T),zero(T),zero(T),one(T))
+Base.zero(::Type{M3x3{T}}) where T <: AbstractFloat = M3x3{T}(zero(T),zero(T),zero(T),zero(T),zero(T),zero(T),zero(T),zero(T),zero(T))
+
+"""
+    struct SU3alg{T} <: Group
+
+\`\`su(3)\`\` elements. The type `T <: AbstractFloat` can be used to define single or double precision elements.
+"""
+struct SU3alg{T} <: Algebra
+    t1::T
+    t2::T
+    t3::T
+    t4::T
+    t5::T
+    t6::T
+    t7::T
+    t8::T
+end
+
+
+Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU3{T}}) where T <: AbstractFloat = exp(SU3alg{T}(randn(rng,T),randn(rng,T),randn(rng,T),randn(rng,T),randn(rng,T),randn(rng,T),randn(rng,T),randn(rng,T)))
+Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU3alg{T}}) where T <: AbstractFloat = SU3alg{T}(randn(rng,T),randn(rng,T),randn(rng,T),randn(rng,T),randn(rng,T),randn(rng,T),randn(rng,T),randn(rng,T))
+
+Base.convert(::Type{M3x3{T}}, a::SU3alg{T}) where T = alg2mat(a)
+Base.convert(::Type{M3x3{T}}, a::SU3{T}) where T = M3x3{T}(a.u11,a.u12,a.u13,
+                                                        a.u21,a.u22,a.u23,
+                                                        conj(a.u12*a.u23 - a.u13*a.u22),
+                                                        conj(a.u13*a.u21 - a.u11*a.u23),
+                                                        conj(a.u11*a.u22 - a.u12*a.u21))
+struct SU3fund{T}
+    t1::Complex{T}
+    t2::Complex{T}
+    t3::Complex{T}
+end
+Base.zero(::Type{SU3fund{T}}) where T <: AbstractFloat = SU3fund{T}(zero(T),zero(T),zero(T))
+Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU3fund{T}}) where T <: AbstractFloat = SU3fund{T}(complex(randn(rng,T),randn(rng,T)),
+                                                                                                      complex(randn(rng,T),randn(rng,T)),
+                                                                                                      complex(randn(rng,T),randn(rng,T)))
+
+

src/LatticeGPU.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    LatticeGPU.jl
+### created: Sat Jul 17 17:19:58 2021
+###                               
+
+
+module LatticeGPU
+
+include("Groups/Groups.jl")
+
+using .Groups
+export Group, Algebra, GMatrix
+export SU2, SU2alg, SU3, SU3alg, M3x3, M2x2, U1, U1alg, SU3fund
+export dot, expm, exp, dag, unitarize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat, dev_one
+
+include("Space/Space.jl")
+
+using .Space
+export SpaceParm
+export up, dw, updw, point_coord, point_index, point_time, point_color
+export BC_PERIODIC, BC_OPEN, BC_SF_AFWB, BC_SF_ORBI
+
+include("Fields/Fields.jl")
+using .Fields
+export vector_field, scalar_field, nscalar_field, scalar_field_point
+
+include("MD/MD.jl")
+using .MD
+export IntrScheme
+export omf4, leapfrog, omf2
+
+include("YM/YM.jl")
+
+using .YM
+export ztwist
+export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zero!, norm2
+export gauge_action, hamiltonian, plaquette, HMC!, OMF4!
+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
+
+include("Spinors/Spinors.jl")
+
+using .Spinors    
+export Spinor, Pgamma
+export imm, mimm
+export pmul, gpmul, gdagpmul, dmul
+
+include("Dirac/Dirac.jl")
+using .Dirac
+export DiracWorkspace
+export Dw!, DwdagDw!
+
+
+include("Solvers/Solvers.jl")
+using .Solvers
+export CG!
+
+end # module

src/MD/MD.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    MD.jl
+### created: Fri Oct  8 21:53:14 2021
+###                               
+
+module MD
+
+# Dalla Brida / Luscher coefficients of
+# OMF integrator
+const r1omf4 =  0.08398315262876693
+const r2omf4 =  0.25397851084105950
+const r3omf4 =  0.68223653357190910
+const r4omf4 = -0.03230286765269967
+const r5omf4 =  0.5-r1omf4-r3omf4
+const r6omf4 =  1.0-2.0*(r2omf4+r4omf4)
+
+const r1omf2 =  0.1931833275037836
+const r2omf2 =  0.5
+const r3omf2 =  1 - 2*r1omf2
+
+struct IntrScheme{N, T}
+    r::NTuple{N, T}
+    eps::T
+    ns::Int64    
+end
+
+
+omf2(::Type{T}, eps, ns) where T = IntrScheme{3,T}((r1omf2,r2omf2,r3omf2), eps, ns)
+omf4(::Type{T}, eps, ns) where T = IntrScheme{6,T}((r1omf4,r2omf4,r3omf4,r4omf4,r5omf4,r6omf4), eps, ns)
+leapfrog(::Type{T}, eps, ns) where T = IntrScheme{2,T}((0.5,1.0), eps, ns)
+
+
+import Base.show
+function Base.show(io::IO, int::IntrScheme{N,T}) where {N,T}
+
+    if N == 2
+        println(io, "LEAPFROG integration scheme")
+    elseif N == 3
+        println(io, "OMF2 integration scheme")
+    elseif N == 6
+        println(io, "OMF4 integration scheme")
+    end
+    println(io, "  - eps: ", int.eps)
+    println(io, "  - ns:  ", int.ns)
+    
+    return nothing
+end
+    
+export IntrScheme, omf4, leapfrog, omf2
+
+
+end

src/Solvers/CG.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    CG.jl
+### created: Tue Nov 30 11:10:57 2021
+###                               
+
+"""
+    function CG!
+
+Solves the linear equation `Ax = si`
+"""
+function CG!(si, U, m0, A, lp::SpaceParm, dws::DiracWorkspace)
+
+    dws.sr .= si
+    dws.sp .= si
+    norm = CUDA.mapreduce(x -> norm2(x), +, si)
+    fill!(si,zero(eltype(so)))
+    err = 0.0
+    
+    tol = eps * norm
+    iterations = 0
+
+    niter = 0
+    for i in 1:maxiter
+        A(dws.sAp, U, dws.sp, am0, dws.st, lp)
+        prod  = CUDA.mapreduce(x -> dot(x[1],x[2]), +, zip(dws.sp, dws.sAp))
+        alpha = norm/prod
+
+        si     .= si     .+ alpha .* dws.sp
+        dws.sr .= dws.sr .- alpha .* dws.sAp
+
+        err = CUDA.mapreduce(x -> norm2(x), +, dws.sr)
+        
+        if err < tol
+	    niter = i
+            break
+        end
+
+        beta = err/norm
+        dws.sp .= dws.sr .+ beta .* dws.sp
+        
+        norm = err;
+    end
+
+    if err > tol
+	error("CG! not converged after $maxiter iterations (Residuals: $err)")
+    end
+    
+    return niter
+end

src/Solvers/Solvers.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    Solvers.jl
+### created: Tue Nov 30 11:09:51 2021
+###                               
+
+module Solvers
+
+using CUDA, TimerOutputs
+using ..Space
+using ..Groups
+using ..Fields
+using ..YM
+using ..Spinors
+using ..Dirac
+
+include("CG.jl")
+export CG!
+
+
+
+
+end

src/Spinors/Spinors.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    Spinor.jl
+### created: Wed Nov 17 13:00:31 2021
+###                               
+
+module Spinors
+
+using ..Groups
+import ..Groups.imm, ..Groups.mimm, ..Groups.norm, ..Groups.norm2, ..Groups.dot
+
+struct Spinor{NS,G}
+    s::NTuple{NS,G}
+end
+Base.zero(::Type{Spinor{NS,G}}) where {NS,G} = Spinor{NS,G}(ntuple(i->zero(G), NS))
+
+"""
+    norm2(a::Spinor)
+
+Returns the norm squared of a fundamental element. Same result as `dot(a,a)`.
+"""
+@generated function norm2(a::Spinor{NS,G}) where {NS,G}
+
+    sum = :(norm2(a.s[1]))
+    for i in 2:NS
+        sum = :($sum + norm2(a.s[$i]))
+    end
+
+    return :($sum)
+end
+
+"""
+    norm(a::Spinor)
+
+Returns the norm of a fundamental element. Same result as sqrt(dot(a,a)).
+"""
+norm(a::Spinor{NS,G}) where {NS,G} = sqrt(norm2(a))
+
+"""
+    dot(a::Spinor,b::Spinor)
+
+Returns the scalar product of two spinors.
+"""
+@generated function dot(a::Spinor{NS,G},b::Spinor{NS,G}) where {NS,G}  
+
+    sum = :(dot(a.s[1],b.s[1]))
+    for i in 2:NS
+        sum = :($sum + norm2(a.s[$i]))
+    end
+
+    return :($sum)
+end
+
+
+"""
+    *(g::SU3{T},b::Spinor)
+
+Returns ga
+"""
+Base.:*(g::S,b::Spinor{NS,G}) where {S <: Group,NS,G} = Spinor{NS,G}(ntuple(i->g*b.s[i], NS))
+
+"""
+    \\(g::SU3{T},b::Spinor{NS,G})
+
+Returns g^+ a
+"""
+Base.:\(g::S,b::Spinor{NS,G}) where {S <: Group,NS,G} = Spinor{NS,G}(ntuple(i->g\b.s[i], NS))
+
+
+Base.:+(a::Spinor{NS,G},b::Spinor{NS,G}) where {NS,G} = Spinor{NS,G}(ntuple(i->a.s[i]+b.s[i], NS))
+Base.:-(a::Spinor{NS,G},b::Spinor{NS,G}) where {NS,G} = Spinor{NS,G}(ntuple(i->a.s[i]-b.s[i], NS))
+Base.:+(a::Spinor{NS,G})                 where {NS,G} = a
+Base.:-(a::Spinor{NS,G})                 where {NS,G} = Spinor{NS,G}(ntuple(i->-b.s[i], NS))
+imm(a::Spinor{NS,G})                     where {NS,G} = Spinor{NS,G}(ntuple(i->imm(b.s[i]), NS))
+mimm(a::Spinor{NS,G})                    where {NS,G} = Spinor{NS,G}(ntuple(i->mimm(b.s[i]), NS))
+
+
+# Operations with numbers
+Base.:*(a::Spinor{NS,G},b::Number) where {NS,G} = Spinor{NS,G}(ntuple(i->b*a.s[i], NS))
+Base.:*(b::Number,a::Spinor{NS,G}) where {NS,G} = Spinor{NS,G}(ntuple(i->b*a.s[i], NS))
+Base.:/(a::Spinor{NS,G},b::Number) where {NS,G} = Spinor{NS,G}(ntuple(i->a.s[i]/b, NS))
+
+
+## 
+# Operations specific for Spinors in 4D
+##
+
+# dummy structs for dispatch:
+# Projectors (1+s\gamma_N)
+struct Pgamma{N,S}
+end
+
+"""
+    pmul(Pgamma{N,S}, a::Spinor)
+
+Returns ``(1+s\\gamma_N)a``.
+"""
+@inline function pmul(::Type{Pgamma{4,1}},  a::Spinor{4,G}) where {NS,G}
+    
+    r1 = a.s[1]+a.s[3]
+    r2 = a.s[2]+a.s[4]
+    return Spinor{4,G}((r1,r2,r1,r2))
+end
+@inline function pmul(::Type{Pgamma{4,-1}}, a::Spinor{4,G}) where {NS,G}
+    
+    r1 = a.s[1]-a.s[3]
+    r2 = a.s[2]-a.s[4] 
+    return Spinor{4,G}((r1,r2,-r1,-r2))
+end
+    
+@inline function pmul(::Type{Pgamma{1,1}},  a::Spinor{4,G}) where {NS,G}
+
+    r1 = a.s[1]+imm(a.s[4])
+    r2 = a.s[2]+imm(a.s[3])
+    return Spinor{4,G}((r1,r2,mimm(r2),mimm(r1)))
+end
+@inline function pmul(::Type{Pgamma{1,-1}}, a::Spinor{4,G}) where {NS,G}
+
+    r1 = a.s[1]-imm(a.s[4])
+    r2 = a.s[2]-imm(a.s[3])
+    return Spinor{4,G}((r1,r2,imm(r2),imm(r1)))
+end
+
+@inline function pmul(::Type{Pgamma{2,1}},  a::Spinor{4,G}) where {NS,G}
+
+    r1 = a.s[1]+a.s[4]
+    r2 = a.s[2]-a.s[3]
+    return Spinor{4,G}((r1,r2,-r2,r1))
+end
+@inline function pmul(::Type{Pgamma{2,-1}}, a::Spinor{4,G}) where {NS,G}
+
+    r1 = a.s[1]-a.s[4]
+    r2 = a.s[2]+a.s[3]
+    return Spinor{4,G}((r1,r2,r2,-r1))
+end
+
+@inline function pmul(::Type{Pgamma{3,1}},  a::Spinor{4,G}) where {NS,G}
+
+    r1 = a.s[1]+imm(a.s[3])
+    r2 = a.s[2]-imm(a.s[4])
+    return Spinor{4,G}((r1,r2,mimm(r1),imm(r2)))
+end
+@inline function pmul(::Type{Pgamma{3,-1}},  a::Spinor{4,G}) where {NS,G}
+
+    r1 = a.s[1]-imm(a.s[3])
+    r2 = a.s[2]+imm(a.s[4])
+    return Spinor{4,G}((r1,r2,imm(r1),mimm(r2)))
+end
+
+
+"""
+    gpmul(pgamma{N,S}, g::G, a::Spinor) G <: Group
+
+Returns ``g(1+s\\gamma_N)a``
+"""
+@inline function gpmul(::Type{Pgamma{4,1}},  g, a::Spinor{4,G}) where {NS,G}
+    
+    r1 = g*(a.s[1]+a.s[3])
+    r2 = g*(a.s[2]+a.s[4])
+    return Spinor{4,G}((r1,r2,r1,r2))
+end
+@inline function gpmul(::Type{Pgamma{4,-1}}, g, a::Spinor{4,G}) where {NS,G}
+    
+    r1 = g*(a.s[1]-a.s[3])
+    r2 = g*(a.s[2]-a.s[4]) 
+    return Spinor{4,G}((r1,r2,-r1,-r2))
+end
+    
+@inline function gpmul(::Type{Pgamma{1,1}},  g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g*(a.s[1]+imm(a.s[4]))
+    r2 = g*(a.s[2]+imm(a.s[3]))
+    return Spinor{4,G}((r1,r2,mimm(r2),mimm(r1)))
+end
+@inline function gpmul(::Type{Pgamma{1,-1}}, g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g*(a.s[1]-imm(a.s[4]))
+    r2 = g*(a.s[2]-imm(a.s[3]))
+    return Spinor{4,G}((r1,r2,imm(r2),imm(r1)))
+end
+
+@inline function gpmul(::Type{Pgamma{2,1}},  g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g*(a.s[1]+a.s[4])
+    r2 = g*(a.s[2]-a.s[3])
+    return Spinor{4,G}((r1,r2,-r2,r1))
+end
+@inline function gpmul(::Type{Pgamma{2,-1}}, g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g*(a.s[1]-a.s[4])
+    r2 = g*(a.s[2]+a.s[3])
+    return Spinor{4,G}((r1,r2,r2,-r1))
+end
+
+@inline function gpmul(::Type{Pgamma{3,1}},  g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g*(a.s[1]+imm(a.s[3]))
+    r2 = g*(a.s[2]-imm(a.s[4]))
+    return Spinor{4,G}((r1,r2,mimm(r1),imm(r2)))
+end
+@inline function gpmul(::Type{Pgamma{3,-1}},  g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g*(a.s[1]-imm(a.s[3]))
+    r2 = g*(a.s[2]+imm(a.s[4]))
+    return Spinor{4,G}((r1,r2,imm(r1),mimm(r2)))
+end
+
+"""
+    gdagpmul(pgamma{N,S}, g::G, a::Spinor) G <: Group
+
+Returns ``g^+ (1+s\\gamma_N)a``
+"""
+@inline function gdagpmul(::Type{Pgamma{4,1}},  g, a::Spinor{4,G}) where {NS,G}
+    
+    r1 = g\(a.s[1]+a.s[3])
+    r2 = g\(a.s[2]+a.s[4])
+    return Spinor{4,G}((r1,r2,r1,r2))
+end
+@inline function gdagpmul(::Type{Pgamma{4,-1}}, g, a::Spinor{4,G}) where {NS,G}
+    
+    r1 = g\(a.s[1]-a.s[3])
+    r2 = g\(a.s[2]-a.s[4]) 
+    return Spinor{4,G}((r1,r2,-r1,-r2))
+end
+    
+@inline function gdagpmul(::Type{Pgamma{1,1}},  g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g\(a.s[1]+imm(a.s[4]))
+    r2 = g\(a.s[2]+imm(a.s[3]))
+    return Spinor{4,G}((r1,r2,mimm(r2),mimm(r1)))
+end
+@inline function gdagpmul(::Type{Pgamma{1,-1}}, g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g\(a.s[1]-imm(a.s[4]))
+    r2 = g\(a.s[2]-imm(a.s[3]))
+    return Spinor{4,G}((r1,r2,imm(r2),imm(r1)))
+end
+
+@inline function gdagpmul(::Type{Pgamma{2,1}},  g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g\(a.s[1]+a.s[4])
+    r2 = g\(a.s[2]-a.s[3])
+    return Spinor{4,G}((r1,r2,-r2,r1))
+end
+@inline function gdagpmul(::Type{Pgamma{2,-1}}, g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g\(a.s[1]-a.s[4])
+    r2 = g\(a.s[2]+a.s[3])
+    return Spinor{4,G}((r1,r2,r2,-r1))
+end
+
+@inline function gdagpmul(::Type{Pgamma{3,1}},  g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g\(a.s[1]+imm(a.s[3]))
+    r2 = g\(a.s[2]-imm(a.s[4]))
+    return Spinor{4,G}((r1,r2,mimm(r1),imm(r2)))
+end
+@inline function gdagpmul(::Type{Pgamma{3,-1}},  g, a::Spinor{4,G}) where {NS,G}
+
+    r1 = g\(a.s[1]-imm(a.s[3]))
+    r2 = g\(a.s[2]+imm(a.s[4]))
+    return Spinor{4,G}((r1,r2,imm(r1),mimm(r2)))
+end
+
+
+# dummy structs for dispatch:
+# Basis of \\Gamma_n
+struct Gamma{N}
+end
+
+"""
+    dmul(n::Int64, a::Spinor)
+
+Returns ``\\Gamma_n a``
+
+indexing for Dirac basis ``\\Gamma_n``:
+
+ 1  gamma1
+ 2  gamma2
+ 3  gamma3
+ 4  gamma0
+ 5  gamma5
+ 6  gamma1 gamma5
+ 7  gamma2 gamma5
+ 8  gamma3 gamma5
+ 9  gamma0 gamma5
+10  sigma01
+11  sigma02
+12  sigma03
+13  sigma12
+14  sigma23
+15  sigma31
+16  identity
+
+"""
+@inline dmul(::Type{Gamma{1}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(mimm(a.s[4]), mimm(a.s[3]), imm(a.s[2]), imm(a.s[1]))
+@inline dmul(::Type{Gamma{2}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[4], a.s[3], a.s[2], -a.s[1])
+@inline dmul(::Type{Gamma{3}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(mimm(a.s[3]), imm(a.s[4]), imm(a.s[1]), mimm(a.s[2]))
+@inline dmul(::Type{Gamma{4}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[3], -a.s[4], -a.s[1], -a.s[2])
+@inline dmul(::Type{Gamma{5}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[1],  a.s[2],  -a.s[3], -a.s[4])
+@inline dmul(::Type{Gamma{6}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( imm(a.s[4]),  imm(a.s[3]),   imm(a.s[2]),  imm(a.s[1]))
+@inline dmul(::Type{Gamma{7}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[4], -a.s[3],   a.s[2], -a.s[1])
+@inline dmul(::Type{Gamma{8}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( imm(a.s[3]), mimm(a.s[4]),   imm(a.s[1]), mimm(a.s[2]))
+@inline dmul(::Type{Gamma{9}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[3],  a.s[4],  -a.s[1], -a.s[2])
+@inline dmul(::Type{Gamma{10}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[2],  a.s[1],  -a.s[4], -a.s[3])
+@inline dmul(::Type{Gamma{11}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(mimm(a.s[2]),  imm(a.s[1]),   imm(a.s[4]), mimm(a.s[3]))
+@inline dmul(::Type{Gamma{12}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[1], -a.s[2],  -a.s[3],  a.s[4])
+@inline dmul(::Type{Gamma{13}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[1],  a.s[2],  -a.s[3],  a.s[4])
+@inline dmul(::Type{Gamma{14}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[2], -a.s[1],  -a.s[4], -a.s[3])
+@inline dmul(::Type{Gamma{15}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( imm(a.s[2]), mimm(a.s[1]),   imm(a.s[4]), mimm(a.s[3]))
+@inline dmul(::Type{Gamma{16}}, a::Spinor{4,G}) where {G} = a
+
+export Spinor, Pgamma
+export norm, norm2, dot, imm, mimm
+export pmul, gpmul, gdagpmul, dmul
+
+end

src/YM/YM.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    YM.jl
+### created: Mon Jul 12 16:23:51 2021
+###                               
+
+
+module YM
+
+using CUDA, Random, TimerOutputs, BDIO
+using ..Space
+using ..Groups
+using ..Fields
+using ..MD
+
+import Base.show
+
+struct GaugeParm{T,G,N}
+    beta::T
+    c0::T
+    cG::NTuple{2,T}
+    ng::Int64
+
+    Ubnd::NTuple{N, G}
+
+    GaugeParm{T1,T2,T3}(a,b,c,d,e) where {T1,T2,T3} = new{T1,T2,T3}(a,b,c,d,e)
+    function GaugeParm{T}(::Type{G}, bt, c0, cG, phi, iL) where {T,G}
+
+        degree(::Type{SU2{T}}) where T <: AbstractFloat = 2
+        degree(::Type{SU3{T}}) where T <: AbstractFloat = 3
+        ng = degree(G)
+        nsd = length(iL)
+        
+        return new{T,G,nsd}(bt, c0, cG, ng, ntuple(id->bndfield(phi[1], phi[2], iL[id]), nsd))
+    end
+    function GaugeParm{T}(::Type{G}, bt, c0) where {T,G}
+        
+        degree(::Type{SU2{T}}) where T <: AbstractFloat = 2
+        degree(::Type{SU3{T}}) where T <: AbstractFloat = 3
+        ng = degree(G)
+
+        return new{T,G,0}(bt, c0, (0.0,0.0), ng, ())
+    end
+end
+export GaugeParm
+function Base.show(io::IO, gp::GaugeParm{T, G, N}) where {T,G,N}
+
+    println(io, "Group:  ", G)
+    println(io, " - beta:              ", gp.beta)
+    println(io, " - c0:                ", gp.c0)
+    println(io, " - cG:                ", gp.cG)
+    if (N > 0)
+        for i in 1:N
+            println(io, "   - Boundary link:     ", gp.Ubnd[i])
+        end
+    end
+    
+    return nothing
+end
+
+struct YMworkspace{T}
+    GRP
+    ALG
+    PRC
+    frc1
+    frc2
+    mom
+    U1
+    cm # complex of volume
+    rm # float   of volume
+    function YMworkspace(::Type{G}, ::Type{T}, lp::SpaceParm) where {G <: Group, T <: AbstractFloat}
+        
+        @timeit "Allocating YMWorkspace" begin
+            if (G == SU2)
+                GRP = SU2
+                ALG = SU2alg
+                f1 = vector_field(SU2alg{T}, lp)
+                f2 = vector_field(SU2alg{T}, lp)
+                mm = vector_field(SU2alg{T}, lp)
+                u1 = vector_field(SU2{T},    lp)
+            end
+            
+            if (G == SU3)
+                GRP = SU3
+                ALG = SU3alg
+                f1 = vector_field(SU3alg{T}, lp)
+                f2 = vector_field(SU3alg{T}, lp)
+                mm = vector_field(SU3alg{T}, lp)
+                u1 = vector_field(SU3{T},    lp)
+            end
+            cs = scalar_field_point(Complex{T}, lp)
+            rs = scalar_field_point(T, lp)
+        end
+            
+        return new{T}(GRP,ALG,T,f1, f2, mm, u1, cs, rs)
+    end
+end
+export YMworkspace
+function Base.show(io::IO, ymws::YMworkspace)
+    
+    println(io, "Workspace for Group:   ", ymws.GRP)
+    println(io, "              Algebra: ", ymws.ALG)
+    println(io, "Precision:             ", ymws.PRC)
+
+    return nothing
+end
+
+
+function ztwist(gp::GaugeParm{T,G}, lp::SpaceParm{N,M,B,D}) where {T,G,N,M,B,D}
+
+    function plnf(ipl)
+        id1, id2 = lp.plidx[ipl]
+        return convert(Complex{T},exp(2im * pi * lp.ntw[ipl]/(gp.ng)))
+    end
+
+    return ntuple(i->plnf(i), M)
+end
+
+function ztwist(gp::GaugeParm{T,G}, lp::SpaceParm{N,M,B,D}, ipl::Int) where {T,G,N,M,B,D}
+
+    id1, id2 = lp.plidx[ipl]
+    return convert(Complex{T},exp(2im * pi * lp.ntw[ipl]/(gp.ng)))
+end
+export ztwist
+
+
+include("YMfields.jl")
+export randomize!, zero!, norm2
+
+include("YMact.jl")
+export krnl_plaq!, force0_wilson!
+
+include("YMhmc.jl")
+export gauge_action, hamiltonian, plaquette, HMC!, OMF4!
+
+include("YMflow.jl")
+export FlowIntr, flw, flw_adapt
+export Eoft_clover, Eoft_plaq, Qtop
+export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3
+
+include("YMsf.jl")
+export sfcoupling, bndfield, setbndfield
+
+include("YMio.jl")
+export import_lex64, import_cern64, import_bsfqcd, save_cnfg, read_cnfg
+
+end

src/YM/YMfields.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    YMfields.jl
+### created: Thu Jul 15 15:16:47 2021
+###                               
+
+function randomize!(f, lp::SpaceParm, ymws::YMworkspace) 
+        
+    if ymws.ALG == SU2alg
+        @timeit "Randomize SU(2) algebra field" begin
+            m = CUDA.randn(ymws.PRC, lp.bsz,lp.ndim,3,lp.rsz)
+            CUDA.@sync begin
+                CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_SU2!(f,m,lp)
+            end
+        end
+        return nothing
+    end
+
+    if ymws.ALG == SU3alg
+        @timeit "Randomize SU(3) algebra field" begin
+            m = CUDA.randn(ymws.PRC, lp.bsz,lp.ndim,8,lp.rsz)
+            CUDA.@sync begin
+                CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_SU3!(f,m,lp)
+            end
+        end
+        return nothing
+    end
+
+    return nothing
+end
+
+function krnl_assign_SU3!(frc::AbstractArray{T}, m, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {T,N,M,D}
+
+    @inbounds begin
+        b = Int64(CUDA.threadIdx().x)
+        r = Int64(CUDA.blockIdx().x)
+        for id in 1:lp.ndim
+            frc[b,id,r] = SU3alg(m[b,id,1,r], m[b,id,2,r], m[b,id,3,r],
+                                 m[b,id,4,r], m[b,id,5,r], m[b,id,6,r],
+                                 m[b,id,7,r], m[b,id,8,r])
+        end
+    end
+    
+    return nothing
+end
+
+function krnl_assign_SU3!(frc::AbstractArray{T}, m, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D}
+
+    @inbounds begin
+        b = Int64(CUDA.threadIdx().x)
+        r = Int64(CUDA.blockIdx().x)
+        it = point_time((b,r), lp)
+
+        if ((B==BC_SF_AFWB)||(B==BC_SF_ORBI))
+            if it == 1
+                for id in 1:lp.ndim-1
+                    frc[b,id,r] = zero(T)
+                end
+                frc[b,N,r] = SU3alg(m[b,N,1,r], m[b,N,2,r], m[b,N,3,r],
+                                    m[b,N,4,r], m[b,N,5,r], m[b,N,6,r],
+                                    m[b,N,7,r], m[b,N,8,r])
+            else
+                for id in 1:lp.ndim
+                    frc[b,id,r] = SU3alg(m[b,id,1,r], m[b,id,2,r], m[b,id,3,r],
+                                         m[b,id,4,r], m[b,id,5,r], m[b,id,6,r],
+                                         m[b,id,7,r], m[b,id,8,r])
+                end
+            end
+        end
+    end
+        
+    return nothing
+end
+
+function krnl_assign_SU2!(frc, m, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {N,M,D}
+
+    @inbounds begin
+        b = Int64(CUDA.threadIdx().x)
+        r = Int64(CUDA.blockIdx().x)
+        for id in 1:lp.ndim
+            frc[b,id,r] = SU2alg(m[b,id,1,r], m[b,id,2,r], m[b,id,3,r])
+        end
+    end
+    
+    return nothing
+end

src/YM/YMhmc.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    YMhmc.jl
+### created: Thu Jul 15 11:27:28 2021
+###                               
+
+"""
+    
+    function gauge_action(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace)
+
+Returns the value of the gauge plaquette action for the configuration U. The parameters `\beta` and `c0` are taken from the `gp` structure. 
+"""
+function gauge_action(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}) where T <: AbstractFloat
+
+    ztw = ztwist(gp, lp)
+    if abs(gp.c0-1) < 1.0E-10
+        @timeit "Wilson gauge action" begin
+            CUDA.@sync begin
+                CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, gp.Ubnd, gp.cG[1], ztw, lp)
+            end
+        end
+    else
+        @timeit "Improved gauge action" begin
+            CUDA.@sync begin
+                CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_impr!(ymws.cm, U, gp.c0, (1-gp.c0)/8, gp.Ubnd, gp.cG[1], ztw, lp)
+            end
+        end
+    end
+    S = gp.beta*( prod(lp.iL)*lp.npls*(gp.c0 + (1-gp.c0)/8) -
+                  CUDA.mapreduce(real, +, ymws.cm)/gp.ng )
+
+    return S
+end
+
+function plaquette(U, lp::SpaceParm{N,M,B,D}, gp::GaugeParm, ymws::YMworkspace) where {N,M,B,D}
+
+    ztw = ztwist(gp, lp)
+    @timeit "Plaquette measurement" begin
+        CUDA.@sync begin
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, gp.Ubnd, one(gp.cG[1]), ztw, lp)
+        end
+    end
+
+    return CUDA.mapreduce(real, +, ymws.cm)/(prod(lp.iL)*lp.npls)
+end
+    
+function hamiltonian(mom, U, lp, gp, ymws)
+    @timeit "Computing Hamiltonian" begin
+        K = CUDA.mapreduce(norm2, +, mom)/2 
+        V = gauge_action(U, lp, gp, ymws)
+    end
+        
+    return K+V
+end
+
+function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false) where T
+
+    @timeit "HMC trayectory" begin
+        
+        ymws.U1 .= U
+        
+        randomize!(ymws.mom, lp, ymws)
+        hini = hamiltonian(ymws.mom, U, lp, gp, ymws)
+        
+        MD!(ymws.mom, U, int, lp, gp, ymws)
+        
+        dh   = hamiltonian(ymws.mom, U, lp, gp, ymws) - hini
+        pacc = exp(-dh)
+        
+        acc = true
+        if (noacc)
+            return dh, acc
+        end
+        
+        if (pacc < 1.0)
+            r = rand()
+            if (pacc < r) 
+                U .= ymws.U1
+                acc = false
+            end
+        end
+
+        U .= unitarize.(U)
+        
+    end
+    return dh, acc
+end
+HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false) where T = HMC!(U, omf4(T, eps, ns), lp, gp, ymws; noacc=noacc)
+
+function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}) where {NI, T <: AbstractFloat}
+    
+    @timeit "MD evolution" begin
+
+        ee = int.eps*gp.beta/gp.ng
+        force_gauge(ymws, U, gp.c0, gp, lp)
+        mom .= mom .+ (int.r[1]*ee) .* ymws.frc1
+        for i in 1:int.ns
+            k   = 2
+            off = 1
+            for j in 1:NI-1
+                U .= expm.(U, mom, int.eps*int.r[k])
+                if k == NI
+                    off = -1
+                end
+                k += off
+                
+                force_gauge(ymws, U, gp.c0, gp, lp)
+                if (i < int.ns) && (k == 1)
+                    mom .= mom .+ (2*int.r[k]*ee) .* ymws.frc1
+                else
+                    mom .= mom .+ (int.r[k]*ee) .* ymws.frc1
+                end
+                if k == NI
+                    off = -1
+                end
+                k += off
+            end
+        end
+    end
+    
+    return nothing
+end

src/YM/YMio.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    YMio.jl
+### created: Wed Nov 10 12:58:27 2021
+###                               
+
+"""
+    read_cnfg(U, lp::SpaceParm{4,M,B,D},)
+
+Reads configuration from file `fname` using the native (BDIO) format.
+"""
+function read_cnfg(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])
+    
+    dtr = (2,3,4,1)
+    assign(id, V, i3) = SU3{Float64}(V[1,dtr[id],i3],V[2,dtr[id],i3],V[3,dtr[id],i3],
+                                     V[4,dtr[id],i3],V[5,dtr[id],i3],V[6,dtr[id],i3])
+
+    while BDIO_get_uinfo(fb) != 8
+        BDIO_seek!(fb)
+    end
+    Ucpu = Array{SU3{Float64}, 3}(undef, lp.bsz, lp.ndim, lp.rsz)
+    V = Array{ComplexF64, 3}(undef, 9, lp.ndim, lp.iL[3])
+    for i4 in 1:lp.iL[4]
+        for i1 in 1:lp.iL[1]
+            for i2 in 1:lp.iL[2]
+                BDIO_read(fb, vec(V))
+                for i3 in 1:lp.iL[3]
+                    b, r = point_index(CartesianIndex(i1,i2,i3,i4), lp)
+                    for id in 1:lp.ndim
+                        Ucpu[b,id,r] = assign(id, V, i3)
+                    end
+                end
+            end
+        end
+    end
+
+    if ibc == BC_SF_AFWB || ibc == BC_SF_ORBI
+        BDIO_read(fb, V)
+        Ubnd = ntuple(i->assign(i, V, 1), 3)
+        BDIO_close!(fb)
+
+        return CuArray(Ucpu), Ubnd
+    else
+        BDIO_close!(fb)
+        return CuArray(Ucpu)
+    end
+end
+
+
+"""
+    save_cnfg(fname, U, lp::SpaceParm, gp::GaugeParm; run::Union{Nothing,String}=nothing)
+
+Saves configuration `U` in the file `fname` using the native (BDIO) format.
+"""
+function save_cnfg(fname::String, U, lp::SpaceParm{4,M,B,D}, gp::GaugeParm; run::Union{Nothing,String}=nothing) where {M,B,D}
+
+    ihdr = [convert(Int32, 1653996111),convert(Int32,2)]
+    UID_HDR = 14
+
+    if isfile(fname)
+        fb = BDIO_open(fname, "a")
+    else
+        fb = BDIO_open(fname, "w", "Configuration of LatticeGPU.jl")
+        BDIO_start_record!(fb, BDIO_BIN_GENERIC, UID_HDR)
+        BDIO_write!(fb, ihdr)
+        if run != nothing
+            BDIO_write!(fb, run*"\0")
+        end
+        BDIO_write_hash!(fb)
+        
+        dfoo = Vector{Float64}(undef, 16)
+        BDIO_start_record!(fb, BDIO_BIN_GENERIC, 1)
+        BDIO_write!(fb, [convert(Int32, 4)])
+        BDIO_write!(fb, [convert(Int32, B)])
+        BDIO_write!(fb, [convert(Int32, gp.ng)])
+        BDIO_write!(fb, [convert(Int32, 0)])
+        BDIO_write!(fb, [convert(Int32, lp.iL[i]) for i in 1:4])
+        BDIO_write!(fb, [convert(Int32, lp.ntw[i]) for i in 1:M])
+        BDIO_write!(fb, [gp.beta, gp.c0, gp.cG[1], gp.cG[2]])
+    end
+    BDIO_write_hash!(fb)
+    
+    dtr = (2,3,4,1)
+
+    function assign!(id, V, i3, M::M3x3{T}) where T
+
+        V[1,dtr[id],i3] = M.u11
+        V[2,dtr[id],i3] = M.u12
+        V[3,dtr[id],i3] = M.u13
+        V[4,dtr[id],i3] = M.u21
+        V[5,dtr[id],i3] = M.u22
+        V[6,dtr[id],i3] = M.u23
+        V[7,dtr[id],i3] = M.u31
+        V[8,dtr[id],i3] = M.u32
+        V[9,dtr[id],i3] = M.u33
+        return nothing
+    end
+    assign!(id, V, i3, g::SU3{T}) where T = assign!(id,V,i3,convert(M3x3{T}, g))
+
+    BDIO_start_record!(fb, BDIO_BIN_F64LE, 8, true)
+    Ucpu = Array(U)
+    V = Array{ComplexF64, 3}(undef, 9, lp.ndim, lp.iL[3])
+    for i4 in 1:lp.iL[4]
+        for i1 in 1:lp.iL[1]
+            for i2 in 1:lp.iL[2]
+                for i3 in 1:lp.iL[3]
+                    b, r = point_index(CartesianIndex(i1,i2,i3,i4), lp)
+                    for id in 1:lp.ndim
+                        assign!(id, V, i3, Ucpu[b,id,r])
+                    end
+                end
+                BDIO_write!(fb, V)
+            end
+        end
+    end
+    
+    if B == BC_SF_AFWB || B == BC_SF_ORBI
+        for i3 in 1:lp.iL[3]
+            for id in 1:lp.ndim-1
+                assign!(id, V, i3, Ubnd[id])
+            end
+            assign!(4, V, i3, zero(M3x3{Float64}))
+        end
+        for i1 in 1:lp.iL[1]
+            for i2 in 1:lp.iL[2]
+                BDIO_write!(fb, V)
+            end
+        end
+    end        
+    BDIO_write_hash!(fb)
+    BDIO_close!(fb)
+    
+    return nothing
+end
+
+
+"""
+    function import_bsfqcd(fname::String, lp::SpaceParm)
+
+import a double precision configuration in bsfqcd format. SF boundary conditions are assummed. 
+"""
+function import_bsfqcd(fname, lp::SpaceParm)
+
+    fp = BDIO_open(fname, "r")
+    while BDIO_seek!(fp)
+        if (BDIO_get_uinfo(fb) == 2)
+            break
+        end
+    end
+
+    
+    dtr = [2,3,4,1]
+
+    assign(id, V, i3) = SU3{Float64}(V[1,dtr[id],i3],V[2,dtr[id],i3],V[3,dtr[id],i3],
+                                     V[4,dtr[id],i3],V[5,dtr[id],i3],V[6,dtr[id],i3])
+    
+    Ucpu = Array{SU3{Float64}, 3}(undef, lp.bsz, lp.ndim, lp.rsz)
+    V = Array{ComplexF64, 3}(undef, 9, lp.ndim, lp.iL[3])
+    for i4 in 1:lp.iL[4]
+        for i1 in 1:lp.iL[1]
+            for i2 in 1:lp.iL[2]
+                BDIO_read(fp, V)
+                for i3 in 1:lp.iL[3]
+                    b, r = point_index(CartesianIndex(i1,i2,i3,i4), lp)
+                    for id in 1:lp.ndim
+                        Ucpu[b,id,r] = assign(id, V, i3)
+                    end
+                end
+            end
+        end
+    end
+
+    BDIO_read(fp, V)
+    Ubnd = ntuple(i->assign(i, V, 1), 3)
+    close(fp)
+
+    return CuArray(Ucpu), Ubnd
+end
+
+"""
+    function import_lex64(fname::String, lp::SpaceParm)
+
+import a double precision configuration in lexicographic format. SF boundary conditions are assummed. 
+"""
+function import_lex64(fname, lp::SpaceParm)
+
+    fp = open(fname, "r")
+
+    dtr = [2,3,4,1]
+
+    assign(id, V, i3) = SU3{Float64}(V[1,dtr[id],i3],V[2,dtr[id],i3],V[3,dtr[id],i3],
+                                     V[4,dtr[id],i3],V[5,dtr[id],i3],V[6,dtr[id],i3])
+    
+    Ucpu = Array{SU3{Float64}, 3}(undef, lp.bsz, lp.ndim, lp.rsz)
+    V = Array{ComplexF64, 3}(undef, 9, lp.ndim, lp.iL[3])
+    for i4 in 1:lp.iL[4]
+        for i1 in 1:lp.iL[1]
+            for i2 in 1:lp.iL[2]
+                read!(fp, V)
+                for i3 in 1:lp.iL[3]
+                    b, r = point_index(CartesianIndex(i1,i2,i3,i4), lp)
+                    for id in 1:lp.ndim
+                        Ucpu[b,id,r] = assign(id, V, i3)
+                    end
+                end
+            end
+        end
+    end
+
+    read!(fp, V)
+    Ubnd = ntuple(i->assign(i, V, 1), 3)
+    close(fp)
+
+    return CuArray(Ucpu), Ubnd
+end
+
+"""
+    function import_cern64(fname::String, ibc, lp::SpaceParm)
+
+import a double precision configuration in CERN format. 
+"""
+function import_cern64(fname, ibc, lp::SpaceParm; log=true)
+
+    fp = open(fname, "r")
+    iL = Vector{Int32}(undef, 4)
+    read!(fp, iL)
+    avgpl = Vector{Float64}(undef, 1)
+    read!(fp, avgpl)
+    if log
+        println("# [import_cern64] Read from conf file: ", iL, " (plaq: ", avgpl, ")")
+    end
+    
+    dtr = [4,1,2,3]
+    assign(V, ic) = SU3{Float64}(V[1,ic],V[2,ic],V[3,ic],V[4,ic],V[5,ic],V[6,ic])
+    
+    Ucpu = Array{SU3{Float64}, 3}(undef, lp.bsz, lp.ndim, lp.rsz)
+    V = Array{ComplexF64, 2}(undef, 9, 2)
+    for i4 in 1:lp.iL[4]
+        for i1 in 1:lp.iL[1]
+            for i2 in 1:lp.iL[2]
+                for i3 in 1:lp.iL[3]
+                    if (mod(i1+i2+i3+i4-4, 2) == 1)
+                        b, r = point_index(CartesianIndex(i1,i2,i3,i4), lp)
+                        for id in 1:lp.ndim
+                            read!(fp, V)
+                            Ucpu[b,dtr[id],r] = assign(V, 1)
+
+                            bd, rd = dw((b,r), dtr[id], lp)
+                            Ucpu[bd,dtr[id],rd] = assign(V, 2)
+                        end
+                    end
+                end
+            end
+        end
+    end
+    close(fp)
+
+    return CuArray(Ucpu)
+end

src/YM/YMsf.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    YMsf.jl
+### created: Tue Oct 26 14:50:55 2021
+###                               
+
+""" 
+    sfcoupling(U, lp::SpaceParm{N,M,B,D}, gp::GaugeParm, ymws::YMworkspace) where {N,M,B,D}
+
+Measures the Schrodinger Functional coupling `ds/d\eta` and `d^2S/d\eta d\nu`. 
+"""
+function sfcoupling(U, lp::SpaceParm{N,M,B,D}, gp::GaugeParm, ymws::YMworkspace) where {N,M,B,D}
+
+    if lp.iL[end] < 4
+        error("Array too small to store partial sums")
+    end
+
+    if !((B==BC_SF_AFWB) || (B==BC_SF_ORBI))
+        error("SF coupling can only be measured with SF boundary conditions")
+    end
+    
+    @timeit "SF coupling measurement" begin
+        T = eltype(ymws.rm)
+        tmp = zeros(T,lp.iL[end])
+        fill!(ymws.rm, zero(T))
+        CUDA.@sync begin
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfcoupling!(ymws.rm, U, gp.Ubnd, lp)
+        end
+
+        tp = ntuple(i->i, N-1)
+        tmp .= reshape(Array(CUDA.reduce(+, ymws.rm;dims=tp)),lp.iL[end])
+        
+        dsdeta = (gp.cG[1]*gp.beta/(2*gp.ng))*(tmp[1] + tmp[end])
+        ddnu   = (gp.cG[1]*gp.beta/(2*gp.ng))*(tmp[2] + tmp[end-1])
+    end
+    
+    return dsdeta, ddnu
+end
+    
+function krnl_sfcoupling!(rm, U::AbstractArray{T}, Ubnd, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D}
+
+    @inbounds begin
+        b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
+        I    = point_coord((b,r), lp)
+        it   = I[N]
+
+        SR3::eltype(rm)   = 1.73205080756887729352744634151
+        SR3x2::eltype(rm) = 3.46410161513775458705489268302
+        
+        if (it == 1)
+            but, rut = up((b,r), N, lp)
+            IU = point_coord((but,rut), lp)
+            for id in 1:N-1
+                bu, ru = up((b,r), id, lp)
+
+                X = projalg(U[b,id,r]*U[bu,N,ru]/(U[b,N,r]*U[but,id,rut]))
+                rm[I]  += (3*X.t7 + SR3   * X.t8)/lp.iL[id]
+                rm[IU] += (2*X.t7 - SR3x2 * X.t8)/lp.iL[id]
+            end
+        elseif (it == lp.iL[end])
+            bdt, rdt = dw((b,r), N, lp)
+            ID = point_coord((bdt,rdt), lp)
+            for id in 1:N-1
+                bu, ru = up((b,r), id, lp)
+
+                X = projalg(Ubnd[id]/(U[b,id,r]*U[bu,N,ru])*U[b,N,r])
+                rm[I]  -= (3*X.t7 + SR3   * X.t8)/lp.iL[id]
+                rm[ID] += (2*X.t7 - SR3x2 * X.t8)/lp.iL[id]
+            end
+        end
+    end
+        
+    return nothing
+end
+
+
+@inline function bndfield(phi1::T, phi2::T, iL) where T <: AbstractFloat
+
+    SR3::T = 1.73205080756887729352744634151
+
+    zt = zero(T)
+    X = SU3alg{T}(zt,zt,zt,zt,zt,zt,(phi1-phi2)/iL,SR3*(phi1+phi2)/iL)
+    
+    return exp(X)
+end
+    
+
+function setbndfield(U, phi, lp::SpaceParm{N,M,B,D}) where {N,M,B,D}
+
+    CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_setbnd_it0!(U, phi[1], phi[2], lp)
+    end
+    
+    return nothing
+end
+
+function krnl_setbnd_it0!(U, phi1, phi2, lp::SpaceParm{N,M,B,D})  where {N,M,B,D}
+
+    @inbounds begin
+        b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
+        it   = point_time((b,r), lp)
+
+        SFBC = (B == BC_SF_AFWB) || (B == BC_SF_ORBI)
+        
+        if (it == 1) && SFBC
+            for id in 1:N-1
+                U[b,id,r] = bndfield(phi1,phi2,lp.iL[id])
+            end
+        end
+    end
+        
+    return nothing
+end

src/main/test.jl

+using CUDA, Logging, Random
+
+CUDA.allowscalar(true)
+import Pkg
+Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU")
+#Pkg.activate("/home/alberto/code/julia/LatticeGPU")
+using LatticeGPU
+
+println(CUDA.device())
+
+
+GRP  = SU3
+ALG  = SU3alg
+PREC = Float64
+lp = SpaceParm{4}((64,64,64,64), (4,4,4,4))
+gp = GaugeParm{PREC}(6.0, 5.0/3.0, (0.0,0.0), 3)
+
+println("Space  Parameters: ", lp)
+println("Gauge  Parameters: ", gp)
+
+println("Seeding CURAND...")
+Random.seed!(CURAND.default_rng(), 1234)
+Random.seed!(1234)
+
+println("Precision: ", PREC)
+
+println("Allocating gauge field")
+U = vector_field(GRP{PREC}, lp)
+fill!(U, one(GRP{PREC}))
+
+println("Take to take the configuration to memory: ")
+@time Ucpu = Array(U)
+
+
+println("Allocating YM workspace")
+ymws = YMworkspace(GRP, PREC, lp)
+
+#CUDA.@sync begin
+#    @device_code_warntype CUDA.@cuda threads=lp.bsz blocks=lp.rsz LatticeGPU.krnl_plaq!(ymws.cm, U, lp)
+#end
+    
+@time S = gauge_action(U, lp, gp, ymws)
+@time S = gauge_action(U, lp, gp, ymws)
+println("Initial Action: ", S)
+
+
+println(" - Randomize momenta: ")
+@time randomize!(ymws.frc1, lp, ymws)
+
+eps = 0.1
+ns  = 10
+ntot = 3
+pl = Vector{Float64}()
+for i in 1:3
+    @time dh, acc = HMC!(U,eps,ns,lp, gp, ymws, noacc=true)
+    println("# HMC: ", acc, " ", dh)
+    push!(pl, plaquette(U,lp, gp, ymws))
+    println("# Plaquette: ", pl[end], "\n")
+end
+for i in 1:ntot
+    CUDA.@profile dh, acc = HMC!(U,eps,ns,lp, gp, ymws)
+    println("# HMC: ", acc, " ", dh)
+    push!(pl, plaquette(U,lp, gp, ymws))
+    println("# Plaquette: ", pl[end], "\n")
+end
+
+@time wfl_rk3(U, 1, 0.01, lp, ymws)
+
+println("Action: ", gauge_action(U, lp, gp, ymws))
+println("Time for 100 steps of RK3 flow integrator: ")
+@time wfl_rk3(U, 100, 0.01, lp, ymws)
+println("Action: ", gauge_action(U, lp, gp, ymws))
+
+
+println("END")
+

src/main/times.jl

+using CUDA, Logging, Random, TimerOutputs
+
+CUDA.allowscalar(false)
+import Pkg
+Pkg.activate("../..")
+#Pkg.activate("/home/alberto/code/julia/LatticeGPU")
+using LatticeGPU
+
+# Set lattice/block size
+ntwist = (0,0,0,0,0,1)
+lp = SpaceParm{4}((32,32,32,32), (4,4,4,4), BC_PERIODIC, ntwist)
+println("Space  Parameters: ", lp)
+
+# Seed RNG
+println("Seeding CURAND...")
+Random.seed!(CURAND.default_rng(), 1234)
+Random.seed!(1234)
+
+# Set group and precision
+GRP  = SU3
+ALG  = SU3alg
+PREC = Float64
+println("Precision: ", PREC)
+
+# MD integrator
+int = omf4(PREC, 1.0/80.0, 80)
+println(int)
+
+println("Allocating YM workspace")
+ymws = YMworkspace(GRP, PREC, lp)
+
+# Main program
+println("Allocating gauge field")
+U = vector_field(GRP{PREC}, lp)
+fill!(U, one(GRP{PREC}))
+
+println("Time to take the configuration to memory: ")
+@time Ucpu = Array(U)
+
+
+# Set gauge parameters
+# FIRST SET: Wilson action/flow
+println("\n## WILSON ACTION/FLOW TIMES")
+#gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 1.0)
+gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 1.6666666666666666)
+println("Gauge  Parameters: ", gp)
+
+flwint = wfl_rk3(PREC, 0.005, 1.0E-6)
+println(flwint)
+
+println("Initial Action: ")
+@time S = gauge_action(U, lp, gp, ymws)
+
+
+HMC!(U,int.eps,1,lp, gp, ymws, noacc=true)
+
+pl = Vector{Float64}()
+for i in 1:4
+    @time dh, acc = HMC!(U,int,lp, gp, ymws)
+    println("# HMC: ", acc, " ", dh)
+    push!(pl, plaquette(U,lp, gp, ymws))
+    println("# Plaquette: ", pl[end], "\n")
+
+#    @time x, y = sfcoupling(U,lp,gp,ymws)
+#    println("SF coupling: ", x, " ", y)
+end
+
+Ucp = Array(U)
+println("Action: ", gauge_action(U, lp, gp, ymws))
+flw(U, flwint, 1, gp, lp, ymws)
+
+println("Time for 100 steps of RK3 flow integrator: ")
+@time flw(U, flwint, 100, gp, lp, ymws)
+eoft = Eoft_plaq(U, gp, lp, ymws)
+eoft = Eoft_clover(U, gp, lp, ymws)
+qtop = Qtop(U, gp, lp, ymws)
+
+@time eoft = Eoft_plaq(U, gp, lp, ymws)
+println("Plaq: ", eoft)
+@time eoft = Eoft_clover(U, gp, lp, ymws)
+println("Clov: ", eoft)
+@time qtop = Qtop(U, gp, lp, ymws)
+println("Qtop: ", qtop)
+println("## END Wilson action/flow measurements")
+
+# Set gauge parameters
+# SECOND SET: Improved action/flow
+println("\n## IMPROVED ACTION/FLOW TIMES")
+gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 5/6)
+println("Gauge  Parameters: ", gp)
+
+flwint = zfl_rk3(PREC, 0.01, 1.0E-6)
+println(flwint)
+
+println("Initial Action: ")
+@time S = gauge_action(U, lp, gp, ymws)
+
+
+HMC!(U,int.eps,1,lp, gp, ymws, noacc=true)
+pl = Vector{Float64}()
+for i in 1:4
+    @time dh, acc = HMC!(U,int,lp, gp, ymws, noacc=true)
+    println("# HMC: ", acc, " ", dh)
+    push!(pl, plaquette(U,lp, gp, ymws))
+    println("# Plaquette: ", pl[end], "\n")
+end
+
+flw(U, flwint, 1, gp, lp, ymws)
+
+println("Action: ", gauge_action(U, lp, gp, ymws))
+println("Time for 100 steps of RK3 flow integrator: ")
+@time flw(U, flwint, 100, gp, lp, ymws)
+println("Action: ", gauge_action(U, lp, gp, ymws))
+println("## END improved action/flow measurements")
+
+print_timer(linechars = :ascii)
+println("END")
+

test/flow/test_adapt.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    test_adapt.jl
+### created: Mon Jun  6 12:01:36 2022
+###                               
+
+using LatticeGPU, Test, CUDA
+
+T = Float64
+lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), BC_PERIODIC, (0,0,0,0,0,0))
+gp = GaugeParm{T}(SU3{T}, 6.1, 1.0)
+ymws = YMworkspace(SU3, T, lp)
+
+randomize!(ymws.mom, lp, ymws)
+U = exp.(ymws.mom)
+
+Ucp = deepcopy(U)
+# First Integrate very precisely up to t=2 (Wilson)
+println(" # Very precise integration ")
+wflw = wfl_rk3(Float64, 0.0004, 1.0E-7)
+flw(U, wflw, 5000, gp, lp, ymws)
+pl_exact = Eoft_plaq(U, gp, lp, ymws) 
+cl_exact = Eoft_clover(U, gp, lp, ymws)
+println("  - Plaq:   ", pl_exact)
+println("  - Clover: ", cl_exact)
+Ufin = deepcopy(U)
+
+
+# Now use Adaptive step size integrator:
+for tol in (1.0E-4, 1.0E-5, 1.0E-6, 1.0E-7, 1.0E-8)
+    local wflw = wfl_rk3(Float64, 0.0001, tol)
+    U .= Ucp
+    ns, eps = flw_adapt(U, wflw, 2.0, gp, lp, ymws)
+    pl = Eoft_plaq(U, gp, lp, ymws) 
+    cl = Eoft_clover(U, gp, lp, ymws)
+
+    println(" # Adaptive integrator (tol=$tol): ", ns, " steps")
+    U .= U ./ Ufin
+    maxd = CUDA.mapreduce(dev_one, max, U, init=0.0)
+    println("  - Plaq:   ", pl," [diff: ", abs(pl-pl_exact), "; ",
+            maxd, "]")
+    println("  - Clover: ", cl, " [diff: ", abs(cl-cl_exact), "; ",
+            maxd, "]")
+end

test/io/test_io.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    test_io.jl
+### created: Thu Jun  2 18:20:42 2022
+###                               
+
+using LatticeGPU, Test
+
+T = Float64
+lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), BC_PERIODIC, (0,0,0,0,0,0))
+gp = GaugeParm{T}(SU3{T}, 6.1, 1.0)
+ymws = YMworkspace(SU3, T, lp)
+
+randomize!(ymws.mom, lp, ymws)
+U = exp.(ymws.mom)
+
+pl1 = plaquette(U, lp, gp, ymws)
+cl1 = Eoft_clover(U, gp, lp, ymws)
+fn = "foo.bdio"
+rm(fn, force=true)
+save_cnfg(fn, U, lp, gp; run="Dumnmy_run")
+
+Ucp = read_cnfg(fn)
+pl2 = plaquette(Ucp, lp, gp, ymws)
+cl2 = Eoft_clover(Ucp, gp, lp, ymws)
+rm(fn, force=true)
+
+@testset "Testing configuration i/o" begin
+    @test isapprox(pl1,pl2)
+    @test isapprox(cl1,cl2)
+end
+

test/runtests.jl

+
+#include("SAD/test_sad.jl")
+include("flow/test_adapt.jl")