diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl index f5a5f4a..78ab6ab 100644 --- a/src/Fields/Fields.jl +++ b/src/Fields/Fields.jl @@ -19,7 +19,9 @@ vector_field(::Type{T}, lp::SpaceParm) where {T} = CuArray{T, 3}(undef, lp.b scalar_field(::Type{T}, lp::SpaceParm) where {T} = CuArray{T, 2}(undef, lp.bsz, lp.rsz) nscalar_field(::Type{T}, n, lp::SpaceParm) where {T} = CuArray{T, 3}(undef, lp.bsz, n, lp.rsz) -export vector_field, scalar_field, nscalar_field +scalar_field_point(::Type{T}, lp::SpaceParm{N,M,D}) where {T,N,M,D} = CuArray{T, N}(undef, lp.iL...) + +export vector_field, scalar_field, nscalar_field, scalar_field_point end diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl index c920257..12c7661 100644 --- a/src/LatticeGPU.jl +++ b/src/LatticeGPU.jl @@ -27,7 +27,7 @@ export up, dw, updw, global_point include("Fields/Fields.jl") using .Fields -export vector_field, scalar_field, nscalar_field +export vector_field, scalar_field, nscalar_field, scalar_field_point include("MD/MD.jl") using .MD diff --git a/src/Space/Space.jl b/src/Space/Space.jl index f029f8b..b13535c 100644 --- a/src/Space/Space.jl +++ b/src/Space/Space.jl @@ -169,25 +169,67 @@ Given a point `x` with index `p`, this routine returns the index of the points return bu, ru, bd, rd end -@inline function point_coord(p::NTuple{2,Int64}, lp::SpaceParm) +@inline cntb(nb, id::Int64, lp::SpaceParm) = mod(div(nb-1,lp.blkS[id]),lp.blk[id]) +@inline cntr(nr, id::Int64, lp::SpaceParm) = mod(div(nr-1,lp.rbkS[id]),lp.rbk[id]) +@inline cnt(nb, nr, id::Int64, lp::SpaceParm) = 1 + cntb(nb,id,lp) + cntr(nr,id,lp)*lp.blk[id] - @inline cntb(nb, id::Int64, lp::SpaceParm) = mod(div(nb-1,lp.blkS[id]),lp.blk[id]) - @inline cntr(nr, id::Int64, lp::SpaceParm) = mod(div(nr-1,lp.rbkS[id]),lp.rbk[id]) +@inline function point_coord(p::NTuple{2,Int64}, lp::SpaceParm{2,M,D}) where {M,D} - @inline cnt(nb, nr, id::Int64, lp::SpaceParm) = 1 + cntb(nb,id,lp) + cntr(nr,id,lp)*lp.blk[id] - - pt = ntuple(i -> cnt(p[1], p[2], i, lp), lp.ndim) - return pt + i1 = cnt(p[1], p[2], 1, lp) + i2 = cnt(p[1], p[2], 2, lp) + +# pt = ntuple(i -> cnt(p[1], p[2], i, lp), lp.ndim) + return CartesianIndex{4}(i1,i2) end -@inline function point_time(p::NTuple{2,Int64}, lp::SpaceParm) +@inline function point_coord(p::NTuple{2,Int64}, lp::SpaceParm{3,M,D}) where {M,D} - @inline cntb(nb, id::Int64, lp::SpaceParm) = mod(div(nb-1,lp.blkS[id]),lp.blk[id]) - @inline cntr(nr, id::Int64, lp::SpaceParm) = mod(div(nr-1,lp.rbkS[id]),lp.rbk[id]) + i1 = cnt(p[1], p[2], 1, lp) + i2 = cnt(p[1], p[2], 2, lp) + i3 = cnt(p[1], p[2], 3, lp) - @inline cnt(nb, nr, id::Int64, lp::SpaceParm) = 1 + cntb(nb,id,lp) + cntr(nr,id,lp)*lp.blk[id] - - return cnt(p[1], p[2], 1, lp) +# pt = ntuple(i -> cnt(p[1], p[2], i, lp), lp.ndim) + return CartesianIndex{4}(i1,i2,i3) +end + +@inline function point_coord(p::NTuple{2,Int64}, lp::SpaceParm{4,M,D}) where {M,D} + + i1 = cnt(p[1], p[2], 1, lp) + i2 = cnt(p[1], p[2], 2, lp) + i3 = cnt(p[1], p[2], 3, lp) + i4 = cnt(p[1], p[2], 4, lp) + +# pt = ntuple(i -> cnt(p[1], p[2], i, lp), lp.ndim) + return CartesianIndex{4}(i1,i2,i3,i4) +end + +@inline function point_coord(p::NTuple{2,Int64}, lp::SpaceParm{5,M,D}) where {M,D} + + i1 = cnt(p[1], p[2], 1, lp) + i2 = cnt(p[1], p[2], 2, lp) + i3 = cnt(p[1], p[2], 3, lp) + i4 = cnt(p[1], p[2], 4, lp) + i5 = cnt(p[1], p[2], 5, lp) + +# pt = ntuple(i -> cnt(p[1], p[2], i, lp), lp.ndim) + return CartesianIndex{4}(i1,i2,i3,i4,i5) +end + +@inline function point_coord(p::NTuple{2,Int64}, lp::SpaceParm{6,M,D}) where {M,D} + + i1 = cnt(p[1], p[2], 1, lp) + i2 = cnt(p[1], p[2], 2, lp) + i3 = cnt(p[1], p[2], 3, lp) + i4 = cnt(p[1], p[2], 4, lp) + i5 = cnt(p[1], p[2], 5, lp) + i6 = cnt(p[1], p[2], 6, lp) + +# pt = ntuple(i -> cnt(p[1], p[2], i, lp), lp.ndim) + return CartesianIndex{4}(i1,i2,i3,i4,i5,i6) +end + +@inline function point_time(p::NTuple{2,Int64}, lp::SpaceParm{N,M,D}) where {N,M,D} + return cnt(p[1], p[2], N, lp) end diff --git a/src/YM/YM.jl b/src/YM/YM.jl index da435d7..e350843 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -66,8 +66,8 @@ struct YMworkspace{T} mm = vector_field(SU3alg{T}, lp) u1 = vector_field(SU3{T}, lp) end - cs = scalar_field(Complex{T}, lp) - rs = scalar_field(T, lp) + cs = scalar_field_point(Complex{T}, lp) + rs = scalar_field_point(T, lp) end return new{T}(GRP,ALG,T,f1, f2, mm, u1, cs, rs) diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl index a2d935b..bab4918 100644 --- a/src/YM/YMact.jl +++ b/src/YM/YMact.jl @@ -15,7 +15,7 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,D}) wher Ush = @cuStaticSharedMem(T, (D,2)) - plx[b,r] = zero(plx[b,r]) + S = zero(eltype(plx)) for id1 in 1:N-1 bu1, ru1 = up((b, r), id1, lp) Ush[b,1] = U[b,id1,r] @@ -85,10 +85,13 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,D}) wher g2 = Ush[b,2]\Ush[b,1] - plx[b,r] += c0*tr(g2*ga/gb) + c1*( tr(g2*h2/gb) + tr(g2*ga/h3)) + S += c0*tr(g2*ga/gb) + c1*( tr(g2*h2/gb) + tr(g2*ga/h3)) end end + I = point_coord((b,r), lp) + plx[I] = S + return nothing end @@ -98,7 +101,7 @@ function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,D}) where {T,N,M Ush = @cuStaticSharedMem(T, (D,2)) - plx[b,r] = zero(plx[b,r]) + S = zero(eltype(plx)) for id1 in 1:N-1 bu1, ru1 = up((b, r), id1, lp) Ush[b,1] = U[b,id1,r] @@ -119,10 +122,13 @@ function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,D}) where {T,N,M gt2 = U[bu2,id1,ru2] end - plx[b,r] += tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) + S += tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) end end + I = point_coord((b,r), lp) + plx[I] = S + return nothing end diff --git a/src/main/times.jl b/src/main/times.jl index 1a999e4..f991e5b 100644 --- a/src/main/times.jl +++ b/src/main/times.jl @@ -39,6 +39,7 @@ println("\n## WILSON ACTION/FLOW TIMES") gp = GaugeParm{PREC}(6.0, 1.0, (0.0,0.0), 3) println("Gauge Parameters: ", gp) + println("Initial Action: ") @time S = gauge_action(U, lp, gp, ymws)