[BUG]: cuda.compute: `reduce_into` / `segmented_reduce` with `d_out` dtype != `h_init` dtype silently write reinterpreted accumulator bytes
- Dominant language
- C++
- Stars
- 2.5k
- Forks
- 487
- Avg merge
- 2d 7h
- Merged PRs (30d)
- 296
Description
### Is this a duplicate?
- [x] I confirmed there appear to be no [duplicate issues](https://github.com/NVIDIA/cccl/issues) for this bug and that I agree to the [Code of Conduct](CODE_OF_CONDUCT.md)
### Type of Bug
Silent Failure
### Component
cuda.compute (Python)
### Describe the bug
`make_reduce_into` / `reduce_into` accepts a `d_out` array whose dtype differs from the `h_init` dtype without any validation. The generated kernel treats the output pointer as a pointer to the accumulator type and performs a raw (bit-reinterpreting) store — no conversion, no error. The call succeeds and the output contains garbage. `segmented_reduce` has the identical bug.
Root cause: `c/parallel/src/reduce.cu` (~line 213) and `c/parallel/src/segmented_reduce.cu` (~line 149) specialize the output iterator with `accum_t` (which is `init.type`, per `get_accumulator_type`), so the final store is effectively `*(accum_t*)d_out_ptr = accumulator`. Nothing in the Python or C layer compares `d_out`'s type to `accum_t`. These are the only two users of the shared `output_iterator_traits` JIT template that pass `accum_t`; `transform.cu` and `three_way_partition.cu` already pass `output_it.value_type` (the correct, converting pattern), and `scan.cu` types pointer outputs with the output's own value type — so scan converts correctly and reduce/segmented_reduce are the outliers.
For output dtypes SMALLER than the accumulator (e.g. float16 out, float32 accumulator) this is also an out-of-bounds write: 4 bytes are stored into a 2-byte allocation (masked in practice by allocator padding; compute-sanitizer should flag it). For `segmented_reduce` it is worse: the kernel strides the output by `sizeof(accum_t)` instead of the element size, so with N segments the stores corrupt *neighboring elements* and run past the end of the allocation (see segmented output table below).
### How to Reproduce
`reduce_into` reproducer:
```python
import numpy as np, cupy as cp
import cuda.compute
d_in = cp.ones(3000, dtype=cp.float16)
h_init = np.zeros((), dtype=np.float32) # accumulator := float32
def run(d_out):
op = cuda.compute.OpKind.PLUS
red = cuda.compute.make_reduce_into(d_in=d_in, d_out=d_out, op=op, h_init=h_init)
n = red(temp_storage=None, d_in=d_in, d_out=d_out,
num_items=d_in.size, op=op, h_init=h_init)
tmp = cp.empty(n, dtype=cp.uint8)
red(temp_storage=tmp, d_in=d_in, d_out=d_out,
num_items=d_in.size, op=op, h_init=h_init)
return d_out[0]
print("float32 out:", run(cp.zeros(1, dtype=cp.float32))) # control
print("float16 out:", run(cp.zeros(1, dtype=cp.float16)))
print("int32 out:", run(cp.zeros(1, dtype=cp.int32)))
print("float64 out:", run(cp.zeros(1, dtype=cp.float64)))
```
`segmented_reduce` reproducer:
```python
import numpy as np, cupy as cp
import cuda.compute
n_segments, segment_size = 4, 750
offsets = cp.arange(n_segments + 1, dtype="int64") * segment_size
d_in = cp.ones(n_segments * segment_size, dtype=cp.float16)
starts, ends = offsets[:-1], offsets[1:]
h_init = np.zeros((), dtype=np.float32) # accumulator := float32
def run(out_dtype):
d_out = cp.zeros(n_segments, dtype=out_dtype)
cuda.compute.segmented_reduce(
d_in=d_in, d_out=d_out, num_segments=n_segments,
start_offsets_in=starts, end_offsets_in=ends,
op=cuda.compute.OpKind.PLUS, h_init=h_init)
return d_out.get()
```
### Expected behavior
For `reduce_into`: the correct sum is 3000.0 (f32 accumulation of 3000 fp16 ones), bit pattern `0x453B8000`.
| case | actual output | expected | what actually happened |
|---|---|---|---|
| float32 out (control) | `3000.0` | `3000.0` | correct |
| float16 out | `-0.0` | `3000.0` (or an error) | low 2 bytes of `0x453B8000` = `0x8000` = fp16 `-0.0`; 4-byte store into a 2-byte buffer (OOB) |
| int32 out | `1161527296` | `3000` (or an error) | `1161527296 == 0x453B8000` — the float bit pattern verbatim |
| float64 out | `5.73870734e-315` | `3000.0` (or an error) | 4 bytes of float bits in the low half of the double, upper half stale → denormal |
For `segmented_reduce`: each segment sums 750 fp16 ones in f32; the correct per-segment result is 750.0 (f32 bit pattern `0x443B8000`, LE bytes `00 80 3B 44`).
| case | actual output | expected | what actually happened |
|---|---|---|---|
| float32 out (control) | `[750. 750. 750. 750.]` | same | correct |
| float16 out | `[-0. 4.23 -0. 4.23]` | `[750. 750. 750. 750.]` (or an error) | each 4-byte store covers TWO fp16 slots (`0x8000` = `-0.0`, `0x443B` = `4.23`); segments 0-1 fill the whole 8-byte buffer, segments 2-3 write 8 bytes past the allocation (OOB) |
| int32 out | `[1144750080 ×4]` | `[750 ×4]` (or an error) | `1144750080 == 0x443B8000` — the f32 bit pattern verbatim (element size matches, so no striding damage) |
| float64 out | `[5.07e+20 5.07e+20 0. 0.]` | `[750. 750. 750. 750.]` (or an error) | 4-byte stores at f32 stride pack two segments' bit patterns into each 8-byte element; elements 2-3 never written (stale) |
No exception or warning is raised in any case.
### Reproduction link
_No response_
### Operating System
_No response_
### nvidia-smi output
_No response_
### NVCC version
_No response_
Contributor guide
Assessment
This issue has not been assessed yet.