First test for improved fermions on Dw!

src/Dirac/Dirac.jl

@@ -45,7 +45,7 @@ struct DiracWorkspace{T}
 end
 export DiracWorkspace, DiracParam
 
-function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}) where {B,D}
+function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}, ymws = "ymws", gp = "gp") where {B,D}
 
     if B == BC_SF_AFWB || B == BC_SF_ORBI
         CUDA.@sync begin
@@ -67,6 +67,22 @@ function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}) where {B,D}
             end
             end
     end 
+
+    if abs(dpar.csw) > 1.0E-10
+        if ymws == "ymws" || gp == "gp"
+            error("YM workspace and Gauge parameters must be included for the improved dirac operator.")
+        else
+        @timeit "Dw_improvement" begin
+            for i in Int(lp.npls/2)
+                z1 = exp(2im * pi * lp.ntw[i]/gp.ng)
+                z2 = exp(2im * pi * lp.ntw[lp.npls-i]/gp.ng)
+                    CUDA.@sync begin
+                        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, i, lp.npls-i, z1, z2, lp)
+                        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw_impr!(so, ymws.frc1, ymws.frc2, i, lp.npls-i, si, csw, lp)
+                    end
+                end
+        end
+    end
     
     return nothing
 end
@@ -212,6 +228,19 @@ function krnl_g5Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
 end
 
 
+function krnl_Dw_impr!(so, F1, F2, id1, id2, si, csw, lp::SpaceParm{4,6,B,D}) where {B,D}
+
+    b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
+
+    @inbounds begin 
+        # asignation between /mu/nu index and sigma index is NOT correct
+        so[b,r] += 0.25*csw*im*( F1*dmul(Gamma{id1},si[b,r]) + F2*dmul(Gamma{id2},si[b,r]))
+        
+    end
+        
+    return nothing
+end
+
 
 ###############################   HMC for fermions   ###################################
 

src/Groups/FundamentalSU3.jl

@@ -60,6 +60,16 @@ Base.:\(g::SU3{T},b::SU3fund{T})     where T <: AbstractFloat = SU3fund{T}(conj(
                                                                            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)
 
+"""
+    *(a::SU3alg{T},b::SU3fund{T})
+
+Returns a*b
+"""
+
+Base.:*(a::SU3alg{T},b::SU3fund{T})     where T <: AbstractFloat = SU3fund{T}(complex(0.0, a.t7/2.0 + a.t8/3.46410161513775458)*b.t1 + complex(a.t4,a.t1)*b.t2/2.0 + complex(a.t5,a.t2)*b.t3/2.0,
+                                                                              complex(-a.t4,a.t1)*b.t1/2.0 + complex(0.0,-a.t7/2.0 + a.t8/3.46410161513775458)*b.t2 + complex(a.t6,a.t3)*b.t3/2.0,
+                                                                              complex(-a.t5,a.t2)*b.t1/2.0 + complex(-a.t6,a.t3)/2.0*b.t2 + complex(0.0,-2.0*a.t8/3.46410161513775458)*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)

src/Spinors/Spinors.jl

@@ -89,6 +89,13 @@ Returns ga
 """
 Base.:*(g::S,b::Spinor{NS,G}) where {S <: Group,NS,G} = Spinor{NS,G}(ntuple(i->g*b.s[i], NS))
 
+"""
+    *(a::SU3alg{T},b::Spinor)
+
+Returns ab
+"""
+Base.:*(a::S,b::Spinor{NS,G}) where {S <: Algebra,NS,G} = Spinor{NS,G}(ntuple(i->a*b.s[i], NS))
+
 """
     \\(g::SU3{T},b::Spinor{NS,G})