Kernel only gives proper result when adding @print()

[SouthEndMusic]

I'm experiencing a weird bug where the kernel I wrote only gives the expected result when I add a print statement in a particular place:

https://github.com/SouthEndMusic/SplineGrids.jl/blob/08f0e91a70c547e807b0a8bd89058b4cb8c06a43/src/spline_dimension.jl#L48-L70

It was suggested on Slack that it could be a synchronization issue, but the different threads are completely independent from eachother.

[vchuravy]

What backend are you executing this one? Can you isolate this into a MWE.

[SouthEndMusic]

MWE:

using KernelAbstractions
using Plots

@kernel function spline_dimension_kernel!(
    eval, @Const(knots_all), @Const(sample_points), @Const(sample_indices), degree)
    l = @index(Global, Linear)
    t = sample_points[l]
    i = sample_indices[l]

    eval[l, 1] = one(eltype(eval))
    for k in 1:degree
        @print()
        B_old = eval[l, 1]
        eval[l, 1] = zero(eltype(eval))
        for k_ in 1:k
            t_min = knots_all[i+k_-k]
            t_max = knots_all[i+k_]
            Δt = t_max - t_min
            frac = B_old / Δt
            B_old = eval[l, k_+1] # Value for next iteration
            # Additions sum to B_old => partition of unity
            eval[l, k_] += frac * (t_max - t)
            eval[l, k_+1] = frac * (t - t_min)
        end
    end
end

function get_index(knot_vector::KnotVector, t, d)
    clamp(searchsortedlast(knot_vector.knots_all, t),
        d + 1, length(knot_vector.knots_all) - d - 1)
end

n_samples = 100
degree = 2
knots_all = [0.0, 0.0, 0.0, 0.25, 0.5, 0.75, 1.0, 1.0, 1.0]
sample_points = collect(range(0.0, 1.0, length=n_samples))
sample_indices = zeros(Int, 100)
sample_indices[1:25] .= 3
sample_indices[26:50] .= 4
sample_indices[51:75] .= 5
sample_indices[76:100] .= 6
eval_ = zeros(n_samples, degree + 1)

backend = get_backend(eval_)
spline_dimension_kernel!(backend)(
    eval_, knots_all, sample_points, sample_indices,
    degree, ndrange=n_samples)
synchronize(backend)

n_basis_functions = length(knots_all) - degree - 1
out = zeros(n_samples, n_basis_functions)

for (l, i) in enumerate(sample_indices)
    out[l, (i-degree):i] .= eval_[l, :]
end

p = plot()

for j in 1:n_basis_functions
    plot!(sample_points, out[:, j])
end

p

Result with print statement:

Result without print statement: