Operator Variations & Capacity Zones

One kernel, four flavors — and a reality check on size

Once vector add works, swapping the operator gets you sub, mul, and div for free. The launch geometry, memory layout, and verification stay identical — only the ALU op changes. But the cost isn't always identical:

Op	Body change	Notes
Add	`C[i] = A[i] + B[i]`	Baseline.
Sub	`+` → `-`	Same throughput as add — single FP unit.
Mul	`+` → `*`	Uses FMA pipelines, same throughput.
Div	`+` → `/`	5–10× slower due to non-pipelined reciprocal. Guard against `B[i] == 0` to avoid NaN/inf.

VRAM Capacity Zones — what happens when you ask for too much

Zone	What happens	Escape hatch
🟢 Green: fits in VRAM	Full-speed, data resident on the GPU	No action needed
🟡 Yellow: UVM paging	10–50× slowdown as Unified Memory migrates pages	Process in chunks
🔴 Red: hard alloc fail	`cudaMalloc` returns `cudaErrorMemoryAllocation`	Slice-and-stream loop
🟣 Crimson: grid limit	`gridDim.x` exceeds 2³¹⁻¹	Outer host loop, base-index arg

Apple Silicon shifts the zones: there's no separate VRAM, so 'Green' extends to the entire unified pool minus what the OS and other apps are holding. The Yellow zone (paging to swap) still exists but kicks in much later.

Code

Generic 1-kernel-per-op pattern (CUDA)·cuda

// Use a template + functor or the preprocessor.
// Functor version is debugger-friendly:

struct AddOp { __device__ float operator()(float a, float b) const { return a + b; } };
struct SubOp { __device__ float operator()(float a, float b) const { return a - b; } };
struct MulOp { __device__ float operator()(float a, float b) const { return a * b; } };
struct DivOp { __device__ float operator()(float a, float b) const {
    return b == 0.f ? 0.f : a / b;   // guarded
} };

template <typename Op>
__global__ void elementwise(const float* A, const float* B, float* C, size_t N, Op op) {
    size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) C[i] = op(A[i], B[i]);
}

// Launch with whichever Op you want; nvcc inlines the call:
// elementwise<<<blocks, threads>>>(dA, dB, dC, N, MulOp{});

Slice-and-stream when data exceeds VRAM·cuda

// Process N elements in chunks of CHUNK to fit in VRAM.
const size_t CHUNK = 1 << 24;        // 16M elements ≈ 64 MB per buffer
for (size_t off = 0; off < N; off += CHUNK) {
    size_t this_n = std::min(CHUNK, N - off);

    cudaMemcpyAsync(dA, hA + off, this_n*sizeof(float),
                    cudaMemcpyHostToDevice, stream);
    cudaMemcpyAsync(dB, hB + off, this_n*sizeof(float),
                    cudaMemcpyHostToDevice, stream);

    int threads = 256;
    int blocks  = (this_n + threads - 1) / threads;
    vec_add<<<blocks, threads, 0, stream>>>(dA, dB, dC, this_n);

    cudaMemcpyAsync(hC + off, dC, this_n*sizeof(float),
                    cudaMemcpyDeviceToHost, stream);
}
cudaStreamSynchronize(stream);

Exercise

Time the four operators (add, sub, mul, div) on the same 64M-element vectors. You should see add ≈ sub ≈ mul, and div noticeably slower (often 3–7× depending on GPU generation). Then pick a vector size 2× your VRAM and watch the kernel either crash with an alloc error (CUDA) or slow to a crawl (Apple Silicon paging). Implement the slice-and-stream version and confirm it runs at the same per-element rate as the in-VRAM case.