[maleadt]

Follow-up to JuliaGPU/CUDA.jl#3446. Base reductions of GPU arrays go through AK.mapreducedim!, and two common layouts of strided reductions into contiguous outputs ran far below bandwidth:

• More outputs than elements per output, e.g. sum(x; dims=2) of a 3072×2048 matrix, took one thread per output. dst_size threads leave the GPU mostly idle below tens of thousands of outputs: 90 µs for that BFloat16 matrix on an RTX 5080, 140 GB/s. Only exactly square shapes took the tiled kernel.
• The middle dimension of a 3-d array, e.g. sum(x; dims=2) of a (features, sequence, batch) array, fell back to one thread or one block per output, since the columns kernel only took contiguous outputs: 46–300 GB/s.

This PR:

1. Picks the kernel for these reductions in one place, _columns_shape:
   • one thread per output only when a row of outputs is at least 64 KiB, or the reduction is shorter than 32;
   • else the tiled kernel (8 outputs per block, one pass) when there are at least as many outputs as elements per output and 32-output tiles would not fill the device;
   • else the columns kernel (32 outputs per block, split across blocks as before).
2. Lets the columns kernel reduce runs of contiguous outputs at the offsets of the other kept dimensions, so the middle dimensions of arrays take it too.

Both are in the spirit of PyTorch's Reduce.cuh: a fixed 2-d block shape for strided reductions, split across blocks only when the outputs alone don't fill the device. No new tuning parameters.

Timings

sum(x; dims=2), base → this PR (RTX 5080 and RTX 6000 Ada: GPU time; Arc A770 and M1: AK.reduce per call, back to back):

	RTX 5080	RTX 6000 Ada	Arc A770	M1
768×1024×8, BFloat16 (Float16 on A770, M1)	42.0 → 19.6 µs	275 → 21.5 µs	420 → 69 µs	338 → 257 µs
768×1024×8, Float32	43.6 → 33.3 µs	141 → 35.8 µs	420 → 113 µs	430 → 456 µs
64×4096×32, Float32	138 → 41 µs		383 → 120 µs	910 → 575 µs
3072×2048, BFloat16 (Float16 on A770, M1)	90 → 16.9 µs		827 → 90 µs	589 → 266 µs
8192×8192, Float32	416 → 295 µs

Geometric mean of per-call time over 26 shapes (Float32 and 16-bit, base = 1): RTX 5080 0.58, RTX 6000 Ada 0.60, Arc A770 0.47, M1 0.80, Ryzen 9950X iGPU 0.71, PoCL 0.67. The other shapes in CUDA.jl#3446 are unchanged.

Known trade-off: on devices with a few compute units (a Ryzen 9950X's iGPU, Jetson Xavier), 2048×768 takes 1.4× as long. One thread per output already fills them with 8 blocks, and making that threshold device-relative didn't seem worth the extra tuning.

Tried and left out

• Narrower unsplit column tiles (8/16 outputs), which avoid the second pass for few outputs: up to 1.5× faster on NVIDIA, but up to 5× slower on AMD and the Jetsons, and four more thresholds.
• Fusing the second pass into the first with a last-block counter, like PyTorch's global_reduce: zeroing the counters costs as much as the pass it saves.

Tested with the reduce tests on CUDA, oneAPI, Metal, AMDGPU and OpenCL (PoCL), plus the full test suite on CUDA.