In warp_reduce_kernel, tail_col starts at (row_size // VEC_STRIDE) * VEC_STRIDE + lid, but the vector loop already covers all elements < vec_limit (including the remainder when row_size < VEC_STRIDE or simd_width == 1). This causes some elements to be reduced twice (e.g. row_size=64, simd_width=8; or row_size=65, simd_width=1), producing incorrect results. Compute the scalar-tail start from the number of full SIMD vectors instead (e.g. tail_start = (row_size // simd_width) * simd_width, then tail_col = tail_start + lid) so only the final row_size % simd_width elements are handled in the tail loop.

The per-thread accumulation in warp_reduce_kernel uses reduce_fn(accum, v) (accumulator as the first operand). Elsewhere in the GPU backend (and the CPU backend) the reduction wrapper is consistently applied as reduce_fn(val, acc) (new value first), so this flips operand order for non-commutative reductions and makes results kernel-dependent. Swap the argument order here (and similarly in the scalar tail / lane-collapse paths) to match the established reduce_fn(val, acc) convention.

row_reduce now manually collapses the SIMD accumulator to a scalar before block_reduce, but the lane fold uses reduce_fn(lane_accum, accum[i][lane]) (accumulator first). Other call sites in this file reduce as reduce_fn(val, acc) (value first), and the std reduction wrapper contract uses (val, acc). To keep semantics consistent across kernels (especially for non-commutative reductions), flip the operand order in this lane-collapse fold.

This doesn't seem right to me. The host side results array should be of out_size. Shouldn't instead res_buffer be set to size out_size?

Don't need to recompute is_contig here. Just use reduce_contig_dim and pull reduce_simd out of the for loop.