FEM first-sample interpolation returns zero gradients with triplet construction

import warp as wp
import warp.fem as fem
import warp.sparse as sparse


@fem.integrand
def scaled_field(s: fem.Sample, u: fem.Field, scale: wp.array[float]):
    return scale[0] * u(s)


@wp.kernel
def sum_values(values: wp.array[float], loss: wp.array[float]):
    wp.atomic_add(loss, 0, values[wp.tid()])


with wp.ScopedDevice("cpu"):
    space = fem.make_polynomial_space(fem.Grid2D(res=wp.vec2i(1)), degree=1)
    restriction = fem.make_space_restriction(space_topology=space.topology)
    matrix = sparse.bsr_zeros(space.node_count(), space.node_count(), float)
    matrix.values.requires_grad = True
    scale = wp.array([2.0], dtype=float, requires_grad=True)
    loss = wp.zeros(1, dtype=float, requires_grad=True)
    with wp.Tape() as tape:
        fem.interpolate(
            scaled_field, dest=matrix, dest_space=space, at=restriction,
            fields={"u": fem.make_trial(space)}, values={"scale": scale},
            reduction="first", bsr_options={"construction": "triplets"},
            kernel_options={"enable_backward": True},
        )
        wp.launch(sum_values, dim=matrix.nnz_sync(), inputs=[matrix.values, loss])
    tape.backward(loss=loss)
    print(f"loss={loss.numpy()[0]}, gradient={scale.grad.numpy()[0]}")
    # Interpolation at these four Q1 nodes is scale * I, so d(sum(matrix))/d(scale) = 4.
    assert loss.numpy()[0] == 8.0
    assert scale.grad.numpy()[0] == 4.0