From cd8a9d4f7cb8af1f5c010afd4c8dbc9d7041d625 Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Tue, 24 Jan 2023 17:57:42 +0100 Subject: [PATCH] Removed Shared memory use. It seems that the use of shared memory hits some compiler bug. See https://discourse.julialang.org/t/bug-with-julia-1-7-1-and-cuda-3-3/74052 --- src/YM/YMact.jl | 279 ++++++++++++----------------------------------- src/YM/YMflow.jl | 94 +++++----------- 2 files changed, 93 insertions(+), 280 deletions(-) diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl index f06e7ce..5cf7063 100644 --- a/src/YM/YMact.jl +++ b/src/YM/YMact.jl @@ -15,101 +15,57 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::NTuple{NB,T}, cG, zt r = Int64(CUDA.blockIdx().x) it = point_time((b, r), lp) - Ush = @cuStaticSharedMem(T, (D,2)) - ipl = 0 S = zero(eltype(plx)) @inbounds begin for id1 in N:-1:1 bu1, ru1 = up((b, r), id1, lp) SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1==N) - Ush[b,1] = U[b,id1,r] for id2 = 1:id1-1 bu2, ru2 = up((b, r), id2, lp) - Ush[b,2] = U[b,id2,r] - sync_threads() ipl = ipl + 1 # 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] + if SFBC && (it == lp.iL[end]-1) + gb = Ubnd[id2] else - if SFBC && (it == lp.iL[end]-1) - gb = Ubnd[id2] - else - gb = U[b2,id2,r2] - end + 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 + h2 = (U[b1,id1,r1]*gb)/U[b2,id1,r2] # 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] + (b3, r3) = up((b1,r1), id1, lp) + if SFBC && (it == lp.iL[end]) + gc = Ubnd[id2] else - if SFBC && (it == lp.iL[end]) - gc = Ubnd[id2] - else - gc = U[b2,id2,r2] - end + gc = U[b3,id2,r3] end - h3 = (ga*gb)/gc + h3 = (U[b1,id2,r1]*U[b2,id1,r2])/gc # END staples - if ru2 == r - gb = Ush[bu2,1] + if SFBC && (it == lp.iL[end]) + ga = Ubnd[id2] else - gb = U[bu2,id1,ru2] - end - if ru1 == r - ga = Ush[bu1,2] - else - if SFBC && (it == lp.iL[end]) - ga = Ubnd[id2] - else - ga = U[bu1,id2,ru1] - end + ga = U[bu1,id2,ru1] end - g2 = Ush[b,2]\Ush[b,1] + g2 = U[b,id2,r]\U[b,id1,r] if (it == lp.iL[end]) && SFBC - S += cG*(c0*tr(g2*ga/gb) + (3*c1/2)*tr(g2*ga/h3)) + S += cG*(c0*tr(g2*ga/U[bu2,id1,ru2]) + (3*c1/2)*tr(g2*ga/h3)) elseif (it == 1) && SFBC - S += cG*(c0*tr(g2*ga/gb) + (3*c1/2)*tr(g2*ga/h3)) + c1*tr(g2*h2/gb) + S += cG*(c0*tr(g2*ga/U[bu2,id1,ru2]) + (3*c1/2)*tr(g2*ga/h3)) + c1*tr(g2*h2/U[bu2,id1,ru2]) else - S += ztw[ipl]*c0*tr(g2*ga/gb) + - (ztw[ipl]^2*c1)*( tr(g2*h2/gb) + tr(g2*ga/h3)) + S += ztw[ipl]*c0*tr(g2*ga/U[bu2,id1,ru2]) + + (ztw[ipl]^2*c1)*( tr(g2*h2/U[bu2,id1,ru2]) + tr(g2*ga/h3)) end end @@ -130,42 +86,27 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd, cG, ztw, lp::SpaceParm{N,M,B r = Int64(CUDA.blockIdx().x) it = point_time((b, r), lp) - Ush = @cuStaticSharedMem(T, (D,2)) - S = zero(eltype(plx)) ipl = 0 for id1 in N:-1:1 bu1, ru1 = up((b, r), id1, lp) - Ush[b,1] = U[b,id1,r] - SFBND = ( ( (B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && ( (it == 1) || (it == lp.iL[end])) ) && (id1 == N) for id2 = 1:id1-1 bu2, ru2 = up((b, r), id2, lp) - Ush[b,2] = U[b,id2,r] - sync_threads() ipl = ipl + 1 - if ru1 == r - gt1 = Ush[bu1,2] + if SFBND && (it == lp.iL[end]) + gt1 = Ubnd[id2] else - if SFBND && (it == lp.iL[end]) - gt1 = Ubnd[id2] - else - gt1 = U[bu1,id2,ru1] - end - end - if ru2 == r - gt2 = Ush[bu2,1] - else - gt2 = U[bu2,id1,ru2] + gt1 = U[bu1,id2,ru1] end if SFBND - S += cG*tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) + S += cG*tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2])) else - S += ztw[ipl]*tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) + S += ztw[ipl]*tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2])) end end end @@ -183,8 +124,6 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd, cG, ztw, r = Int64(CUDA.blockIdx().x) it = point_time((b,r), lp) - Ush = @cuStaticSharedMem(T, (D,2)) - @inbounds begin id1, id2 = lp.plidx[ipl] bu1, ru1 = up((b, r), id1, lp) @@ -192,42 +131,29 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd, cG, ztw, SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1 == N) - Ush[b,1] = U[b,id1,r] - Ush[b,2] = U[b,id2,r] - sync_threads() - - if ru2 == r - gt2 = Ush[bu2,1] + if SFBC && (it == lp.iL[end]) + gt1 = Ubnd[id2] else - gt2 = U[bu2,id1,ru2] - end - if ru1 == r - gt1 = Ush[bu1,2] - else - if SFBC && (it == lp.iL[end]) - gt1 = Ubnd[id2] - else - gt1 = U[bu1,id2,ru1] - end + gt1 = U[bu1,id2,ru1] end - g1 = gt1/gt2 - g2 = Ush[b,2]\Ush[b,1] + g1 = gt1/U[bu2,id1,ru2] + g2 = U[b,id2,r]\U[b,id1,r] if SFBC && (it == 1) - X = cG*projalg(ztw,Ush[b,1]*g1/Ush[b,2]) + X = cG*projalg(ztw,U[b,id1,r]*g1/U[b,id2,r]) frc1[b ,id1, r ] -= X frc2[bu1,id2,ru1] -= cG*projalg(ztw,g1*g2) frc2[bu2,id1,ru2] += cG*projalg(ztw,g2*g1) elseif SFBC && (it == lp.iL[end]) - X = cG*projalg(ztw,Ush[b,1]*g1/Ush[b,2]) + X = cG*projalg(ztw,U[b,id1,r]*g1/U[b,id2,r]) frc1[b ,id1, r ] -= X frc1[b ,id2, r ] += X frc2[bu2,id1,ru2] += cG*projalg(ztw,g2*g1) else - X = projalg(ztw,Ush[b,1]*g1/Ush[b,2]) + X = projalg(ztw,U[b,id1,r]*g1/U[b,id2,r]) frc1[b ,id1, r ] -= X frc1[b ,id2, r ] += X @@ -245,8 +171,6 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, cG, r = Int64(CUDA.blockIdx().x) it = point_time((b, r), lp) - Ush = @cuStaticSharedMem(T, (D,2)) - @inbounds begin id1, id2 = lp.plidx[ipl] bu1, ru1 = up((b, r), id1, lp) @@ -254,156 +178,89 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, cG, SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1 == N) - Ush[b,1] = U[b,id1,r] - Ush[b,2] = U[b,id2,r] - sync_threads() - # H1 staple (b1, r1) = dw((b,r), id2, lp) - if r1 == r - ga = Ush[b1,2] - gb = Ush[b1,1] - else - ga = U[b1,id2,r1] - gb = U[b1,id1,r1] - end - (b2, r2) = up((b1,r1), id1, lp) - if r2 == r - gc = Ush[b2,2] + if SFBC && (it == lp.iL[end]) + gc = Ubnd[id2] else - if SFBC && (it == lp.iL[end]) - gc = Ubnd[id2] - else - gc = U[b2,id2,r2] - end + gc = U[b2,id2,r2] end - h1 = (ga\gb)*gc + h1 = (U[b1,id2,r1]\U[b1,id1,r1])*gc # H2 staple (b1, r1) = up((b,r), id1, lp) - if r1 == r - ga = Ush[b1,1] - else - ga = U[b1,id1,r1] - end + ga = (b2, r2) = up((b1,r1), id1, lp) - if r2 == r - gb = Ush[b2,2] + if SFBC && (it == lp.iL[end]-1) + gb = Ubnd[id2] else - if SFBC && (it == lp.iL[end]-1) - gb = Ubnd[id2] - else - gb = U[b2,id2,r2] - end + 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 + h2 = (U[b1,id1,r1]*gb)/U[b2,id1,r2] # 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] + (b3, r3) = up((b1,r1), id1, lp) + if SFBC && (it == lp.iL[end]) + gc = Ubnd[id2] else - gb = U[b2,id1,r2] + gc = U[b3,id2,r3] end - - (b2, r2) = up((b1,r1), id1, lp) - if r2 == r - gc = Ush[b2,2] - else - if SFBC && (it == lp.iL[end]) - gc = Ubnd[id2] - else - gc = U[b2,id2,r2] - end - end - h3 = (ga*gb)/gc + h3 = (U[b1,id2,r1]*U[b2,id1,r2])/gc # H4 staple (b1, r1) = dw((b,r), id1, lp) - if r1 == r - ga = Ush[b1,1] - gb = Ush[b1,2] - else - ga = U[b1,id1,r1] - gb = U[b1,id2,r1] - end - (b2, r2) = up((b1,r1), id2, lp) - if r2 == r - gc = Ush[b2,1] - else - gc = U[b2,id1,r2] - end - h4 = (ga\gb)*gc + h4 = (U[b1,id1,r1]\U[b1,id2,r1])*U[b2,id1,r2] # END staples - if ru2 == r - gb = Ush[bu2,1] + if SFBC && (it == lp.iL[end]) + ga = Ubnd[id2] else - gb = U[bu2,id1,ru2] - end - if ru1 == r - ga = Ush[bu1,2] - else - if SFBC && (it == lp.iL[end]) - ga = Ubnd[id2] - else - ga = U[bu1,id2,ru1] - end + ga = U[bu1,id2,ru1] end - g1 = ga/gb - g2 = Ush[b,2]\Ush[b,1] + g1 = ga/U[bu2,id1,ru2] + g2 = U[b,id2,r]\U[b,id1,r] if SFBC && (it == 1) - X = (cG*c0)*projalg(Ush[b,1]*g1/Ush[b,2]) + c1*projalg(Ush[b,1]*h2/(Ush[b,2]*gb)) + - (3*c1*cG/2)*projalg(Ush[b,1]*ga/(Ush[b,2]*h3)) + X = (cG*c0)*projalg(U[b,id1,r]*g1/U[b,id2,r]) + c1*projalg(U[b,id1,r]*h2/(U[b,id2,r]*U[bu2,id1,ru2])) + + (3*c1*cG/2)*projalg(U[b,id1,r]*ga/(U[b,id2,r]*h3)) frc1[b,id1,r] -= X frc2[bu1,id2,ru1] -= (cG*c0)*projalg(g1*g2) + (3*c1*cG/2)*projalg((ga/h3)*g2) + - (3*c1*cG/2)*projalg((g1/Ush[b,2])*h1) + (3*c1*cG/2)*projalg((g1/U[b,id2,r])*h1) - frc2[bu2,id1,ru2] += (cG*c0)*projalg(g2*g1) + (3*c1*cG/2) * projalg((Ush[b,2]\h1)*g1) + - c1*projalg(g2*h2/gb) + frc2[bu2,id1,ru2] += (cG*c0)*projalg(g2*g1) + (3*c1*cG/2) * projalg((U[b,id2,r]\h1)*g1) + + c1*projalg(g2*h2/U[bu2,id1,ru2]) elseif SFBC && (it == lp.iL[end]) - X = (cG*c0)*projalg(Ush[b,1]*g1/Ush[b,2]) + - (3*c1*cG/2) * (projalg(Ush[b,1]*ga/(Ush[b,2]*h3))) + X = (cG*c0)*projalg(U[b,id1,r]*g1/U[b,id2,r]) + + (3*c1*cG/2) * (projalg(U[b,id1,r]*ga/(U[b,id2,r]*h3))) - frc1[b,id1,r] -= X + c1*projalg(Ush[b,1]*g1/h4) - frc1[b,id2,r] += X + (3*c1*cG/2)*projalg(h1*g1/Ush[b,2]) + frc1[b,id1,r] -= X + c1*projalg(U[b,id1,r]*g1/h4) + frc1[b,id2,r] += X + (3*c1*cG/2)*projalg(h1*g1/U[b,id2,r]) - frc2[bu2,id1,ru2] += (cG*c0)*projalg(g2*g1) + (3*c1*cG/2) * projalg((Ush[b,2]\h1)*g1) + - c1 * projalg(h4\Ush[b,1]*g1) + frc2[bu2,id1,ru2] += (cG*c0)*projalg(g2*g1) + (3*c1*cG/2) * projalg((U[b,id2,r]\h1)*g1) + + c1 * projalg(h4\U[b,id1,r]*g1) else zsq = ztw[ipl]^2 - X = projalg(c0*ztw[ipl],Ush[b,1]*g1/Ush[b,2]) + projalg(zsq*c1,Ush[b,1]*h2/(Ush[b,2]*gb)) + - projalg(zsq*c1,Ush[b,1]*ga/(Ush[b,2]*h3)) + X = projalg(c0*ztw[ipl],U[b,id1,r]*g1/U[b,id2,r]) + projalg(zsq*c1,U[b,id1,r]*h2/(U[b,id2,r]*U[bu2,id1,ru2])) + + projalg(zsq*c1,U[b,id1,r]*ga/(U[b,id2,r]*h3)) - frc1[b,id1,r] -= X + projalg(zsq*c1,Ush[b,1]*g1/h4) - frc1[b,id2,r] += X + projalg(zsq*c1,h1*g1/Ush[b,2]) + frc1[b,id1,r] -= X + projalg(zsq*c1,U[b,id1,r]*g1/h4) + frc1[b,id2,r] += X + projalg(zsq*c1,h1*g1/U[b,id2,r]) frc2[bu1,id2,ru1] -= projalg(c0*ztw[ipl],g1*g2) + projalg(zsq*c1,(ga/h3)*g2) + - projalg(zsq*c1,(g1/h4)*Ush[b,1]) + projalg(zsq*c1,(g1/Ush[b,2])*h1) + projalg(zsq*c1,(g1/h4)*U[b,id1,r]) + projalg(zsq*c1,(g1/U[b,id2,r])*h1) - frc2[bu2,id1,ru2] += projalg(c0*ztw[ipl],g2*g1) + projalg(zsq*c1,(Ush[b,2]\h1)*g1) + - projalg(zsq*c1,g2*h2/gb) + projalg(zsq*c1,h4\Ush[b,1]*g1) + frc2[bu2,id1,ru2] += projalg(c0*ztw[ipl],g2*g1) + projalg(zsq*c1,(U[b,id2,r]\h1)*g1) + + projalg(zsq*c1,g2*h2/U[bu2,id1,ru2]) + projalg(zsq*c1,h4\U[b,id1,r]*g1) end end diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl index 3d5aa40..4a0b374 100644 --- a/src/YM/YMflow.jl +++ b/src/YM/YMflow.jl @@ -96,41 +96,25 @@ function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::S SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) - Ush = @cuStaticSharedMem(TG, D) - Fsh = @cuStaticSharedMem(TA, D) - @inbounds for id in 1:N - Ush[b] = U[b,id,r] - Fsh[b] = frc[b,id,r] - sync_threads() - bu, ru = up((b,r), id, lp) bd, rd = dw((b,r), id, lp) - if ru == r - X = Fsh[bu] - else - X = frc[bu,id,ru] - end - if rd == r - Y = Fsh[bd] - Ud = Ush[bd] - else - Y = frc[bd,id,rd] - Ud = U[bd,id,rd] - end + X = frc[bu,id,ru] + Y = frc[bd,id,rd] + Ud = U[bd,id,rd] if SFBC if (it > 1) && (it < lp.iL[end]) - frc2[b,id,r] = (5/6)*Fsh[b] + (1/6)*(projalg(Ud\Y*Ud) + - projalg(Ush[b]*X/Ush[b])) + frc2[b,id,r] = (5/6)*frc[b,id,r] + (1/6)*(projalg(Ud\Y*Ud) + + projalg(U[b,id,r]*X/U[b,id,r])) elseif (it == lp.iL[end]) && (id < N) - frc2[b,id,r] = (5/6)*Fsh[b] + (1/6)*(projalg(Ud\Y*Ud) + - projalg(Ush[b]*X/Ush[b])) + frc2[b,id,r] = (5/6)*frc[b,id,r] + (1/6)*(projalg(Ud\Y*Ud) + + projalg(U[b,id,r]*X/U[b,id,r])) end else - frc2[b,id,r] = (5/6)*Fsh[b] + (1/6)*(projalg(Ud\Y*Ud) + - projalg(Ush[b]*X/Ush[b])) + frc2[b,id,r] = (5/6)*frc[b,id,r] + (1/6)*(projalg(Ud\Y*Ud) + + projalg(U[b,id,r]*X/U[b,id,r])) end end end @@ -428,87 +412,59 @@ function krnl_field_tensor!(frc1::AbstractArray{TA}, frc2, U::AbstractArray{T}, it = point_time((b,r), lp) SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) - Ush = @cuStaticSharedMem(T, (D,2)) - #First plane id1, id2 = lp.plidx[ipl1] - Ush[b,1] = U[b,id1,r] - Ush[b,2] = U[b,id2,r] - sync_threads() - SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1 == 4) bu1, ru1 = up((b, r), id1, lp) bu2, ru2 = up((b, r), id2, lp) bd, rd = up((bu1, ru1), id2, lp) - if ru1 == r - gt1 = Ush[bu1,2] + if SFBC && (it == lp.iL[end]) + gt1 = Ubnd[id2] else - if SFBC && (it == lp.iL[end]) - gt1 = Ubnd[id2] - else - gt1 = U[bu1,id2,ru1] - end - end - if ru2 == r - gt2 = Ush[bu2,1] - else - gt2 = U[bu2,id1,ru2] + gt1 = U[bu1,id2,ru1] end - l1 = gt1/gt2 - l2 = Ush[b,2]\Ush[b,1] + l1 = gt1/U[bu2,id1,ru2] + l2 = U[b,id2,r]\U[b,id1,r] if SFBC && (it == lp.iL[end]) - frc1[b,1,r] = projalg(Ush[b,1]*l1/Ush[b,2]) + frc1[b,1,r] = projalg(U[b,id1,r]*l1/U[b,id2,r]) frc1[bu1,2,ru1] = zero(TA) frc1[bd,3,rd] = zero(TA) frc1[bu2,4,ru2] = projalg(l2*l1) else - frc1[b,1,r] = projalg(ztw1, Ush[b,1]*l1/Ush[b,2]) + frc1[b,1,r] = projalg(ztw1, U[b,id1,r]*l1/U[b,id2,r]) frc1[bu1,2,ru1] = projalg(ztw1, l1*l2) - frc1[bd,3,rd] = projalg(ztw1, gt2\(l2*gt1)) + frc1[bd,3,rd] = projalg(ztw1, U[bu2,id1,ru2]\(l2*gt1)) frc1[bu2,4,ru2] = projalg(ztw1, l2*l1) end # Second plane id1, id2 = lp.plidx[ipl2] - Ush[b,1] = U[b,id1,r] - Ush[b,2] = U[b,id2,r] - sync_threads() - SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1 == 4) bu1, ru1 = up((b, r), id1, lp) bu2, ru2 = up((b, r), id2, lp) bd, rd = up((bu1, ru1), id2, lp) - if ru1 == r - gt1 = Ush[bu1,2] + if SFBC && (it == lp.iL[end]) + gt1 = Ubnd[id2] else - if SFBC && (it == lp.iL[end]) - gt1 = Ubnd[id2] - else - gt1 = U[bu1,id2,ru1] - end - end - if ru2 == r - gt2 = Ush[bu2,1] - else - gt2 = U[bu2,id1,ru2] + gt1 = U[bu1,id2,ru1] end - l1 = gt1/gt2 - l2 = Ush[b,2]\Ush[b,1] + l1 = gt1/U[bu2,id1,ru2] + l2 = U[b,id2,r]\U[b,id1,r] if SFBC && (it == lp.iL[end]) - frc2[b,1,r] = projalg(Ush[b,1]*l1/Ush[b,2]) + frc2[b,1,r] = projalg(U[b,id1,r]*l1/U[b,id2,r]) frc2[bu1,2,ru1] = zero(TA) frc2[bd,3,rd] = zero(TA) frc2[bu2,4,ru2] = projalg(l2*l1) else - frc2[b,1,r] = projalg(ztw2, Ush[b,1]*l1/Ush[b,2]) + frc2[b,1,r] = projalg(ztw2, U[b,id1,r]*l1/U[b,id2,r]) frc2[bu1,2,ru1] = projalg(ztw2, l1*l2) - frc2[bd,3,rd] = projalg(ztw2, gt2\(l2*gt1)) + frc2[bd,3,rd] = projalg(ztw2, U[bu2,id1,ru2]\(l2*gt1)) frc2[bu2,4,ru2] = projalg(ztw2, l2*l1) end end