Change in Csw: now a matrix in dws

parent e2a0507f
......@@ -50,6 +50,7 @@ struct DiracWorkspace{T}
sp
sAp
st
csw
function DiracWorkspace(::Type{G}, ::Type{T}, lp::SpaceParm{4,6,B,D}) where {G,T <: AbstractFloat, B,D}
......@@ -57,12 +58,14 @@ struct DiracWorkspace{T}
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)
csw = tensor_field(M3x3{T},lp)
return new{T}(sr,sp,sAp,st,csw)
end
end
export DiracWorkspace, DiracParam
function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}, ymws::YMworkspace) where {B,D}
function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}) where {B,D}
if B == BC_SF_AFWB || B == BC_SF_ORBI
......@@ -95,7 +98,7 @@ function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}, ymws::YMworksp
if abs(dpar.csw) > 1.0E-10
@timeit "Dw_improvement" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw_impr!(so, ymws.csw, dpar.csw, si, lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw_impr!(so, dws.csw, dpar.csw, si, lp)
end
end
end
......@@ -103,7 +106,7 @@ function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}, ymws::YMworksp
return nothing
end
function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}, ymws::YMworkspace) where {B,D}
function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}) where {B,D}
if B == BC_SF_AFWB || B == BC_SF_ORBI
@timeit "DwdagDw" begin
......@@ -129,7 +132,7 @@ function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}, ymws:
if abs(dpar.csw) > 1.0E-10
@timeit "Dw_improvement" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw_impr!(st, ymws.csw, dpar.csw, si, lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw_impr!(st, dws.csw, dpar.csw, si, lp)
end
end
end
......@@ -151,7 +154,7 @@ function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}, ymws:
if abs(dpar.csw) > 1.0E-10
@timeit "Dw_improvement" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw_impr!(so, ymws.csw, dpar.csw, st, lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw_impr!(so, dws.csw, dpar.csw, st, lp)
end
end
end
......@@ -168,7 +171,7 @@ function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}, ymws:
if abs(dpar.csw) > 1.0E-10
@timeit "Dw_improvement" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw_impr!(st, ymws.csw, dpar.csw, si, lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw_impr!(st, dws.csw, dpar.csw, si, lp)
end
end
end
......@@ -180,7 +183,7 @@ function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}, ymws:
if abs(dpar.csw) > 1.0E-10
@timeit "Dw_improvement" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw_impr!(so, ymws.csw, dpar.csw, st, lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw_impr!(so, dws.csw, dpar.csw, st, lp)
end
end
end
......@@ -192,7 +195,7 @@ function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}, ymws:
end
function g5Dw!(so, U, si, dpar, lp::SpaceParm{4,6,B,D}, ymws::YMworkspace) where {B,D}
function g5Dw!(so, U, si, dpar, lp::SpaceParm{4,6,B,D}) where {B,D}
if B == BC_SF_AFWB || B == BC_SF_ORBI
@timeit "SF fix" begin
......@@ -222,7 +225,7 @@ function g5Dw!(so, U, si, dpar, lp::SpaceParm{4,6,B,D}, ymws::YMworkspace) where
if abs(dpar.csw) > 1.0E-10
@timeit "Dw_improvement" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw_impr!(so, ymws.csw, dpar.csw, si, lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw_impr!(so, dws.csw, dpar.csw, si, lp)
end
end
end
......@@ -232,14 +235,14 @@ function g5Dw!(so, U, si, dpar, lp::SpaceParm{4,6,B,D}, ymws::YMworkspace) where
end
function Csw!(ymws, U, gp, lp::SpaceParm{4,6,B,D}) where {B,D}
function Csw!(dws, U, gp, lp::SpaceParm{4,6,B,D}) where {B,D}
@timeit "Csw computation" begin
for i in 1:Int(lp.npls)
z = exp(2im * pi * lp.ntw[i]/gp.ng)
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_csw!(ymws.csw, U, gp.Ubnd, i, z, lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_csw!(dws.csw, U, gp.Ubnd, i, z, lp)
end
end
end
......@@ -374,14 +377,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])
if !(SFBC && (it == 1))
if TWP
csw[b,ipl,r] = projalg(ztw,U[b,id1,r]*gt1/(U[b,id2,r]*U[bu2,id1,ru2]) + (U[bd2,id2,rd2]\(U[bd2,id1,rd2]*gt2))/U[b,id1,r] +
(U[bdd,id2,rdd]*U[bd1,id1,rd1])\(U[bdd,id1,rdd]*U[bd2,id2,rd2]) + (U[b,id2,r]/(U[bd1,id2,rd1]*U[bdu,id1,rdu]))*U[bd1,id1,rd1])
csw[b,ipl,r] = ztw*(M1 + M2 + M3 + M4 - dag(M1 + M2 + M3 + M4))
else
csw[b,ipl,r] = projalg(U[b,id1,r]*gt1/(U[b,id2,r]*U[bu2,id1,ru2]) + (U[bd2,id2,rd2]\(U[bd2,id1,rd2]*gt2))/U[b,id1,r] +
(U[bdd,id2,rdd]*U[bd1,id1,rd1])\(U[bdd,id1,rdd]*U[bd2,id2,rd2]) + (U[b,id2,r]/(U[bd1,id2,rd1]*U[bdu,id1,rdu]))*U[bd1,id1,rd1])
csw[b,ipl,r] = M1 + M2 + M3 + M4 - dag(M1 + M2 + M3 + M4)
end
end
......
......@@ -86,3 +86,6 @@ Base.:*(a::SU3fund{T},b::Number) where T <: AbstractFloat = SU3fund{T}(b*a.t1,b*
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)
Base.:*(a::M3x3{T}),b::SU3fund{T} where T <: AbstractFloat = SU3fund{T}(a.u11*b.t1 + a.u12*b.t2 + a.u13*b.t3,
a.u21*b.t1 + a.u22*b.t2 + a.u23*b.t3,
a.u31*b.t1 + a.u32*b.t2 + a.u33*b.t3)
......@@ -90,4 +90,3 @@ function isgroup(a::SU3{T}) where T <: AbstractFloat
end
end
Base.:+(a::SU3{T},b::SU3{T}) where T <: AbstractFloat = SU3{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)
......@@ -69,7 +69,6 @@ struct YMworkspace{T}
PRC
frc1
frc2
csw
mom
U1
cm # complex of volume
......@@ -84,7 +83,6 @@ struct YMworkspace{T}
f2 = vector_field(SU2alg{T}, lp)
mm = vector_field(SU2alg{T}, lp)
u1 = vector_field(SU2{T}, lp)
csw = tensor_field(SU2alg{T},lp)
end
if (G == SU3)
......@@ -94,13 +92,12 @@ struct YMworkspace{T}
f2 = vector_field(SU3alg{T}, lp)
mm = vector_field(SU3alg{T}, lp)
u1 = vector_field(SU3{T}, lp)
csw = tensor_field(SU3alg{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,csw, mm, u1, cs, rs)
return new{T}(GRP,ALG,T,f1, f2, mm, u1, cs, rs)
end
end
export YMworkspace
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment