Skip to content

saxpy

saxpy computes out = alpha * x + y for vectors x, y and out of length N and a scalar alpha. It is the smallest kernel in the repository, and the first one phobos-kbench runs.

On an RTX 2080 SUPER at n = 2^25, phobos-kbench measured this kernel moving 430 GB/s against 433 GB/s for cublasSaxpy_v2, 99% of it. saxpy is memory-bound, so bandwidth, not GFLOP/s, is its measure. From results/results.csv, committed 2026-10-08; all kernels are on the front page.

@autotune(BLOCK in [16, 8192])
kernel saxpy(x: tensor<f32>[N], y: tensor<f32>[N], out: tensor<f32>[N], alpha: f32) {
let base = ctx.block.x * BLOCK
out[base :+ BLOCK] = alpha * x[base :+ BLOCK] + y[base :+ BLOCK]
}

examples/saxpy_fp32.phOpen full view

The saxpy kernel, stepped through as the compiler built it for sm_75A 4 grid of CTAs with 256 threads each. The drawing shows the tensors x, y, out in global memory with the picked CTA's current boxes, shared memory as a byte ruler of 128 bytes with each buffer at its planned offset, the barriers kept and elided, and how the CTA's threads and warp 0's lanes take part. The readout below describes each step.

CTA 0
  • global memory
  • shared memory
  • threads and lanes
  • extra 128-byte lines, masked elements

Your first kernel reads it line by line. This page looks at what happens when it runs.

The kernel is launched as a one-dimensional grid. Every CTA runs the same two lines, and ctx.block.x tells each one which part of the vectors is its own: CTA i handles elements i * BLOCK up to (i + 1) * BLOCK. The launcher starts N / BLOCK CTAs, rounded up, so the slices tile the vectors with no gaps and no overlap.

x: [ CTA 0: 0 .. BLOCK ) [ CTA 1: BLOCK .. 2*BLOCK ) ... [ last CTA, may be short )

Inside a CTA, the BLOCK elements are shared among its threads, 256 by default (see @launch). With BLOCK = 4096, each thread handles 16 elements. You never write that split; the compiler derives it from the slice length and the thread count.

The kernel declares no tile: the expression reads slices of x and y and the assignment writes a slice of out. Expressions over slices are covered in Elementwise math and reductions.

@autotune(BLOCK in [16, 8192]) makes BLOCK a tuning parameter: 16, 32, 64 and so on up to 8192, doubling at each step. A small BLOCK gives many small CTAs; a large one gives fewer CTAs that each do more. Which one is fastest depends on the card, so phobos-kbench measures them all and keeps the best. See Autotuning and @autotune.

saxpy does one multiply and one add per element, against two loads and a store. Its speed is set by how fast the card moves memory, not by how fast it computes.

N is only known at launch, and it need not be a multiple of BLOCK. The last CTA’s slices can then run past the end of the vectors. Because the slice offset comes from ctx.block and the extent N is dynamic, the compiler masks the access against the runtime length: elements past the end read as zero, and their stores are skipped. Tensor extents are assumed to be multiples of 4.

A kernel that can promise whole slices removes the mask with @aligned, as the gemm example does. The rules are in Slicing, bounds and masking and the spec’s Bounds and masking.

Terminal window
phobos-compile examples/saxpy_fp32.ph -o saxpy.ptx
cargo run -r -p phobos-kbench -- --bench saxpy_fp32

The first needs no GPU. The second autotunes BLOCK on the GPU and prints the throughput of the winner.