dim-reduce: implement the 7 dimension-reduction kernels + kernel backlog

KERNEL_BACKLOG.md

@@ -0,0 +1,95 @@
+# Kernel backlog — remaining CUDA kernels to implement
+
+Ranked easiest→hardest by tensara difficulty (from `Problems.mhtml`: EASY/MEDIUM/
+HARD), then clustered by similarity to each other and to already-implemented
+kernels/frameworks. Reuse the existing frameworks — most clusters are "one
+framework, many kernels."
+
+**Existing frameworks to reuse**
+- Dimension reductions: `kernel-implementation/dim-reduce.cuh`
+- Global scalar reductions (losses): `block-reduce.cuh` + `loss.cuh` (`SmemTreeReduce`)
+- Row-wise norms: `reduction.cuh`
+- Activations / scalar ops: `activation.cuh`, `scalar-ops.cuh`
+- Matmul + epilogue fusion: `gemm-epilogue.cuh`
+- Matmul + prologue (dequant): `gemm-prologue.cuh`
+- Pooling: `pooling.cuh`
+
+> Note: `cosine-similarity`, `frobenius-norm`, `triplet-margin` are implemented in
+> PR #38 (unmerged) — stubs on `main` but effectively done.
+
+---
+
+## EASY
+
+### ✅ Dimension reductions — DONE (`dim-reduce.cuh`)
+`sum-dim`, `mean-dim`, `max-dim`, `min-dim`, `product-dim`, `argmax`, `argmin`
+
+### Global-reduction norms/stats — (frobenius/cosine in #38)
+`frobenius-norm`, `cosine-similarity` — reuse `block-reduce.cuh` / `reduction.cuh`.
+
+### Elementwise / image stencils — like `grayscale`, `threshold`, `conv-1d`
+`matrix-scalar` (elementwise ×s), `box-blur` (2-D box stencil), `edge-detect`
+(Sobel 2-D stencil), `histogram` (atomic scatter into bins)
+
+### Simple / misc
+`diagonal-matmul` (row-scale ≈ elementwise), `running-sum-1d` (1-D window ≈ `conv-1d`),
+`ecc-point-negation` (batched modular negate ≈ elementwise)
+
+### Low-precision dequant — prologue of `int8-weight-gemm` (`gemm-prologue.cuh`)
+`mxfp4-dequantize`, `mxfp8-dequantize`
+
+---
+
+## MEDIUM
+
+### Core matmul — `gemm-epilogue.cuh`, CS4803 lab2
+`matrix-multiplication`, `square-matmul`, `symmetric-matmul`,
+`upper-trig-matmul`, `lower-trig-matmul` (masked GEMM), `matrix-power` (repeated GEMM)
+
+### Fused GEMM/conv + epilogue — directly `gemm-epilogue.cuh` + `scalar-ops.cuh`
+`matmul-swish`, `matmul-sigmoid-sum`, `gemm-multiply-leakyrelu`, `conv2d-relu-hardswish`
+
+### 2-D conv & pooling — dimensionality-up from 1-D (`pooling.cuh`, CS4803 lab3)
+`conv-2d`, `avg-pool-2d`, `max-pool-2d`
+
+### Per-row normalize with mean+var — extends `rms-norm`/`mean-subtract`
+`layer-norm`, `batch-norm`
+
+### Scan — new pattern (`running-sum-1d` is the warm-up)
+`cumsum`, `cumprod`
+
+### Low-precision quantize — inverse of dequant
+`mxfp4-quantize`, `mxfp8-quantize`, `nvfp4-quantize`, `nvfp4-dequantize`
+
+### Graphs — new pattern (iterative relaxation)
+`all-pairs-shortest-path` (Floyd–Warshall), `shortest-path` (Bellman-Ford/Dijkstra),
+`min-spanning-tree` (Prim)
+
+### Finite-field / crypto — new pattern (needs modulus spec; currently test-uncovered)
+`poly-multiply-ff`, `vector-multiply-ff`
+
+### (`triplet-margin` — in #38, two-level reduce)
+
+---
+
+## HARD
+
+### 3-D conv/pool — one more dimension than the 2-D cluster
+`conv-square-3d`, `avg-pool-3d`, `max-pool-3d`
+
+### Batched/high-dim matmul
+`matmul-3d`, `matmul-4d`
+
+### Low-precision GEMM — prologue-dequant matmul (`gemm-prologue.cuh`)
+`mxfp4-gemm`, `mxfp8-gemm`, `nvfp4-gemm`, `nvfp4-gemv`
+
+### Attention — matmul→softmax→matmul (uses `softmax` + GEMM; FlashAttention territory)
+`scaled-dot-attention`
+
+---
+
+## Recommended order (max framework reuse first)
+1. ✅ dim reductions → 2. norm/stats (merge #38; then layer-norm/batch-norm) →
+3. fused GEMM epilogues → 4. core matmul → 5. 2-D conv/pool →
+6. dequant → quantize → fp GEMM (one lineage) → 7. scan →
+8. 3-D conv/pool + 3-D/4-D matmul → 9. graphs & finite-field → 10. attention.

kernel-implementation/dim-reduce.cuh

@@ -0,0 +1,74 @@
+#pragma once
+// Generic reduction along one axis `dim` of an N-D tensor described by shape/ndim.
+// A tensor factors as [outer, L, inner] around the reduced axis (L = shape[dim]),
+// so output[o, i] = reduce over l of input[o*L*inner + l*inner + i]. One thread
+// per output element. `shape` is a device pointer, copied to host to compute the
+// three extents. Covers sum/mean/max/min/product (value) and argmax/argmin (index).
+#include <cuda_runtime.h>
+#include <cfloat>
+
+// ---- value ops: init / combine / finalize(acc, L) --------------------------
+struct SumOp  { __device__ static float init(){return 0.0f;}      __device__ static float combine(float a,float b){return a+b;}        __device__ static float finalize(float a,long){return a;} };
+struct MeanOp { __device__ static float init(){return 0.0f;}      __device__ static float combine(float a,float b){return a+b;}        __device__ static float finalize(float a,long L){return a/(float)L;} };
+struct MaxOp  { __device__ static float init(){return -FLT_MAX;}  __device__ static float combine(float a,float b){return fmaxf(a,b);}  __device__ static float finalize(float a,long){return a;} };
+struct MinOp  { __device__ static float init(){return  FLT_MAX;}  __device__ static float combine(float a,float b){return fminf(a,b);}  __device__ static float finalize(float a,long){return a;} };
+struct ProdOp { __device__ static float init(){return 1.0f;}      __device__ static float combine(float a,float b){return a*b;}        __device__ static float finalize(float a,long){return a;} };
+
+template <class Op>
+__global__ void dimreduce_kernel(const float* __restrict__ in, float* __restrict__ out,
+                                 long outer, long L, long inner) {
+    long idx = (long)blockIdx.x * blockDim.x + threadIdx.x;   // output element
+    if (idx >= outer * inner) return;
+    long o = idx / inner, i = idx % inner;
+    const float* base = in + o * L * inner + i;
+    float acc = Op::init();
+    for (long l = 0; l < L; ++l) acc = Op::combine(acc, base[l * inner]);
+    out[idx] = Op::finalize(acc, L);
+}
+
+// ---- argmax/argmin: SIGN=+1 -> max (v>best), SIGN=-1 -> min (v<best) --------
+template <int SIGN>
+__global__ void argreduce_kernel(const float* __restrict__ in, int* __restrict__ out,
+                                 long outer, long L, long inner) {
+    long idx = (long)blockIdx.x * blockDim.x + threadIdx.x;
+    if (idx >= outer * inner) return;
+    long o = idx / inner, i = idx % inner;
+    const float* base = in + o * L * inner + i;
+    float best = base[0]; int bi = 0;
+    for (long l = 1; l < L; ++l) {
+        float v = base[l * inner];
+        if (SIGN * v > SIGN * best) { best = v; bi = (int)l; }   // strict => first index wins
+    }
+    out[idx] = bi;
+}
+
+// ---- extent helper: outer/L/inner from a device shape (T = size_t or int) --
+template <class T>
+static inline void dimreduce_extents(const T* shape_dev, int ndim, int dim,
+                                     long& outer, long& L, long& inner) {
+    T hshape[16];
+    cudaMemcpy(hshape, shape_dev, (size_t)ndim * sizeof(T), cudaMemcpyDeviceToHost);
+    outer = 1; inner = 1; L = (long)hshape[dim];
+    for (int a = 0; a < dim; ++a)       outer *= (long)hshape[a];
+    for (int a = dim + 1; a < ndim; ++a) inner *= (long)hshape[a];
+}
+
+template <class Op>
+static inline void launch_dimreduce(const float* in, float* out,
+                                    const size_t* shape_dev, size_t ndim, int dim) {
+    long outer, L, inner;
+    dimreduce_extents<size_t>(shape_dev, (int)ndim, dim, outer, L, inner);
+    long nout = outer * inner;
+    int BS = 256; long grid = (nout + BS - 1) / BS;
+    dimreduce_kernel<Op><<<grid, BS>>>(in, out, outer, L, inner);
+}
+
+template <int SIGN>
+static inline void launch_argreduce(const float* in, int* out,
+                                    const int* shape_dev, int ndim, int dim) {
+    long outer, L, inner;
+    dimreduce_extents<int>(shape_dev, ndim, dim, outer, L, inner);
+    long nout = outer * inner;
+    int BS = 256; long grid = (nout + BS - 1) / BS;
+    argreduce_kernel<SIGN><<<grid, BS>>>(in, out, outer, L, inner);
+}

solutions-cuda/argmax.cu

@@ -1,16 +1,6 @@
-// Solution stub for "argmax".
-// TODO: implement the body below. The signature is derived from
-//       kernel-harnesses/argmax.cu and must stay in sync with it.
-//
-// Build it (the build system auto-picks this file):
-//   make build/bin/argmax.exe
-//   ./build/bin/argmax.exe
-#include <cuda_runtime.h>
-#include <cuda_fp16.h>
-#include <cuda_fp8.h>
-#include <cstdint>
+// argmax — index of max along dim of an N-D tensor. Uses the shared dim-reduce framework.
+#include "../kernel-implementation/dim-reduce.cuh"
 
-// Note: all pointer arguments are device pointers.
 extern "C" void solution(const float* input, int dim, int* output, const int* shape, int ndim) {
-    // TODO: implement argmax
+    launch_argreduce<1>(input, output, shape, ndim, dim);
 }

solutions-cuda/argmin.cu

@@ -1,16 +1,6 @@
-// Solution stub for "argmin".
-// TODO: implement the body below. The signature is derived from
-//       kernel-harnesses/argmin.cu and must stay in sync with it.
-//
-// Build it (the build system auto-picks this file):
-//   make build/bin/argmin.exe
-//   ./build/bin/argmin.exe
-#include <cuda_runtime.h>
-#include <cuda_fp16.h>
-#include <cuda_fp8.h>
-#include <cstdint>
+// argmin — index of min along dim of an N-D tensor. Uses the shared dim-reduce framework.
+#include "../kernel-implementation/dim-reduce.cuh"
 
-// Note: all pointer arguments are device pointers.
 extern "C" void solution(const float* input, int dim, int* output, const int* shape, int ndim) {
-    // TODO: implement argmin
+    launch_argreduce<-1>(input, output, shape, ndim, dim);
 }

solutions-cuda/max-dim.cu

@@ -1,16 +1,6 @@
-// Solution stub for "max-dim".
-// TODO: implement the body below. The signature is derived from
-//       kernel-harnesses/max-dim.cu and must stay in sync with it.
-//
-// Build it (the build system auto-picks this file):
-//   make build/bin/max-dim.exe
-//   ./build/bin/max-dim.exe
-#include <cuda_runtime.h>
-#include <cuda_fp16.h>
-#include <cuda_fp8.h>
-#include <cstdint>
+// max-dim — reduce along dim of an N-D tensor (shape/ndim). Uses the shared dim-reduce framework.
+#include "../kernel-implementation/dim-reduce.cuh"
 
-// Note: all pointer arguments are device pointers.
 extern "C" void solution(const float* input, int dim, float* output, const size_t* shape, size_t ndim) {
-    // TODO: implement max-dim
+    launch_dimreduce<MaxOp>(input, output, shape, ndim, dim);
 }

solutions-cuda/mean-dim.cu

@@ -1,16 +1,6 @@
-// Solution stub for "mean-dim".
-// TODO: implement the body below. The signature is derived from
-//       kernel-harnesses/mean-dim.cu and must stay in sync with it.
-//
-// Build it (the build system auto-picks this file):
-//   make build/bin/mean-dim.exe
-//   ./build/bin/mean-dim.exe
-#include <cuda_runtime.h>
-#include <cuda_fp16.h>
-#include <cuda_fp8.h>
-#include <cstdint>
+// mean-dim — reduce along dim of an N-D tensor (shape/ndim). Uses the shared dim-reduce framework.
+#include "../kernel-implementation/dim-reduce.cuh"
 
-// Note: all pointer arguments are device pointers.
 extern "C" void solution(const float* input, int dim, float* output, const size_t* shape, size_t ndim) {
-    // TODO: implement mean-dim
+    launch_dimreduce<MeanOp>(input, output, shape, ndim, dim);
 }

solutions-cuda/min-dim.cu

@@ -1,16 +1,6 @@
-// Solution stub for "min-dim".
-// TODO: implement the body below. The signature is derived from
-//       kernel-harnesses/min-dim.cu and must stay in sync with it.
-//
-// Build it (the build system auto-picks this file):
-//   make build/bin/min-dim.exe
-//   ./build/bin/min-dim.exe
-#include <cuda_runtime.h>
-#include <cuda_fp16.h>
-#include <cuda_fp8.h>
-#include <cstdint>
+// min-dim — reduce along dim of an N-D tensor (shape/ndim). Uses the shared dim-reduce framework.
+#include "../kernel-implementation/dim-reduce.cuh"
 
-// Note: all pointer arguments are device pointers.
 extern "C" void solution(const float* input, int dim, float* output, const size_t* shape, size_t ndim) {
-    // TODO: implement min-dim
+    launch_dimreduce<MinOp>(input, output, shape, ndim, dim);
 }

solutions-cuda/product-dim.cu

@@ -1,16 +1,6 @@
-// Solution stub for "product-dim".
-// TODO: implement the body below. The signature is derived from
-//       kernel-harnesses/product-dim.cu and must stay in sync with it.
-//
-// Build it (the build system auto-picks this file):
-//   make build/bin/product-dim.exe
-//   ./build/bin/product-dim.exe
-#include <cuda_runtime.h>
-#include <cuda_fp16.h>
-#include <cuda_fp8.h>
-#include <cstdint>
+// product-dim — reduce along dim of an N-D tensor (shape/ndim). Uses the shared dim-reduce framework.
+#include "../kernel-implementation/dim-reduce.cuh"
 
-// Note: all pointer arguments are device pointers.
 extern "C" void solution(const float* input, int dim, float* output, const size_t* shape, size_t ndim) {
-    // TODO: implement product-dim
+    launch_dimreduce<ProdOp>(input, output, shape, ndim, dim);
 }

solutions-cuda/sum-dim.cu

@@ -1,16 +1,6 @@
-// Solution stub for "sum-dim".
-// TODO: implement the body below. The signature is derived from
-//       kernel-harnesses/sum-dim.cu and must stay in sync with it.
-//
-// Build it (the build system auto-picks this file):
-//   make build/bin/sum-dim.exe
-//   ./build/bin/sum-dim.exe
-#include <cuda_runtime.h>
-#include <cuda_fp16.h>
-#include <cuda_fp8.h>
-#include <cstdint>
+// sum-dim — reduce along dim of an N-D tensor (shape/ndim). Uses the shared dim-reduce framework.
+#include "../kernel-implementation/dim-reduce.cuh"
 
-// Note: all pointer arguments are device pointers.
 extern "C" void solution(const float* input, int dim, float* output, const size_t* shape, size_t ndim) {
-    // TODO: implement sum-dim
+    launch_dimreduce<SumOp>(input, output, shape, ndim, dim);
 }