Added clover type stable function

src/Dirac/Dirac.jl

@@ -234,6 +234,18 @@ function g5Dw!(so, U, si, dpar, lp::SpaceParm{4,6,B,D}) where {B,D}
     return nothing
 end
 
+function clover(g1::SU3{T},g2::SU3{T},g3::SU3{T},g4::SU3{T}) where {T}
+    u11 = g1.u11 + g2.u11 + g3.u11 + g4.u11
+    u12 = g1.u12 + g2.u12 + g3.u12 + g4.u12
+    u13 = g1.u13 + g2.u13 + g3.u13 + g4.u13
+    u21 = g1.u21 + g2.u21 + g3.u21 + g4.u21
+    u22 = g1.u22 + g2.u22 + g3.u22 + g4.u22
+    u23 = g1.u23 + g2.u23 + g3.u23 + g4.u23
+    u31 = conj(g1.u12*g1.u23 - g1.u13*g1.u22) + conj(g2.u12*g2.u23 - g2.u13*g2.u22) + conj(g3.u12*g3.u23 - g3.u13*g3.u22) + conj(g4.u12*g4.u23 - g4.u13*g4.u22)
+    u32 = conj(g1.u13*g1.u21 - g1.u11*g1.u23) + conj(g2.u13*g2.u21 - g2.u11*g2.u23) + conj(g3.u13*g3.u21 - g3.u11*g3.u23) + conj(g4.u13*g4.u21 - g4.u11*g4.u23)
+    u33 = conj(g1.u11*g1.u22 - g1.u12*g1.u21) + conj(g2.u11*g2.u22 - g2.u12*g2.u21) + conj(g3.u11*g3.u22 - g3.u12*g3.u21) + conj(g4.u11*g4.u22 - g4.u12*g4.u21)
+    return M3x3{T}(u11-conj(u11),u12-conj(u21),u13-conj(u31),u21-conj(u12),u22-conj(u22),u23-conj(u32),u31-conj(u13),u32-conj(u23),u33-conj(u33))
+end
 
 function Csw!(dws, U, gp, lp::SpaceParm{4,6,B,D}) where {B,D}
 
@@ -377,17 +389,17 @@ function krnl_csw!(csw::AbstractArray{T}, U, Ubnd, ipl, ztw, lp::SpaceParm{4,M,B
             gt2 = U[bud,id2,rud]
         end
 
-        M1 = convert(M3x3{T}, U[b,id1,r]*gt1/(U[b,id2,r]*U[bu2,id1,ru2]))
-        M2 = convert(M3x3{T}, (U[bd2,id2,rd2]\(U[bd2,id1,rd2]*gt2))/U[b,id1,r])
-        M3 = convert(M3x3{T}, (U[bdd,id2,rdd]*U[bd1,id1,rd1])\(U[bdd,id1,rdd]*U[bd2,id2,rd2]))
-        M4 = convert(M3x3{T}, (U[b,id2,r]/(U[bd1,id2,rd1]*U[bdu,id1,rdu]))*U[bd1,id1,rd1])
+        M1 = U[b,id1,r]*gt1/(U[b,id2,r]*U[bu2,id1,ru2])
+        M2 = (U[bd2,id2,rd2]\(U[bd2,id1,rd2]*gt2))/U[b,id1,r]
+        M3 = (U[bdd,id2,rdd]*U[bd1,id1,rd1])\(U[bdd,id1,rdd]*U[bd2,id2,rd2])
+        M4 =  (U[b,id2,r]/(U[bd1,id2,rd1]*U[bdu,id1,rdu]))*U[bd1,id1,rd1]
 
 
         if !(SFBC && (it == 1))
             if TWP
-                csw[b,ipl,r]  = ztw*(M1 + M2 + M3 + M4 - dag(M1 + M2 + M3 + M4))
+                csw[b,ipl,r]  = ztw*clover(M1,M2,M3,M4)
             else
-                csw[b,ipl,r]  = M1 + M2 + M3 + M4 - dag(M1 + M2 + M3 + M4)
+                csw[b,ipl,r]  = clover(M1,M2,M3,M4)
             end
         end