Added tests and support for U(1) group

src/Groups/GroupU1.jl

@@ -0,0 +1,84 @@
+###
+### "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:    GroupU1.jl
+### created: Tue Sep 21 09:33:44 2021
+###                               
+
+using CUDA, Random
+
+import Base.:*, Base.:+, Base.:-,Base.:/,Base.:\,Base.exp,Base.zero,Base.one
+import Random.rand
+struct U1{T} <: Group
+    t1::T
+    t2::T
+end
+U1(a::T)          where T <: AbstractFloat = U1{T}(a,zero(T))
+inverse(b::U1{T}) where T <: AbstractFloat = U1{T}(b.t1,-b.t2)
+dag(a::U1{T})    where T <: AbstractFloat = inverse(a)
+norm(a::U1{T})    where T <: AbstractFloat = sqrt(a.t1^2+a.t2^2)
+norm2(a::U1{T})   where T <: AbstractFloat = a.t1^2+a.t2^2
+tr(g::U1{T})      where T <: AbstractFloat = complex(a.t1)
+Base.one(::Type{U1{T}}) where T <: AbstractFloat = U1{T}(one(T), zero(T))
+function Random.rand(rng::AbstractRNG, ::Random.SamplerType{U1{T}}) where T <: AbstractFloat
+    r = randn(rng,T)
+    return U1{T}(CUDA.cos(r),CUDA.sin(r))
+end
+
+"""
+    function normalize(a::U1)
+
+Return a normalized element of `SU(2)`
+"""
+function normalize(a::U1{T}) where T <: AbstractFloat
+    dr = norm(a)
+    return U1{T}(a.t1/dr, a.t2/dr)
+end
+
+Base.:*(a::U1{T},b::U1{T}) where T <: AbstractFloat = U1{T}(a.t1*b.t1-a.t2*b.t2, a.t1*b.t2+a.t2*b.t1)
+Base.:/(a::U1{T},b::U1{T}) where T <: AbstractFloat = U1{T}(a.t1*b.t1+a.t2*b.t2, -a.t1*b.t2+a.t2*b.t1)
+Base.:\(a::U1{T},b::U1{T}) where T <: AbstractFloat = U1{T}(a.t1*b.t1+a.t2*b.t2, a.t1*b.t2-a.t2*b.t1)
+
+struct U1alg{T} <: Algebra
+    t::T
+end
+projalg(g::U1{T})             where T <: AbstractFloat = U1alg{T}(g.t2)
+dot(a::U1alg{T}, b::U1alg{T}) where T <: AbstractFloat = a.t*b.t
+norm(a::U1alg{T})             where T <: AbstractFloat = abs(a.t)
+norm2(a::U1alg{T})            where T <: AbstractFloat = a.t^2
+Base.zero(::Type{U1alg{T}})   where T <: AbstractFloat = U1alg{T}(zero(T))
+Random.rand(rng::AbstractRNG, ::Random.SamplerType{U1alg{T}}) where T <: AbstractFloat = U1alg{T}(randn(rng,T))
+
+Base.:+(a::U1alg{T})             where T <: AbstractFloat = U1alg{T}(a.t)
+Base.:-(a::U1alg{T})             where T <: AbstractFloat = U1alg{T}(-a.t)
+Base.:+(a::U1alg{T},b::U1alg{T}) where T <: AbstractFloat = U1alg{T}(a.t+b.t)
+Base.:-(a::U1alg{T},b::U1alg{T}) where T <: AbstractFloat = U1alg{T}(a.t-b.t)
+
+Base.:*(a::U1alg{T},b::Number)   where T <: AbstractFloat = U1alg{T}(a.t*b)
+Base.:*(b::Number,a::U1alg{T})   where T <: AbstractFloat = U1alg{T}(a.t*b)
+Base.:/(a::U1alg{T},b::Number)   where T <: AbstractFloat = U1alg{T}(a.t/b)
+
+isgroup(a::U1{T}) where T <: AbstractFloat = (abs(a.t) -1.0) < 1.0E-10
+
+"""
+    function Base.exp(a::U1alg, t::Number=1)
+
+Computes `exp(a)`
+"""
+Base.exp(a::U1alg{T}) where T <: AbstractFloat = U1{T}(CUDA.cos(a.t), CUDA.sin(a.t))
+Base.exp(a::U1alg{T}, t::T) where T <: AbstractFloat  = U1{T}(CUDA.cos(t*a.t), CUDA.sin(t*a.t))
+
+"""
+    function expm(g::U1, a::U1alg; t=1)
+
+Computes `exp(a)*g`
+
+"""
+expm(g::U1{T}, a::U1alg{T}) where T <: AbstractFloat = U1{T}(CUDA.cos(a.t), CUDA.sin(a.t))*g
+expm(g::U1{T}, a::U1alg{T}, t::T) where T <: AbstractFloat = U1{T}(CUDA.cos(t*a.t), CUDA.sin(t*a.t))*g
+
+export U1, U1alg, inverse, dag, tr, projalg, expm, exp, norm, norm2, isgroup

src/Groups/Groups.jl

@@ -23,6 +23,10 @@ export SU2, SU2alg
 include("GroupSU3.jl")
 export SU3, SU3alg
 
+include("GroupU1.jl")
+export U1, U1alg
+
+
 export dot, expm, exp, dag, normalize, inverse, tr, projalg, norm, norm2, isgroup
 
 

src/LatticeGPU.jl

@@ -16,7 +16,7 @@ include("Groups/Groups.jl")
 
 using .Groups
 export Group, Algebra
-export SU2, SU2alg, SU3, SU3alg
+export SU2, SU2alg, SU3, SU3alg, U1, U1alg
 export dot, expm, exp, dag, normalize, inverse, tr, projalg, norm, norm2, isgroup
 
 include("Space/Space.jl")

tests/testU1.jl

@@ -0,0 +1,40 @@
+using LinearAlgebra, Random
+
+import Pkg
+#Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU")
+Pkg.activate("/home/alberto/code/julia/LatticeGPU")
+using LatticeGPU
+
+
+T = Float32
+
+b = rand(U1{T})
+println(b)
+
+ba = rand(U1alg{T})
+println("Ba:        ", ba)
+b = exp(ba)
+println("B:         ", b)
+c = exp(ba, convert(T,-1))
+println(typeof(norm2(ba)))
+d = b*c
+println("Test:      ", d)
+
+c = inverse(b)
+println("Inverse B: ", c)
+
+d = b*c
+println("Test:      ", d)
+
+println("B:         ", b)
+println("Ba:        ", ba)
+b = expm(b, ba, convert(T,-1))
+println("Test:      ", b)
+
+
+Ma = Array{U1{T}}(undef, 2)
+rand!(Ma)
+println(Ma)
+
+fill!(Ma, one(eltype(Ma)))
+println(Ma)

tests/test_SU2.jl

@@ -0,0 +1,31 @@
+using LinearAlgebra, Random
+
+import Pkg
+#Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU")
+Pkg.activate("/home/alberto/code/julia/LatticeGPU")
+using LatticeGPU
+
+
+T = Float32
+
+b = rand(SU2{T})
+println(b)
+
+ba = rand(SU2alg{T})
+println("Ba:        ", ba)
+b = exp(ba)
+println("B:         ", b)
+println(typeof(norm2(ba)))
+
+c = inverse(b)
+println("Inverse B: ", c)
+
+d = b*c
+println("Test:      ", d)
+
+Ma = Array{SU2{T}}(undef, 100)
+rand!(Ma)
+println(Ma)
+
+fill!(Ma, one(eltype(Ma)))
+println(Ma)

tests/test_space.jl

@@ -0,0 +1,41 @@
+
+import Pkg
+#Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU")
+Pkg.activate("/home/alberto/code/julia/LatticeGPU")
+using LatticeGPU, BenchmarkTools
+
+
+lp = SpaceParm{4}((8,12,6,6), (8,2,2,3))
+
+function test_point(pt::NTuple{2,Int64}, lp::SpaceParm)
+    ok = true
+    println("Global point: ", global_point(pt, lp))
+    for id in 1:lp.ndim
+        ua, ub = up(pt, id, lp)
+        println("  - UP in id $id: ", global_point((ua,ub), lp))
+        
+        da, db = dw(pt, id, lp)
+        println("  - DW in id $id: ", global_point((da,db), lp), "\n")
+        
+        ua2, ub2, da2, db2 = updw(pt, id, lp)
+        ok = ok && (ua == ua2)
+        ok = ok && (ub == ub2)
+        ok = ok && (da == da2)
+        ok = ok && (db == db2)
+    end
+    return ok
+end
+
+global ok = true
+for i in 1:lp.bsz, j in 1:lp.rsz
+    global ok = ok && test_point((i,j), lp)
+end
+
+if ok
+    println("ALL tests passed")
+else
+    println("ERROR in test")
+end
+    
+println(lp)
+

tests/test_su3.jl

@@ -0,0 +1,131 @@
+using CUDA, LinearAlgebra
+
+import Pkg
+#Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU")
+Pkg.activate("/home/alberto/code/julia/LatticeGPU")
+using LatticeGPU
+
+function g2mat(g::SU3)
+
+    M = Array{ComplexF64, 2}(undef, 3,3)
+
+    M[1,1] = g.u11
+    M[1,2] = g.u12
+    M[1,3] = g.u13
+
+    M[2,1] = g.u21
+    M[2,2] = g.u22
+    M[2,3] = g.u23
+
+    M[3,1] = conj(g.u12*g.u23 - g.u13*g.u22)
+    M[3,2] = conj(g.u13*g.u21 - g.u11*g.u23)
+    M[3,3] = conj(g.u11*g.u22 - g.u12*g.u21)
+
+    return M
+end
+    
+
+
+# 0.40284714488721746 + 0.2704272209422031im -0.029482825024553627 - 0.8247329455356851im 0.28771631112777535 + 0.027366985901323956im; -0.08478364480998268 + 0.8226014762207954im -0.4790638417896126 + 0.24301903735299646im -0.022591091614522323 + 0.16452285690920823im; 0.28083864951126214 + 0.04302898862961919im 0.0066864552013863165 - 0.17418727240313508im -0.939634663641523 + 0.07732362776719631im
+
+T  = Float64
+a = rand(SU3alg{T})
+println("Random algebra: ", a)
+g1 = exp(a, 0.2)
+g2 = exp(a, -0.2)
+
+g = expm(g1, a, -0.2)
+println(g)
+
+
+
+
+
+M = [0.40284714488721746 + 0.2704272209422031im -0.029482825024553627 - 0.8247329455356851im 0.28771631112777535 + 0.027366985901323956im; -0.08478364480998268 + 0.8226014762207954im -0.4790638417896126 + 0.24301903735299646im -0.022591091614522323 + 0.16452285690920823im; 0.28083864951126214 + 0.04302898862961919im 0.0066864552013863165 - 0.17418727240313508im -0.939634663641523 + 0.07732362776719631im]
+
+println(det(M))
+
+
+g1 = SU3(M[1,1],M[1,2],M[1,3], M[2,1],M[2,2],M[2,3])
+println("dir: ", g1)
+g2 = exp(a)
+println("exp: ", g2)
+println("dif: ", g2mat(g1)-g2mat(g2))
+
+g3 = g1/g2
+println(g3)
+
+println("END")
+
+ftest(g::Group) = LatticeGPU.tr(g)
+
+#println(" ## SU(2)")
+#asu2 = SU2alg(0.23, 1.23, -0.34)
+#gsu2 = exp(asu2) 
+#
+#eps = 0.001
+#h   = SU2alg(eps,0.0,0.0)
+#fp = ftest(exp(h)*gsu2)
+#fm = ftest(exp(h,-1.0)*gsu2)
+#println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+#h   = SU2alg(0.0,eps,0.0)
+#fp = ftest(exp(h)*gsu2)
+#fm = ftest(exp(h,-1.0)*gsu2)
+#println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+#h   = SU2alg(0.0,0.0,eps)
+#fp = ftest(exp(h)*gsu2)
+#fm = ftest(exp(h,-1.0)*gsu2)
+#println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+#println("Exact     derivative: ", -projalg(gsu2))
+
+
+
+println("\n\n ## SU(3)")
+asu3 = SU3alg{T}(0.23, 1.23, -0.34, 2.34, -0.23, 0.23, -1.34, 1.34)
+gsu3 = exp(asu3) 
+
+eps = 0.001
+h   = SU3alg{T}(eps,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
+fp = ftest(exp(h)*gsu3)
+fm = ftest(exp(h,-1.0)*gsu3)
+println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+h   = SU3alg{T}(0.0,eps,0.0,0.0,0.0,0.0,0.0,0.0)
+fp = ftest(exp(h)*gsu3)
+fm = ftest(exp(h,-1.0)*gsu3)
+println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+h   = SU3alg{T}(0.0,0.0,eps,0.0,0.0,0.0,0.0,0.0)
+fp = ftest(exp(h)*gsu3)
+fm = ftest(exp(h,-1.0)*gsu3)
+println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+h   = SU3alg{T}(0.0,0.0,0.0,eps,0.0,0.0,0.0,0.0)
+fp = ftest(exp(h)*gsu3)
+fm = ftest(exp(h,-1.0)*gsu3)
+println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+h   = SU3alg{T}(0.0,0.0,0.0,0.0,eps,0.0,0.0,0.0)
+fp = ftest(exp(h)*gsu3)
+fm = ftest(exp(h,-1.0)*gsu3)
+println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+h   = SU3alg{T}(0.0,0.0,0.0,0.0,0.0,eps,0.0,0.0)
+fp = ftest(exp(h)*gsu3)
+fm = ftest(exp(h,-1.0)*gsu3)
+println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+h   = SU3alg{T}(0.0,0.0,0.0,0.0,0.0,0.0,eps,0.0)
+fp = ftest(exp(h)*gsu3)
+fm = ftest(exp(h,-1.0)*gsu3)
+println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+h   = SU3alg{T}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,eps)
+fp = ftest(exp(h)*gsu3)
+fm = ftest(exp(h,-1.0)*gsu3)
+println("Numerical derivative: ", (fp-fm)/(2.0*eps))
+println("Exact     derivative: ", -projalg(gsu3))
+
+println("\n # Mutiplications")
+
+g1 = exp(SU3alg{T}(0.23, 1.23, -0.34, 2.34, -0.23, 0.23, -1.34, 1.34))
+g2 = exp(SU3alg{T}(1.23, -0.23, -0.14, 0.4, -1.23, -0.8, -0.34, 0.34))
+
+a = g1/(g2*g1)
+b = g2*a
+println("b is one: ", b)
+
+