improvement term added for g5Dw

src/Dirac/Dirac.jl

@@ -88,7 +88,7 @@ function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}, ymws = "ymws",
     return nothing
 end
 
-function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}) where {B,D}
+function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}, ymws = "ymws", gp = "gp") where {B,D}
 
     
 
@@ -129,7 +129,7 @@ function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}) where
 end
 
 
-function g5Dw!(so, U, si, dpar, lp::SpaceParm{4,6,B,D}) where {B,D}
+function g5Dw!(so, U, si, dpar, 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
@@ -152,6 +152,24 @@ function g5Dw!(so, U, si, dpar, lp::SpaceParm{4,6,B,D}) where {B,D}
             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 1: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_g5Dw_impr!(so, ymws.frc1, ymws.frc2, Gamma{i+9}, Gamma{lp.npls-i+8}, si, dpar.csw, lp)
+                    end
+            end
+        end
+        end
+    end
+    
     return nothing
 end
 
@@ -234,7 +252,7 @@ function krnl_Dw_impr!(so, F1, F2, Gid1, Gid2, si, csw, lp::SpaceParm{4,6,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 : sigma 31 -> sigma13 in Spinors end
+        
         so[b,r] += 0.25*csw*im*( F1[b,r]*dmul(Gid1,si[b,r]) - F2[b,r]*dmul(Gid2,si[b,r]))
         
     end
@@ -242,6 +260,19 @@ function krnl_Dw_impr!(so, F1, F2, Gid1, Gid2, si, csw, lp::SpaceParm{4,6,B,D})
     return nothing
 end
 
+function krnl_g5Dw_impr!(so, F1, F2, Gid1, Gid2, si, csw, lp::SpaceParm{4,6,B,D}) where {B,D}
+
+    b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
+
+    @inbounds begin 
+        
+        so[b,r] += 0.25*csw*im*dmul(Gamma{5},( F1[b,r]*dmul(Gid1,si[b,r]) - F2[b,r]*dmul(Gid2,si[b,r])))
+
+    end
+        
+    return nothing
+end
+
 
 ###############################   HMC for fermions   ###################################