mirror of
https://igit.ific.uv.es/alramos/latticegpu.jl.git
synced 2025-05-15 03:33:42 +02:00
First commit with layout for Latttice GPU code
This commit is contained in:
Alberto Ramos
commit
5bb4f28c8b
9 changed files with 747 additions and 0 deletions
146
src/Groups/GroupSU2.jl
Normal file
146
src/Groups/GroupSU2.jl
Normal file
|
|
@ -0,0 +1,146 @@
|
|||
###
|
||||
### "THE BEER-WARE LICENSE":
|
||||
### Alberto Ramos wrote this file. As long as you retain this
|
||||
### notice you can do whatever you want with this stuff. If we meet some
|
||||
### day, and you think this stuff is worth it, you can buy me a beer in
|
||||
### return. <alberto.ramos@cern.ch>
|
||||
###
|
||||
### file: GroupSU2.jl
|
||||
### created: Sun Jul 11 17:23:12 2021
|
||||
###
|
||||
|
||||
#
|
||||
# SU(2) group elements represented trough Cayley-Dickson
|
||||
# construction
|
||||
# https://en.wikipedia.org/wiki/Cayley%E2%80%93Dickson_construction
|
||||
|
||||
import Base.:*, Base.:+, Base.:-,Base.:/,Base.:\
|
||||
struct SU2 <: Group
|
||||
t1::ComplexF64
|
||||
t2::ComplexF64
|
||||
end
|
||||
SU2() = SU2(1.0, 0.0)
|
||||
inverse(b::SU2) = SU2(conj(b.t1), -b.t2)
|
||||
dag(a::SU2) = inverse(a)
|
||||
norm(a::SU2) = sqrt(abs2(a.t1) + abs2(a.t2))
|
||||
tr(g::SU2) = 2.0*real(a.t1)
|
||||
|
||||
"""
|
||||
function normalize(a::SU2)
|
||||
|
||||
Return a normalized element of `SU(2)`
|
||||
"""
|
||||
function normalize(a::SU2)
|
||||
dr = sqrt(abs2(a.t1) + abs2(a.t2))
|
||||
if (dr == 0.0)
|
||||
return SU2(0.0)
|
||||
end
|
||||
return SU2(a.t1/dr,a.t2/dr)
|
||||
end
|
||||
|
||||
Base.:+(a::SU2,b::SU2) = SU2(a.t1+b.t1,a.t2+b.t2)
|
||||
Base.:-(a::SU2,b::SU2) = SU2(a.t1-b.t1,a.t2-b.t2)
|
||||
Base.:*(a::SU2,b::SU2) = SU2(a.t1*b.t1-a.t2*conj(b.t2),a.t1*b.t2+a.t2*conj(b.t1))
|
||||
Base.:/(a::SU2,b::SU2) = SU2(a.t1*conj(b.t1)+a.t2*conj(b.t2),-a.t1*b.t2+a.t2*b.t1)
|
||||
Base.:\(a::SU2,b::SU2) = SU2(conj(a.t1)*b.t1+a.t2*conj(b.t2),conj(a.t1)*b.t2-a.t2*conj(b.t1))
|
||||
Base.:+(a::SU2) = SU2(a.t1,a.t2)
|
||||
Base.:-(a::SU2) = SU2(-a.t1,-a.t2)
|
||||
|
||||
struct SU2alg <: Algebra
|
||||
t1::Float64
|
||||
t2::Float64
|
||||
t3::Float64
|
||||
end
|
||||
SU2alg(x::Real) = SU2alg(x,0.0,0.0)
|
||||
SU2alg(v::Vector) = SU2alg(v[1],v[2],v[3])
|
||||
projalg(g::SU2) = SU2alg(imag(g.t1), real(g.t2), imag(g.t2))
|
||||
dot(a::SU2alg, b::SU2alg) = a.t1*b.t1 + a.t2*b.t2 + a.t3*b.t3
|
||||
Base.:+(a::SU2alg) = SU2alg(a.t1,a.t2,a.t3)
|
||||
Base.:-(a::SU2alg) = SU2alg(-a.t1,-a.t2,-a.t3)
|
||||
Base.:+(a::SU2alg,b::SU2alg) = SU2alg(a.t1+b.t1,a.t2+b.t2,a.t3+b.t3)
|
||||
Base.:-(a::SU2alg,b::SU2alg) = SU2alg(a.t1-b.t1,a.t2-b.t2,a.t3-b.t3)
|
||||
|
||||
Base.:*(a::SU2alg,b::Real) = SU2alg(a.t1*b,a.t2*b,a.t3*b)
|
||||
Base.:*(b::Real,a::SU2alg) = SU2alg(a.t1*b,a.t2*b,a.t3*b)
|
||||
Base.:/(a::SU2alg,b::Real) = SU2alg(a.t1/b,a.t2/b,a.t3/b)
|
||||
|
||||
|
||||
"""
|
||||
function Base.exp(a::SU2alg, t::Number=1)
|
||||
|
||||
Computes `exp(a)`
|
||||
"""
|
||||
function Base.exp(a::SU2alg)
|
||||
|
||||
rm = sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0
|
||||
if (abs(rm) < 0.05)
|
||||
rms = rm^2/2.0
|
||||
ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0))
|
||||
sa = 0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0))
|
||||
else
|
||||
ca = cos(rm)
|
||||
sa = sin(rm)/(2.0*rm)
|
||||
end
|
||||
|
||||
return SU2(complex(ca,sa*a.t1),complex(sa*a.t2,sa*a.t3))
|
||||
end
|
||||
|
||||
function Base.exp(a::SU2alg, t::Number)
|
||||
|
||||
rm = t*sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0
|
||||
if (abs(rm) < 0.05)
|
||||
rms = rm^2/2.0
|
||||
ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0))
|
||||
sa = t*(0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0)))
|
||||
else
|
||||
ca = cos(rm)
|
||||
sa = t*sin(rm)/(2.0*rm)
|
||||
end
|
||||
|
||||
return SU2(complex(ca,sa*a.t1),complex(sa*a.t2,sa*a.t3))
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
function expm(g::SU2, a::SU2alg)
|
||||
|
||||
Computes `exp(a)*g`
|
||||
|
||||
"""
|
||||
function expm(g::SU2, a::SU2alg)
|
||||
|
||||
rm = sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0
|
||||
if (abs(rm) < 0.05)
|
||||
rms = rm^2/2.0
|
||||
ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0))
|
||||
sa = 0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0))
|
||||
else
|
||||
ca = cos(rm)
|
||||
sa = sin(rm)/(2.0*rm)
|
||||
end
|
||||
|
||||
return SU2(complex(ca,sa*a.t1)*g.t1-complex(sa*a.t2,sa*a.t3)*conj(g.t2),
|
||||
complex(ca,sa*a.t1)*g.t2+complex(sa*a.t2,sa*a.t3)*conj(g.t1))
|
||||
end
|
||||
|
||||
"""
|
||||
function expm(g::SU2, a::SU2alg, t::Float64)
|
||||
|
||||
Computes `exp(t*a)*g`
|
||||
|
||||
"""
|
||||
function expm(g::SU2, a::SU2alg, t::Float64)
|
||||
|
||||
rm = t*sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0
|
||||
if (abs(rm) < 0.05)
|
||||
rms = rm^2/2.0
|
||||
ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0))
|
||||
sa = t*(0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0)))
|
||||
else
|
||||
ca = cos(rm)
|
||||
sa = t*sin(rm)/(2.0*rm)
|
||||
end
|
||||
|
||||
return SU2(complex(ca,sa*a.t1)*g.t1-complex(sa*a.t2,sa*a.t3)*conj(g.t2),
|
||||
complex(ca,sa*a.t1)*g.t2+complex(sa*a.t2,sa*a.t3)*conj(g.t1))
|
||||
end
|
||||
48
src/Groups/GroupSU3.jl
Normal file
48
src/Groups/GroupSU3.jl
Normal file
|
|
@ -0,0 +1,48 @@
|
|||
###
|
||||
### "THE BEER-WARE LICENSE":
|
||||
### Alberto Ramos wrote this file. As long as you retain this
|
||||
### notice you can do whatever you want with this stuff. If we meet some
|
||||
### day, and you think this stuff is worth it, you can buy me a beer in
|
||||
### return. <alberto.ramos@cern.ch>
|
||||
###
|
||||
### file: GroupSU3.jl
|
||||
### created: Sun Jul 11 17:23:02 2021
|
||||
###
|
||||
|
||||
#
|
||||
# Use memory efficient representation: Only store
|
||||
# first two rows
|
||||
#
|
||||
# a.u31 = a.u12*conj(a.u23) - a.u13*conj(a.u22)
|
||||
# a.u32 = a.u13*conj(a.u21) - a.u11*conj(a.u23)
|
||||
# a.u33 = a.u11*conj(a.u22) - a.u12*conj(a.u21)
|
||||
#
|
||||
|
||||
import Base.:*, Base.:+, Base.:-,Base.:/,Base.:\
|
||||
struct SU3 <: Group
|
||||
u11::ComplexF64
|
||||
u12::ComplexF64
|
||||
u13::ComplexF64
|
||||
u21::ComplexF64
|
||||
u22::ComplexF64
|
||||
u23::ComplexF64
|
||||
end
|
||||
SU3() = SU3(1,0,0,0,1,0)
|
||||
inverse(a::SU3) = SU3(conj(a.u11),conj(a.u21),conj(a.u12*conj(a.u23) - a.u13*conj(a.u22)),
|
||||
conj(a.u12),conj(a.u22),conj(a.u13*conj(a.u21) - a.u11*conj(a.u23)))
|
||||
dag(a::SU3) = inverse(a)
|
||||
tr(g::SU3) = a.u11+a.u22+a.u11*conj(a.u22)-a.u12*conj(a.u21)
|
||||
function Base.:*(a::SU3,b::SU3)
|
||||
|
||||
bu31 = (b.u12*conj(a.u23) - b.u13*conj(b.u22))
|
||||
bu32 = (b.u13*conj(b.u21) - b.u11*conj(b.u23))
|
||||
bu33 = (b.u11*conj(b.u22) - b.u12*conj(b.u21))
|
||||
|
||||
return SU3(a.u11*b.u11 + a.u12*b.u21 + a.u13*bu31,
|
||||
a.u11*b.u12 + a.u12*b.u22 + a.u13*bu32,
|
||||
a.u11*b.u13 + a.u12*b.u23 + a.u13*bu33,
|
||||
a.u21*b.u11 + a.u22*b.u21 + a.u23*bu31,
|
||||
a.u21*b.u12 + a.u22*b.u22 + a.u23*bu32,
|
||||
a.u21*b.u13 + a.u22*b.u23 + a.u23*bu33)
|
||||
end
|
||||
|
||||
28
src/Groups/Groups.jl
Normal file
28
src/Groups/Groups.jl
Normal file
|
|
@ -0,0 +1,28 @@
|
|||
###
|
||||
### "THE BEER-WARE LICENSE":
|
||||
### Alberto Ramos wrote this file. As long as you retain this
|
||||
### notice you can do whatever you want with this stuff. If we meet some
|
||||
### day, and you think this stuff is worth it, you can buy me a beer in
|
||||
### return. <alberto.ramos@cern.ch>
|
||||
###
|
||||
### file: Groups.jl
|
||||
### created: Sun Jul 11 18:02:16 2021
|
||||
###
|
||||
|
||||
|
||||
module Groups
|
||||
|
||||
abstract type Group end
|
||||
abstract type Algebra end
|
||||
|
||||
include("GroupSU2.jl")
|
||||
|
||||
export SU2, SU2alg, dag, normalize, inverse, tr, projalg, norm
|
||||
export dot, expm, exp
|
||||
|
||||
include("GroupSU3.jl")
|
||||
|
||||
|
||||
|
||||
|
||||
end # module
|
||||
24
src/LatticeGPU.jl
Normal file
24
src/LatticeGPU.jl
Normal file
|
|
@ -0,0 +1,24 @@
|
|||
module LatticeGPU
|
||||
|
||||
using CUDA
|
||||
|
||||
include("Groups/Groups.jl")
|
||||
|
||||
using .Groups
|
||||
export Group, Algebra
|
||||
export SU2, SU2alg, dag, normalize, inverse, tr, projalg, norm
|
||||
export dot, expm, exp
|
||||
|
||||
include("Space/Space.jl")
|
||||
|
||||
using .Space
|
||||
export SpaceParm, KernelParm
|
||||
export map2latt, up, dw, shift
|
||||
|
||||
|
||||
include("YM/YM.jl")
|
||||
|
||||
using .YM
|
||||
|
||||
|
||||
end # module
|
||||
133
src/Space/Space.jl
Normal file
133
src/Space/Space.jl
Normal file
|
|
@ -0,0 +1,133 @@
|
|||
###
|
||||
### "THE BEER-WARE LICENSE":
|
||||
### Alberto Ramos wrote this file. As long as you retain this
|
||||
### notice you can do whatever you want with this stuff. If we meet some
|
||||
### day, and you think this stuff is worth it, you can buy me a beer in
|
||||
### return. <alberto.ramos@cern.ch>
|
||||
###
|
||||
### file: Space.jl
|
||||
### created: Mon Jul 12 16:44:35 2021
|
||||
###
|
||||
|
||||
|
||||
module Space
|
||||
|
||||
struct SpaceParm{N,M}
|
||||
ndim::Int64
|
||||
iL::NTuple{N,Int64}
|
||||
npls::Int64
|
||||
plidx::NTuple{M, Tuple{Int64, Int64}}
|
||||
|
||||
function SpaceParm{N}(x, bt, c=(0.0,0.0)) where {N}
|
||||
M = convert(Int64, round(N*(N-1)/2))
|
||||
N == length(x) || throw(ArgumentError("Tuple of incorrect length for dimension $N"))
|
||||
|
||||
pls = Vector{Tuple{Int64, Int64}}()
|
||||
for i in 1:N
|
||||
for j in i+1:N
|
||||
push!(pls, (i,j))
|
||||
end
|
||||
end
|
||||
return new{N,M}(N, bt, c, x, M, tuple(pls...))
|
||||
end
|
||||
end
|
||||
export SpaceParm
|
||||
|
||||
struct KernelParm
|
||||
threads::Tuple{Int64,Int64,Int64}
|
||||
blocks::Tuple{Int64,Int64,Int64}
|
||||
end
|
||||
export KernelParm
|
||||
|
||||
@inline shift(p::CartesianIndex{4}, sh::CartesianIndex{4}, lp::SpaceParm{4}) = CartesianIndex(mod1(p[1]+sh[1], lp.iL[1]),
|
||||
mod1(p[2]+sh[2], lp.iL[2]),
|
||||
mod1(p[3]+sh[3], lp.iL[3]),
|
||||
mod1(p[4]+sh[4], lp.iL[4]))
|
||||
@inline shift(p::CartesianIndex{3}, sh::CartesianIndex{3}, lp::SpaceParm{3}) = CartesianIndex(mod1(p[1]+sh[1], lp.iL[1]),
|
||||
mod1(p[2]+sh[2], lp.iL[2]),
|
||||
mod1(p[3]+sh[3], lp.iL[3]))
|
||||
@inline shift(p::CartesianIndex{2}, sh::CartesianIndex{2}, lp::SpaceParm{2}) = CartesianIndex(mod1(p[1]+sh[1], lp.iL[1]),
|
||||
mod1(p[2]+sh[2], lp.iL[2]))
|
||||
@inline function dw(p::CartesianIndex{4}, id, lp::SpaceParm{4})
|
||||
|
||||
if (id == 1)
|
||||
s = CartesianIndex(mod1(p[1]-1, lp.iL[1]), p[2], p[3], p[4])
|
||||
elseif (id == 2)
|
||||
s = CartesianIndex(p[1], mod1(p[2]-1, lp.iL[2]), p[3], p[4])
|
||||
elseif (id == 3)
|
||||
s = CartesianIndex(p[1], p[2], mod1(p[3]-1, lp.iL[3]), p[4])
|
||||
elseif (id == 4)
|
||||
s = CartesianIndex(p[1], p[2], p[3], mod1(p[4]-1, lp.iL[4]))
|
||||
end
|
||||
|
||||
return s
|
||||
end
|
||||
|
||||
@inline function dw(p::CartesianIndex{3}, id, lp::SpaceParm{3})
|
||||
|
||||
if (id == 1)
|
||||
s = CartesianIndex(mod1(p[1]-1, lp.iL[1]), p[2], p[3])
|
||||
elseif (id == 2)
|
||||
s = CartesianIndex(p[1], mod1(p[2]-1, lp.iL[2]), p[3])
|
||||
elseif (id == 3)
|
||||
s = CartesianIndex(p[1], p[2], mod1(p[3]-1, lp.iL[3]))
|
||||
end
|
||||
|
||||
return s
|
||||
end
|
||||
|
||||
@inline function dw(p::CartesianIndex{2}, id, lp::SpaceParm{2})
|
||||
|
||||
if (id == 1)
|
||||
s = CartesianIndex(mod1(p[1]-1, lp.iL[1]), p[2])
|
||||
elseif (id == 2)
|
||||
s = CartesianIndex(p[1], mod1(p[2]-1, lp.iL[2]))
|
||||
end
|
||||
|
||||
return s
|
||||
end
|
||||
|
||||
@inline function up(p::CartesianIndex{4}, id, lp::SpaceParm{4})
|
||||
|
||||
if (id == 1)
|
||||
s = CartesianIndex(mod1(p[1]+1, lp.iL[1]), p[2], p[3], p[4])
|
||||
elseif (id == 2)
|
||||
s = CartesianIndex(p[1], mod1(p[2]+1, lp.iL[2]), p[3], p[4])
|
||||
elseif (id == 3)
|
||||
s = CartesianIndex(p[1], p[2], mod1(p[3]+1, lp.iL[3]), p[4])
|
||||
elseif (id == 4)
|
||||
s = CartesianIndex(p[1], p[2], p[3], mod1(p[4]+1, lp.iL[4]))
|
||||
end
|
||||
|
||||
return s
|
||||
end
|
||||
|
||||
@inline function up(p::CartesianIndex{3}, id, lp::SpaceParm{3})
|
||||
|
||||
if (id == 1)
|
||||
s = CartesianIndex(mod1(p[1]+1, lp.iL[1]), p[2], p[3])
|
||||
elseif (id == 2)
|
||||
s = CartesianIndex(p[1], mod1(p[2]+1, lp.iL[2]), p[3])
|
||||
elseif (id == 3)
|
||||
s = CartesianIndex(p[1], p[2], mod1(p[3]+1, lp.iL[3]))
|
||||
end
|
||||
|
||||
return s
|
||||
end
|
||||
|
||||
@inline function up(p::CartesianIndex{2}, id, lp::SpaceParm{2})
|
||||
|
||||
if (id == 1)
|
||||
s = CartesianIndex(mod1(p[1]+1, lp.iL[1]), p[2])
|
||||
elseif (id == 2)
|
||||
s = CartesianIndex(p[1], mod1(p[2]+1, lp.iL[2]))
|
||||
end
|
||||
|
||||
return s
|
||||
end
|
||||
|
||||
@inline map2latt(th::NTuple{3,Int64},bl::NTuple{3,Int64}) = CartesianIndex(th[1],bl[1],bl[2],bl[2])
|
||||
|
||||
export map2latt, up, dw, shift
|
||||
|
||||
end
|
||||
25
src/YM/YM.jl
Normal file
25
src/YM/YM.jl
Normal file
|
|
@ -0,0 +1,25 @@
|
|||
###
|
||||
### "THE BEER-WARE LICENSE":
|
||||
### Alberto Ramos wrote this file. As long as you retain this
|
||||
### notice you can do whatever you want with this stuff. If we meet some
|
||||
### day, and you think this stuff is worth it, you can buy me a beer in
|
||||
### return. <alberto.ramos@cern.ch>
|
||||
###
|
||||
### file: YM.jl
|
||||
### created: Mon Jul 12 16:23:51 2021
|
||||
###
|
||||
|
||||
|
||||
module YM
|
||||
|
||||
struct GaugeParm
|
||||
beta::Float64
|
||||
cG::Tuple{Float64,Float64}
|
||||
end
|
||||
export GaugeParm
|
||||
|
||||
include("YMact.jl")
|
||||
export krnl_plaq!
|
||||
|
||||
|
||||
end
|
||||
41
src/YM/YMact.jl
Normal file
41
src/YM/YMact.jl
Normal file
|
|
@ -0,0 +1,41 @@
|
|||
###
|
||||
### "THE BEER-WARE LICENSE":
|
||||
### Alberto Ramos wrote this file. As long as you retain this
|
||||
### notice you can do whatever you want with this stuff. If we meet some
|
||||
### day, and you think this stuff is worth it, you can buy me a beer in
|
||||
### return. <alberto.ramos@cern.ch>
|
||||
###
|
||||
### file: YMact.jl
|
||||
### created: Mon Jul 12 18:31:19 2021
|
||||
###
|
||||
|
||||
function krnl_plaq!(plx, U, ipl, lp::SpaceParm)
|
||||
|
||||
id1, id2 = lp.plidx(ipl)
|
||||
X = map2latt((CUDA.threadIdx().x,CUDA.threadIdx().y,CUDA.threadIdx().z),
|
||||
(CUDA.blockIdx().x,CUDA.blockIdx().y,CUDA.blockIdx().z))
|
||||
Xu1 = up(X, id1)
|
||||
Xu2 = up(X, id2)
|
||||
|
||||
plx[X] = tr(U[X, id1]*U[Xu1, id2] / (U[X, id2]*U[Xu2, id1]))
|
||||
|
||||
return nothing
|
||||
end
|
||||
|
||||
function krnl_plaq!(plx, U, lp::SpaceParm)
|
||||
|
||||
X = map2latt((CUDA.threadIdx().x,CUDA.threadIdx().y,CUDA.threadIdx().z),
|
||||
(CUDA.blockIdx().x,CUDA.blockIdx().y,CUDA.blockIdx().z))
|
||||
|
||||
plx[X] = 0.0
|
||||
for ipl in 1:lp.npls
|
||||
id1, id2 = lp.plidx(ipl)
|
||||
Xu1 = up(X, id1)
|
||||
Xu2 = up(X, id2)
|
||||
|
||||
plx[X] += tr(U[X, id1]*U[Xu1, id2] / (U[X, id2]*U[Xu2, id1]))
|
||||
end
|
||||
plx[X] = plx[X]/lp.npls
|
||||
|
||||
return nothing
|
||||
end
|
||||
Loading…
Add table
Add a link
Reference in a new issue