diff --git a/src/Dirac/Diracfields.jl b/src/Dirac/Diracfields.jl
index c38aefe1437aa8cdd236be5e34a2e1755fbe652c..9f4c46cd683ad03b05232c66b28fea8a8a8be44f 100644
--- a/src/Dirac/Diracfields.jl
+++ b/src/Dirac/Diracfields.jl
@@ -131,7 +131,7 @@ end
Randomizes the SU2fund / SU3fund fermion field. If the argument t is present, it only randomizes that time-slice.
"""
-function pfrandomize!(f::AbstractArray{Spinor{4, SU3fund{T,2}}}, lp::SpaceParm{4,6,BC_PERIODIC,D}, t::Int64 = 0) where {T,D}
+function pfrandomize!(f::AbstractArray{Spinor{4, SU3fund{T,N}}}, lp::SpaceParm{4,6,BC_PERIODIC,D}, t::Int64 = 0) where {T,D,N}
@timeit "Randomize pseudofermion field" begin
p = ntuple(i->CUDA.randn(Complex{T}, lp.bsz, 3, lp.rsz),4) # complex generation not suported for Julia 1.5.4
@@ -143,7 +143,7 @@ function pfrandomize!(f::AbstractArray{Spinor{4, SU3fund{T,2}}}, lp::SpaceParm{4
return nothing
end
-function pfrandomize!(f::AbstractArray{Spinor{4, SU3fund{T,2}}}, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D},SpaceParm{4,6,BC_OPEN,D}}, t::Int64 = 0) where {T,D}
+function pfrandomize!(f::AbstractArray{Spinor{4, SU3fund{T,N}}}, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D},SpaceParm{4,6,BC_OPEN,D}}, t::Int64 = 0) where {T,D,N}
@timeit "Randomize pseudofermion field" begin
p = ntuple(i->CUDA.randn(Complex{T}, lp.bsz, 3, lp.rsz),4) # complex generation not suported for Julia 1.5.4
@@ -156,20 +156,20 @@ function pfrandomize!(f::AbstractArray{Spinor{4, SU3fund{T,2}}}, lp::Union{Space
return nothing
end
-function krnl_assign_pf_su3!(f::AbstractArray{Spinor{4, SU3fund{T,2}}}, p , lp::SpaceParm, t::Int64) where {T}
+function krnl_assign_pf_su3!(f::AbstractArray{Spinor{4, SU3fund{T,N}}}, p, lp::SpaceParm, t::Int64) where {T,N}
@inbounds begin
b = Int64(CUDA.threadIdx().x)
r = Int64(CUDA.blockIdx().x)
if t == 0
- f[b,r] = Spinor(map(x->SU3fund{Float64,2}(Series{Complex{T},2}((x[b,1,r],zero(T))),
- Series{Complex{T},2}((x[b,2,r],zero(T))),
- Series{Complex{T},2}((x[b,3,r],zero(T)))),p))
+ f[b,r] = Spinor(map(x->SU3fund{Float64,N}(Series{Complex{T},N}(ntuple(i -> (i==1)*x[b,1,r],N)),
+ Series{Complex{T},N}(ntuple(i -> (i==1)*x[b,2,r],N)),
+ Series{Complex{T},N}(ntuple(i -> (i==1)*x[b,2,r],N))),p))
elseif point_time((b,r),lp) == t
- f[b,r] = Spinor(map(x->SU3fund{Float64,2}(Series{Complex{T},2}((x[b,1,r],zero(T))),
- Series{Complex{T},2}((x[b,2,r],zero(T))),
- Series{Complex{T},2}((x[b,3,r],zero(T)))),p))
+ f[b,r] = Spinor(map(x->SU3fund{Float64,N}(Series{Complex{T},N}(ntuple(i -> (i==1)*x[b,1,r],N)),
+ Series{Complex{T},N}(ntuple(i -> (i==1)*x[b,2,r],N)),
+ Series{Complex{T},N}(ntuple(i -> (i==1)*x[b,2,r],N))),p))
else
f[b,r] = 0.0*f[b,r]
end
diff --git a/src/Groups/SU3Types.jl b/src/Groups/SU3Types.jl
index cf78be3cfdff0de1f03667a825f9f9ebcf7d0128..260a9198c351a5d523f6a4ebd99979fbd0eceeda 100644
--- a/src/Groups/SU3Types.jl
+++ b/src/Groups/SU3Types.jl
@@ -86,7 +86,10 @@ struct SU3fund{T,N}
t3::Series{Complex{T},N}
end
Base.zero(::Type{SU3fund{T,1}}) where T <: AbstractFloat = SU3fund{T,1}(zero(T),zero(T),zero(T))
-Base.zero(::Type{SU3fund{T,2}}) where T <: AbstractFloat = SU3fund{T,2}(Series{Complex{T},2}((zero(Complex{T}),zero(Complex{T}))),Series{Complex{T},2}((zero(Complex{T}),zero(Complex{T}))),Series{Complex{T},2}((zero(Complex{T}),zero(Complex{T}))))
+# Base.zero(::Type{SU3fund{T,2}}) where T <: AbstractFloat = SU3fund{T,2}(Series{Complex{T},2}((zero(Complex{T}),zero(Complex{T}))),Series{Complex{T},2}((zero(Complex{T}),zero(Complex{T}))),Series{Complex{T},2}((zero(Complex{T}),zero(Complex{T}))))
+Base.zero(::Type{SU3fund{T,N}}) where {T <: AbstractFloat,N} = SU3fund{T,N}(Series{Complex{T},N}(ntuple(i -> zero(Complex{T}),N)),
+ Series{Complex{T},N}(ntuple(i -> zero(Complex{T}),N)),
+ Series{Complex{T},N}(ntuple(i -> zero(Complex{T}),N)))
Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU3fund{T,1}}) where T <: AbstractFloat = SU3fund{T,1}(complex(randn(rng,T),randn(rng,T)),
complex(randn(rng,T),randn(rng,T)),
complex(randn(rng,T),randn(rng,T)))
diff --git a/src/Solvers/CG.jl b/src/Solvers/CG.jl
index 84061a7728fc9f7e4199f0906ecb8bda38ada592..822132732c186d3c2bb5dc4593cdf3457433fd16 100644
--- a/src/Solvers/CG.jl
+++ b/src/Solvers/CG.jl
@@ -33,7 +33,7 @@ end
Solves the linear equation `Ax = si`
"""
-function CG!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T}, maxiter::Int64 = 10, tol=1.0) where {T}
+function CG!(si, U, A, dpar::DiracParam{T,R,N} where R, lp::SpaceParm, dws::DiracWorkspace{T}, maxiter::Int64 = 10, tol=1.0) where {T,N}
dws.sr .= si
dws.sp .= si
@@ -44,7 +44,7 @@ function CG!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T},
tol = abs((tol*norm).c[1])
iterations = 0
- sumf = scalar_field(Series{Complex{T},2}, lp)
+ sumf = scalar_field(Series{Complex{T},N}, lp)
niter = 0
for i in 1:maxiter
diff --git a/src/Solvers/Propagators.jl b/src/Solvers/Propagators.jl
index b4bc9dd3f307710e768fabd6c0c413af78401fc0..726696989617bea589b04a4af726725c5a9362f1 100644
--- a/src/Solvers/Propagators.jl
+++ b/src/Solvers/Propagators.jl
@@ -18,7 +18,7 @@ Saves the fermionic progapator in pro for a source at point `y` with color `c` a
by a random source in spin and color at t = `time`. Returns the number of iterations.
"""
-function propagator!(pro, U, dpar::DiracParam{T,SU3fund{T,2}}, dws::DiracWorkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, y::NTuple{4,Int64}, c::Int64, s::Int64) where {T}
+function propagator!(pro, U, dpar::DiracParam{T,SU3fund{T,N},N}, dws::DiracWorkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, y::NTuple{4,Int64}, c::Int64, s::Int64) where {T,N}
function krnlg5!(src)
b=Int64(CUDA.threadIdx().x)
@@ -29,9 +29,9 @@ function propagator!(pro, U, dpar::DiracParam{T,SU3fund{T,2}}, dws::DiracWorkspa
@timeit "Propagator computation" begin
- fill!(dws.sp,zero(Spinor{4,SU3fund{T,2}}))
+ fill!(dws.sp,zero(Spinor{4,SU3fund{T,N}}))
- CUDA.@allowscalar dws.sp[point_index(CartesianIndex{lp.ndim}(y),lp)...] = Spinor{4,SU3fund{Float64,2}}(ntuple(i -> (i==s)*SU3fund{T,2}(ntuple(j -> (j==c)*Series{Complex{T},2}((complex(1.0,0.0),0.0)),3)...),4))
+ CUDA.@allowscalar dws.sp[point_index(CartesianIndex{lp.ndim}(y),lp)...] = Spinor{4,SU3fund{Float64,N}}(ntuple(i -> (i==s)*SU3fund{T,N}(ntuple(j -> (j==c)*Series{Complex{T},N}((complex(1.0,0.0),zeros(T,N-1)...)),3)...),4))
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp)
@@ -45,7 +45,7 @@ function propagator!(pro, U, dpar::DiracParam{T,SU3fund{T,2}}, dws::DiracWorkspa
return niter
end
-function propagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, time::Int64) where {T}
+function propagator!(pro, U, dpar::DiracParam{T,SU3fund{T,N},N}, dws::DiracWorkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, time::Int64) where {T,N}
function krnlg5!(src)
b=Int64(CUDA.threadIdx().x)
@@ -55,7 +55,7 @@ function propagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Space
end
@timeit "Propagator computation" begin
- fill!(dws.sp,zero(Spinor{4,SU3fund{T,2}}))
+ fill!(dws.sp,zero(Spinor{4,SU3fund{T,N}}))
pfrandomize!(dws.sp,lp,time)
@@ -79,7 +79,7 @@ Saves the propagator from the t=0 boundary to the bulk for the SF boundary condi
For the propagator from T to the bulk, use the function Tbndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64). Returns the number of iterations.
"""
-function bndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D}
+function bndpropagator!(pro, U, dpar::DiracParam{T,SU3fund{T,N},N}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D,N}
function krnlg5!(src)
b=Int64(CUDA.threadIdx().x)
@@ -94,7 +94,7 @@ function bndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Sp
if (point_time((b,r),lp) == 2)
bd4, rd4 = dw((b,r), 4, lp)
- src[b,r] = gdagpmul(Pgamma{4,1},U[bd4,4,rd4],Spinor{4,SU3fund{Float64,2}}(ntuple(i -> (i==s)*SU3fund{T,2}(ntuple(j -> (j==c)*Series{Complex{Float64},2}((complex(1.0,0.0),0.0)),3)...),4)))/2
+ src[b,r] = gdagpmul(Pgamma{4,1},U[bd4,4,rd4],Spinor{4,SU3fund{Float64,N}}(ntuple(i -> (i==s)*SU3fund{T,N}(ntuple(j -> (j==c)*Series{Complex{Float64},N}(ntuple(k -> (k==1)*complex(1.0,0.0),N)),3)...),4)))/2
end
return nothing
@@ -102,7 +102,7 @@ function bndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Sp
@timeit "Propagator computation" begin
SF_bndfix!(pro,lp)
- fill!(dws.sp,zero(Spinor{4,SU3fund{T,2}}))
+ fill!(dws.sp,zero(Spinor{4,SU3fund{T,N}}))
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_bndsrc!(dws.sp, U, lp, c, s)
@@ -128,7 +128,7 @@ Returns the propagator from the t=T boundary to the bulk for the SF boundary con
For the propagator from t=0 to the bulk, use the function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64). Returns the number of iterations.
"""
-function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D}
+function Tbndpropagator!(pro, U, dpar::DiracParam{T,SU3fund{T,N},N}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D,N}
function krnlg5!(src)
b=Int64(CUDA.threadIdx().x)
@@ -142,14 +142,14 @@ function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::S
r=Int64(CUDA.blockIdx().x)
if (point_time((b,r),lp) == lp.iL[end])
- src[b,r] = gpmul(Pgamma{4,-1},U[b,4,r],Spinor{4,SU3fund{Float64,2}}(ntuple(i -> (i==s)*SU3fund{T,2}(ntuple(j -> (j==c)*Series{Complex{T},2}((complex(1.0,0.0),0.0)),3)...),4)))/2
+ src[b,r] = gpmul(Pgamma{4,-1},U[b,4,r],Spinor{4,SU3fund{Float64,N}}(ntuple(i -> (i==s)*SU3fund{T,N}(ntuple(j -> (j==c)*Series{Complex{T},N}(ntuple(k -> (k==1)*complex(1.0,0.0),N)),3)...),4)))/2
end
return nothing
end
@timeit "Propagator computation" begin
- fill!(dws.sp,zero(Spinor{4,SU3fund{T,2}}))
+ fill!(dws.sp,zero(Spinor{4,SU3fund{T,N}}))
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_bndsrc!(dws.sp, U, lp, c, s)