diff --git a/src/YM/YM.jl b/src/YM/YM.jl index d5438b5..b5672e5 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -20,6 +20,7 @@ import Base.show struct GaugeParm{T} beta::T + c0::T cG::NTuple{2,T} ng::Int64 end @@ -27,6 +28,7 @@ export GaugeParm function Base.show(io::IO, gp::GaugeParm) println(io, "beta: ", gp.beta) + println(io, "c0: ", gp.c0) println(io, "Ngauge: ", gp.ng) return nothing @@ -38,23 +40,17 @@ struct YMworkspace{T} PRC frc1 frc2 - fpln mom U1 cm # complex of volume rm # float of volume - function YMworkspace(::Type{G}, ::Type{T}, lp::SpaceParm; save_mem = false) where {G <: Group, T <: AbstractFloat} + function YMworkspace(::Type{G}, ::Type{T}, lp::SpaceParm) where {G <: Group, T <: AbstractFloat} if (G == SU2) GRP = SU2 ALG = SU2alg f1 = field(SU2alg{T}, lp) f2 = field(SU2alg{T}, lp) - if save_mem - fp = nothing - else - fp = field_pln(SU2alg{T}, lp) - end mm = field(SU2alg{T}, lp) u1 = field(SU2{T}, lp) end @@ -64,18 +60,13 @@ struct YMworkspace{T} ALG = SU3alg f1 = field(SU3alg{T}, lp) f2 = field(SU3alg{T}, lp) - if save_mem - fp = nothing - else - fp = field_pln(SU2alg{T}, lp) - end mm = field(SU3alg{T}, lp) u1 = field(SU3{T}, lp) end cs = CuArray{Complex{T}, 2}(undef, lp.bsz,lp.rsz) rs = CuArray{T, 2}(undef, lp.bsz,lp.rsz) - return new{T}(GRP,ALG,T,f1, f2, fp, mm, u1, cs, rs) + return new{T}(GRP,ALG,T,f1, f2, mm, u1, cs, rs) end end export YMworkspace diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl index 26fa4ad..75105a2 100644 --- a/src/YM/YMact.jl +++ b/src/YM/YMact.jl @@ -9,6 +9,88 @@ ### created: Mon Jul 12 18:31:19 2021 ### +function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,D}) where {T,N,M,D} + + b, r = CUDA.threadIdx().x, CUDA.blockIdx().x + + Ush = @cuStaticSharedMem(T, (D,2)) + + plx[b,r] = zero(plx[b,r]) + for id1 in 1:N-1 + bu1, ru1 = up((b, r), id1, lp) + Ush[b,1] = U[b,id1,r] + + for id2 = id1+1:N + bu2, ru2 = up((b, r), id2, lp) + Ush[b,2] = U[b,id2,r] + sync_threads() + + # H2 staple + (b1, r1) = up((b,r), id1, lp) + if r1 == r + ga = Ush[b1,1] + else + ga = U[b1,id1,r1] + end + + (b2, r2) = up((b1,r1), id1, lp) + if r2 == r + gb = Ush[b2,2] + else + gb = U[b2,id2,r2] + end + + (b2, r2) = up((b1,r1), id2, lp) + if r2 == r + gc = Ush[b2,1] + else + gc = U[b2,id1,r2] + end + h2 = (ga*gb)/gc + + # H3 staple + (b1, r1) = up((b,r), id2, lp) + if r1 == r + ga = Ush[b1,2] + else + ga = U[b1,id2,r1] + end + + (b2, r2) = up((b1,r1), id2, lp) + if r2 == r + gb = Ush[b2,1] + else + gb = U[b2,id1,r2] + end + + (b2, r2) = up((b1,r1), id1, lp) + if r2 == r + gc = Ush[b2,2] + else + gc = U[b2,id2,r2] + end + h3 = (ga*gb)/gc + + if ru2 == r + gb = Ush[bu2,1] + else + gb = U[bu2,id1,ru2] + end + if ru1 == r + ga = Ush[bu1,2] + else + ga = U[bu1,id2,ru1] + end + + g2 = Ush[b,2]\Ush[b,1] + + plx[b,r] += c0*tr(g2*ga/gb) + c1*( tr(g2*h2/gb) + tr(g2*ga/h3)) + end + end + + return nothing +end + function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,D}) where {T,N,M,D} b, r = CUDA.threadIdx().x, CUDA.blockIdx().x @@ -121,14 +203,14 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp:: if r1 == r ga = Ush[b1,1] else - ga = U[b1,id2,r1] + ga = U[b1,id1,r1] end (b2, r2) = up((b1,r1), id1, lp) if r2 == r - gc = Ush[b2,2] + gb = Ush[b2,2] else - gc = U[b2,id2,r2] + gb = U[b2,id2,r2] end (b2, r2) = up((b1,r1), id2, lp) @@ -149,9 +231,9 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp:: (b2, r2) = up((b1,r1), id2, lp) if r2 == r - gc = Ush[b2,1] + gb = Ush[b2,1] else - gc = U[b2,id1,r2] + gb = U[b2,id1,r2] end (b2, r2) = up((b1,r1), id1, lp) @@ -160,7 +242,7 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp:: else gc = U[b2,id2,r2] end - h3 = (ga\gb)/gc + h3 = (ga*gb)/gc # H4 staple (b1, r1) = dw((b,r), id1, lp) @@ -207,67 +289,7 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp:: frc2[bu2,id1,ru2] += c0*projalg(g2*g1) + c1*( projalg((Ush[b,2]\h1)*g1) + projalg(g2*h2/gb) + - projalg(h4*Ush[b,1]*g1) ) - end - - return nothing -end - -function krnl_force_wilson_nw!(fpl, U::AbstractArray{T}, lp::SpaceParm{N,M,D}) where {T,N,M,D} - - b, r = CUDA.threadIdx().x, CUDA.blockIdx().x - - Ush = @cuStaticSharedMem(T, (D,N)) - - for id in 1:N - @inbounds Ush[b,id] = U[b,id,r] - end - sync_threads() - - @inbounds for ipl in 1:lp.npls - id1, id2 = lp.plidx[ipl] - - bu1, ru1 = up((b, r), id1, lp) - bu2, ru2 = up((b, r), id2, lp) - - if ru2 == r - gt2 = Ush[bu2,id1] - else - gt2 = U[bu2,id1,ru2] - end - if ru1 == r - gt1 = Ush[bu1,id2] - else - gt1 = U[bu1,id2,ru1] - end - - g1 = gt1/gt2 - g2 = Ush[b,id2]\Ush[b,id1] - - fpl[b ,ipl, r ,1] = projalg(Ush[b,id1]*g1/Ush[b,id2]) - fpl[bu1,ipl, ru1,2] = projalg(g1*g2) - fpl[bu2,ipl, ru2,3] = projalg(g2*g1) - - end - - return nothing -end - -function krnl_add_force_plns!(frc::AbstractArray{T}, fpl, lp::SpaceParm{N,M,D}) where {T,N,M,D} - - b, r = CUDA.threadIdx().x, CUDA.blockIdx().x - - Fsh = @cuStaticSharedMem(T, (D,M)) - @inbounds for j in 1:M - Fsh[b,j] = fpl[b,j,r,1] - end - sync_threads() - - @inbounds for ipl in 1:lp.npls - id1, id2 = lp.plidx[ipl] - - frc[b,id1,r] += -Fsh[b,ipl] + fpl[b,ipl,r,3] - frc[b,id2,r] += Fsh[b,ipl] - fpl[b,ipl,r,2] + projalg(h4\Ush[b,1]*g1) ) end return nothing @@ -279,25 +301,34 @@ end Computes the force deriving from the Wilson plaquette action, without the prefactor 1/g0^2, and assign it to the workspace force `ymws.frc1` """ -function force_wilson(ymws::YMworkspace, U, lp::SpaceParm) +function force_gauge(ymws::YMworkspace, U, c0, lp::SpaceParm) - if ymws.fpln == nothing + if abs(c0-1) < 1.0E-10 force_wilson_pln!(ymws.frc1, ymws.frc2, U, lp::SpaceParm) else - force_wilson_nw!(ymws.frc1, ymws.fpln, U, lp::SpaceParm) + force_wilson_pln!(ymws.frc1, ymws.frc2, U, lp::SpaceParm, c0) end return nothing end +force_wilson(ymws::YMworkspace, U, lp::SpaceParm) = force_gauge(ymws, U, 1, lp) -function force_wilson_pln!(frc1, ftmp, U, lp::SpaceParm) + +function force_wilson_pln!(frc1, ftmp, U, lp::SpaceParm, c0=1) fill!(frc1, zero(eltype(frc1))) fill!(ftmp, zero(eltype(ftmp))) - for i in 1:lp.npls - CUDA.@sync begin -# CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_wilson_pln!(frc1,ftmp,U,i,lp) - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_impr_pln!(frc1,ftmp,U,1.0,0.0,i,lp) + if c0 == 1 + for i in 1:lp.npls + CUDA.@sync begin + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_wilson_pln!(frc1,ftmp,U,i,lp) + end + end + else + for i in 1:lp.npls + CUDA.@sync begin + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_impr_pln!(frc1,ftmp,U,c0,(1-c0)/8,i,lp) + end end end frc1 .= frc1 .+ ftmp @@ -305,15 +336,3 @@ function force_wilson_pln!(frc1, ftmp, U, lp::SpaceParm) return nothing end -function force_wilson_nw!(frc1, fpln, U, lp::SpaceParm) - - CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_wilson_nw!(fpln,U,lp) - end - fill!(frc1, zero(eltype(frc1))) - CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_force_plns!(frc1, fpln,lp) - end - - return nothing -end diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl index 8b02100..0a3f4ca 100644 --- a/src/YM/YMhmc.jl +++ b/src/YM/YMhmc.jl @@ -16,11 +16,17 @@ Returns the value of the Wilson plaquette action for the configuration U. """ function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: AbstractFloat - - CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, lp) + + if abs(gp.c0-1) < 1.0E-10 + CUDA.@sync begin + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, lp) + end + else + CUDA.@sync begin + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_impr!(ymws.cm, U, gp.c0, (1-gp.c0)/8, lp) + end end - S = gp.beta*( prod(lp.iL)*lp.npls - + S = gp.beta*( prod(lp.iL)*lp.npls*(gp.c0 + (1-gp.c0)/8) - CUDA.mapreduce(real, +, ymws.cm)/gp.ng ) return S @@ -80,28 +86,28 @@ function OMF4!(mom, U, eps, ns, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspa r6::T = 1.0-2.0*(r2+r4) ee = eps*gp.beta/gp.ng - force_wilson(ymws, U, lp) + force_gauge(ymws, U, gp.c0, lp) for i in 1:ns mom .= mom .+ (r1*ee) .* ymws.frc1 U .= expm.(U, mom, eps*r2) - force_wilson(ymws, U, lp) + force_gauge(ymws, U, gp.c0, lp) mom .= mom .+ (r3*ee) .* ymws.frc1 U .= expm.(U, mom, eps*r4) - force_wilson(ymws, U, lp) + force_gauge(ymws, U, gp.c0, lp) mom .= mom .+ (r5*ee) .* ymws.frc1 U .= expm.(U, mom, eps*r6) - force_wilson(ymws, U, lp) + force_gauge(ymws, U, gp.c0, lp) mom .= mom .+ (r5*ee) .* ymws.frc1 U .= expm.(U, mom, eps*r4) - force_wilson(ymws, U, lp) + force_gauge(ymws, U, gp.c0, lp) mom .= mom .+ (r3*ee) .* ymws.frc1 U .= expm.(U, mom, eps*r2) - force_wilson(ymws, U, lp) + force_gauge(ymws, U, gp.c0, lp) mom .= mom .+ (r1*ee) .* ymws.frc1 end diff --git a/src/main/test.jl b/src/main/test.jl index 2290676..1268753 100644 --- a/src/main/test.jl +++ b/src/main/test.jl @@ -13,7 +13,7 @@ GRP = SU3 ALG = SU3alg PREC = Float64 lp = SpaceParm{4}((32,32,32,32), (4,4,4,4)) -gp = GaugeParm{PREC}(6.0, (0.0,0.0), 3) +gp = GaugeParm{PREC}(6.0, 5.0/3.0, (0.0,0.0), 3) println("Space Parameters: ", lp) println("Gauge Parameters: ", gp) @@ -33,7 +33,7 @@ println("Take to take the configuration to memory: ") println("Allocating YM workspace") -ymws = YMworkspace(GRP, PREC, lp, save_mem = true) +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)