tileinterp

A host interpreter for the Tile IR kernels packages/tileir records: the kernel body itself runs on a machine with no GPU, behind std/gpu's fake device.

std/gpu's fake device answers a launch with a host reference you register for it; it never runs the kernel. This package runs it. interp_ref(prog, grid) lowers the recorded TileProg with tileir's lower and answers a WideRefFn, the function the fake device's table holds, so std/gpu and the host program need no change. The interpreter executes the lowered instruction table, the same table the text renderer and the bytecode writer spell, block by block in order, x fastest.

use std/dtype.{F64}
use std/gpu.{Gpu, with_gpu_fake, alloc, copy_from_host, copy_to_host, launch, buffer_arg, erase}
use std/list
use std/map
use tileir/dev.{Dev, Param, load_cell, store_cell, add, DYN_DIM}
use tileir/prog.{trace3, cells, In, Out}
use tileinterp/interp.{interp_ref}

fn vadd(a: Param[Float], b: Param[Float], out: Param[Float]) -> Unit !Dev =
  store_cell(out, add(load_cell(a), load_cell(b)))

fn run(xs: List[Float], ys: List[Float]) -> Result[List[Float], ForeignError] !Gpu = {
  let a = alloc(F64, len(xs))?
  let b = alloc(F64, len(xs))?
  let out = alloc(F64, len(xs))?
  copy_from_host(a, xs)?
  copy_from_host(b, ys)?
  launch("vadd", len(xs) / 128, [erase(buffer_arg(a)), erase(buffer_arg(b)), erase(buffer_arg(out))])?
  copy_to_host(out)
}

let g = cells([DYN_DIM], [128])
let (prog, _entry) = trace3("vadd", In(F64, g), In(F64, g), Out(F64, g), vadd)
let xs = list.map(range(0, 256), i => to_float(i))
# two blocks of 128: the grid is fixed when the function is made, because
# the fake device does not read the launch's grid
let got = with_gpu_fake(map.from([("vadd", (3, interp_ref(prog, [2])))]), () => run(xs, xs))
# got == Ok([0.0, 2.0, 4.0, ...])

When the grid follows the data, device.with_gpu_interp(progs, body) reads it off every launch instead: install it inside with_gpu_fake (which keeps the buffers) and give it a lookup from kernel name to TileProg.

What it interprets

The integer formats (i1, i4, i8, i16, i32, i64) and every float format std/narrow models: f64, f32, tf32, bf16, f16, f8E5M2, f8E4M3FN, f4E2M1FN and f8E8M0FNU. Element-wise arithmetic and comparisons, select, the shape operations (reshape, broadcast, extract, insert, cat, permute, iota), reductions and scans (left fold), for, loop and if with break, continue and return, tensor and partition views with their five padding values, pointer loads and stores with mask and padding, assume, the conversions between those formats (ftof with its five roundings, itof, ftoi, bitcast, pack / unpack), fpowi, the transcendental functions, fma and the matrix products mmaf, mmaf_scaled and mmai.

A float lane of a narrow format is a Float that is a value of the format, and every operation ends by rounding into the format with std/narrow's own functions, so the interpreter and std/narrow cannot disagree about a format. Add, subtract, multiply, divide, square root, the directed roundings and fma are exact (one rounding of the exact result); the transcendental functions are written here (not taken from packages/tileref, which is what CI holds this package against) and are good to about 1e-13 relative, the tolerance tier's side of docs/tile-backend-design.md 6.6.

A matrix product promises a product and nothing about how the sum over k is grouped, so the interpreter's accumulation is a model: a binary64 accumulator takes one fused step per k, ascending; a narrower one sums the exact products in binary64 and rounds once at the end. Comparisons against it are tolerance-tier ones (docs/tileinterp-design.md 7.1).

Everything else is refused before anything runs, with interp.unsupported_op and the instruction's name: atomics, fences and PDL, strided and gather views, module globals, alloca, assert, print_tko, the pointer casts and the one format std/narrow does not model (f8E5M3FNU) are not planned, because a sequential run gives one legal answer where a device has several and could not check the rest.

Undefined behaviour

The checks are on by default (ub_checks: true), and a launch that reaches undefined behaviour fails with interp.ub rather than answering: an access outside its buffer, a broken inbounds promise, an extract or insert out of range, an integer division by zero, a shift outside the width, a broken overflow assumption, a non-saturating ftoi of NaN or of a value out of range, a false assume, and the use (a store, a branch, a bound, an index) of a lane a masked load or an out-of-bounds view read left unspecified. ub_checks: false runs on with NaN or 0 in such lanes.

A refused launch panics with <kind>: <message> under interp_ref, whose WideRefFn type has no error channel; interp_kernel answers the same as Err(InterpError), and with_gpu_interp as the launch's Err.

How it is held

scripts/tileinterp-diff runs every layer-2 program of scripts/tile-gpu-diff with this interpreter in the device's chair, against the same packages/tileref references and tiers, in CI. The kernels it covers are listed in scripts/tileinterp-diff/covered.txt. The GPU stays the authority: the device ledgers are recorded on the machines that have one. Design and measurements: `docs/tileinterp-design.md` (in Chinese).

Changes between versions: CHANGELOG.md.

API

tileinterp/device

A Gpu handler that answers launch by interpreting the kernel, at the grid the launch names, and hands every other operation to the handler outside it.

device.with_gpu_interp_log

fn with_gpu_interp_log[T](progs: fn(String) -> Option[TileProg], body: fn() -> T !Gpu, ub_checks: Bool = true) -> (T, List[String]) !Gpu

with_gpu_interp, and beside its answer one line per launch in the order they ran: the kernel's name, then ok or the kind of the refusal (interp.unsupported_op, interp.ub, ...). scripts/tileinterp-diff reads these to know which kernels a case actually interpreted.

device.with_gpu_interp

fn with_gpu_interp[T](progs: fn(String) -> Option[TileProg], body: fn() -> T !Gpu, ub_checks: Bool = true) -> T !Gpu

Run body with every launch interpreted (interp.interp_kernel) at the grid the launch names, the kernel found by name through progs. Every other Gpu operation goes to the handler outside, which is normally std/gpu's with_gpu_fake:

with_gpu_fake(map.empty(), () => with_gpu_interp(name => map.get(progs, name), body))

A refused launch answers Err with the interpreter's kind and message, and a name progs does not know answers gpu.no_kernel.

tileinterp/formats

The float formats a tile may hold besides binary64: how a lane of each is held, rounded, flushed and spelled in bits, all through std/narrow.

formats.is_float_format

fn is_float_format(name: String) -> Bool

Whether name is a float format the interpreter models.

formats.is_unmodelled_float

fn is_unmodelled_float(name: String) -> Bool

Whether name is a float format that exists and is not modelled.

formats.float_width

fn float_width(name: String) -> Int

The width of a float format in bits.

formats.mode_of

fn mode_of(attr: String) -> Rounding

The rounding a dialect attribute spells (nearest_even, nearest_away, zero, negative_inf, positive_inf).

formats.round_into

fn round_into(fmt: String, mode: Rounding, x: Float) -> Float

x rounded into format fmt in direction mode. The fp8 and fp4 formats take no direction (their ftof allows nearest-even only, and f8E8M0FNU's zero is the one function std/narrow has for it): the direction is ignored there, because tileir refuses a pair that disagrees before a kernel can be recorded.

formats.nearest

fn nearest(fmt: String, x: Float) -> Float

x rounded into fmt to nearest, ties to even: what every plain operation of the format does to its exact result.

formats.to_bits

fn to_bits(fmt: String, x: Float) -> Int

The pattern of x, a value of fmt, in the low float_width(fmt) bits of an Int (never negative, but for binary64's own sign bit).

formats.of_bits

fn of_bits(fmt: String, b: Int) -> Float

The value of fmt whose pattern is the low bits of b; a NaN pattern keeps its payload (see nan_bits).

formats.flush_in

fn flush_in(fmt: String, x: Float) -> Float

x with the subnormals of fmt replaced by the zero of their sign: the dialect's flush_to_zero. A format with no subnormal range of its own (the fp8 and fp4 ones, which the attribute never applies to) is left.

tileinterp/interp

A host interpreter for the Tile IR kernels packages/tileir records: the lowered instruction table (lower.Kernel) executed block by block over buffers in host memory, as a WideRefFn the fake device can call.

本模块登记的 impl(全程序可见,不必 use):Show[InterpError]

InterpError

type InterpError = { kind: String, message: String }

Why a launch failed: kind is one of interp.unsupported_op (an instruction or format this subset does not model, refused before anything ran), interp.ub (undefined behaviour the kernel reached: an out-of-bounds access, a use of an undefined lane, a false assume, an integer division by zero, a non-saturating ftoi of a value out of range, a broken inbounds promise, a loop past the step limit), interp.bad_launch (the launch does not fit the kernel: buffer or scalar counts, a buffer whose format is not the parameter's) or interp.internal (a defect of this package). message names the instruction.

interp.STEP_LIMIT

const STEP_LIMIT: Int

How many iterations one loop or for may run before the launch is refused as interp.ub: a loop whose exit never comes does not fail on the device, it runs until something kills it.

interp.int_width

fn int_width(name: String) -> Int

The bit width of an integer format, or 0 for any other name.

interp.wrap

fn wrap(w: Int, v: Int) -> Int

v wrapped into the signed range of w bits (an i1 is 0 or 1).

interp.unsupported_in

fn unsupported_in(k: Kernel) -> Option[String]

The first instruction or format of k this subset does not interpret, said as the refusal will say it, or None when every instruction and every format in it is modelled. The interpreter refuses a launch with exactly this message and kind interp.unsupported_op before anything runs, and raises that kind nowhere else.

interp.interp_kernel

fn interp_kernel(k: Kernel, grid: List[Int], dtypes: List[String], ins: List[List[Float]], scalars: List[Float], ub_checks: Bool = true) -> Result[List[(Int, List[Float])], InterpError]

Run k over the grid grid ([gx], [gx, gy] or [gx, gy, gz], missing axes 1) on the buffers ins of formats dtypes and the by-value arguments scalars, in the conventions of std/gpu's WideRefFn. Blocks run one after another, x fastest. The answer is each buffer a store reached, by its position among the buffers, with its whole contents. ub_checks (default on) refuses undefined behaviour the kernel reaches; off, an undefined lane reads NaN or 0 and the run goes on.

interp.interp_ref

fn interp_ref(prog: TileProg, grid: List[Int], ub_checks: Bool = true) -> fn(List[String], List[List[Float]], List[Float]) -> List[(Int, List[Float])]

prog as the fake device runs a kernel: std/gpu's WideRefFn, lowered once here and interpreted over grid on every call. with_gpu_fake does not consult a launch's grid, so the grid is fixed when the function is made (with_gpu_interp in the device module reads it from each launch instead). A refused launch panics with <kind>: <message>, kinds as InterpError has them, because a WideRefFn has no error channel.

tileinterp/math

The binary64 arithmetic the interpreter needs beyond + - * /: the correctly rounded square root, the directed roundings, flush to zero, the integer-valued roundings and thirteen transcendental functions.

Dir

type Dir = ToNearest | Downward | Upward

A rounding direction for the operations that take one.

math.INF

const INF: Float

The positive infinity.

math.NAN

const NAN: Float

A quiet NaN.

math.SIGN_BIT

const SIGN_BIT: Int

The sign bit of a binary64 pattern.

math.LN_2

const LN_2: Float

ln 2 to the nearest Float.

math.PI_F

const PI_F: Float

pi to the nearest Float.

math.bits_of

fn bits_of(x: Float) -> Int

The bit pattern of x (binary64), and its inverse.

math.float_of

fn float_of(b: Int) -> Float

The binary64 value whose bit pattern is b.

math.is_nan

fn is_nan(x: Float) -> Bool

x is a NaN.

math.is_finite

fn is_finite(x: Float) -> Bool

x is neither infinite nor a NaN.

math.sign_bit

fn sign_bit(x: Float) -> Bool

The sign bit of x is set (so -0.0 and negative NaNs count).

math.is_subnormal

fn is_subnormal(x: Float) -> Bool

x is a nonzero subnormal number.

math.signed_zero

fn signed_zero(x: Float) -> Float

The zero with the sign of x.

math.flush

fn flush(x: Float) -> Float

x with subnormals replaced by the zero of their sign: the dialect's flush_to_zero, which applies to the inputs and to the result.

math.next_up

fn next_up(x: Float) -> Float

The next binary64 value above x (toward +inf). +inf and NaN are their own successor; the successor of -0.0 and of +0.0 is the least positive subnormal.

math.next_down

fn next_down(x: Float) -> Float

The next binary64 value below x (toward -inf).

math.pow2

fn pow2(k: Int) -> Float

2^k as a Float, for any k binary64 can hold (subnormal results included).

math.scale

fn scale(x: Float, k: Int) -> Float

x * 2^k by at most three exact steps, so that a result in the subnormal range is rounded once.

math.exponent_of

fn exponent_of(x: Float) -> Int

The binary exponent of a finite nonzero x: x = m * 2^e with m in [1, 2).

math.two_prod

fn two_prod(a: Float, b: Float) -> (Float, Float)

Dekker's product: (p, e) with p = fl(a * b) and a * b = p + e EXACTLY, provided neither the product nor the splits overflow and the product does not underflow. Callers scale into range first.

math.two_sum

fn two_sum(a: Float, b: Float) -> (Float, Float)

Knuth's TwoSum: (s, e) with s = fl(a + b) and a + b = s + e exactly, for finite operands whose sum does not overflow.

math.odd_of

fn odd_of(s: Float, e: Float) -> Float

The ROUND-TO-ODD value of the exact sum s + e, where s = fl(x + y) and e is the error TwoSum reports: s itself when the sum was exact or its significand is already odd, otherwise the neighbour of s on the side of e, whose significand is odd. Rounding that value to nearest, or in any direction, at a format with at most 51 significand bits gives what rounding the EXACT sum once would (Boldo and Melquiond, "Emulation of FMA and correctly rounded sums: proved algorithms using rounding to odd", IEEE Trans. Computers 2008). It is the one device that lets binary64 stand in for a wider format in the fused operations below.

math.add_odd

fn add_odd(x: Float, y: Float) -> Float

x + y rounded to odd (see odd_of); x and y finite.

math.fma_f64

fn fma_f64(a: Float, b: Float, c: Float) -> Float

a * b + c with ONE rounding to nearest, for binary64 operands: the fused multiply-add the dialect's fma performs. Boldo and Melquiond's algorithm: the exact product as uh + ul (Dekker), the sum c + uh as th + tl (TwoSum), v = RO(tl + ul), and RN(th + v). It needs the product where Dekker's split cannot overflow and ul cannot underflow (|a| and |b| below 2^500, |c| below 2^1000, |a b| between 2^-900 and 2^900). Outside that range, and for an infinity, a NaN or a zero product, it answers `a b + c`: exact where an operand is special, two roundings in the extreme-range case (docs/tileinterp-design.md 7.3).

math.sqrt_rn

fn sqrt_rn(x: Float) -> Float

The correctly rounded square root of x (ties cannot occur). NaN for a negative or NaN operand, x itself for either zero and for +inf.

x = m * 4^k with m in [1, 4), Newton's iteration for sqrt(m), then the exact test that picks the right neighbour: y is sqrt(m) rounded to nearest iff pred(y) * y < m <= y * succ(y), both products formed exactly by two_prod (the quarter-ulp-squared term that separates the true midpoint square from these products is below the granularity both sides of each comparison share, so it cannot change the answer).

math.directed

fn directed(r: Float, err: Int, dir: Dir) -> Float

The value r, the round-to-nearest result of an exact operation, moved to the directed result: err is the sign of exact - r (-1, 0 or 1).

math.add_dir

fn add_dir(a: Float, b: Float, dir: Dir) -> Float

a + b rounded in direction dir.

math.mul_dir

fn mul_dir(a: Float, b: Float, dir: Dir) -> Float

a * b rounded in direction dir.

math.div_dir

fn div_dir(a: Float, b: Float, dir: Dir) -> Float

a / b rounded in direction dir.

math.floor_f

fn floor_f(x: Float) -> Float

floor(x): the largest integer-valued Float not above x.

math.ceil_f

fn ceil_f(x: Float) -> Float

ceil(x): the smallest integer-valued Float not below x.

math.rem_f

fn rem_f(a: Float, b: Float) -> Float

The IEEE remainder of truncated division, fmod: a - n * b with n the quotient truncated toward zero, computed exactly by long division on the binary exponents (the result is always representable).

math.exp_f

fn exp_f(x: Float) -> Float

e^x. The argument is divided by 2^10, the Taylor series of the small result summed to fourteen terms, and the answer squared ten times. Each squaring doubles the relative error, so the result is good to about 1024 ulps; that is 2e-13 relative, inside the tolerance tier by eight orders of magnitude.

math.exp2_f

fn exp2_f(x: Float) -> Float

2^x, exact at integer arguments.

math.log_f

fn log_f(x: Float) -> Float

ln x: x = m * 2^e with m in [1, 2), then Halley's iteration for the y with e^y = m, which converges cubically from the two-term guess.

math.log2_f

fn log2_f(x: Float) -> Float

log2 x, exact at powers of two.

math.rsqrt_f

fn rsqrt_f(x: Float) -> Float

1 / sqrt(x).

math.tanh_f

fn tanh_f(x: Float) -> Float

tanh x as 1 - 2 / (e^(2|x|) + 1) with the sign put back; past 20 the answer is 1 to the last bit. Small arguments lose relative accuracy to the subtraction but keep an absolute error near 1e-16, which is what the tolerance tier measures.

math.sinh_f

fn sinh_f(x: Float) -> Float

sinh x.

math.cosh_f

fn cosh_f(x: Float) -> Float

cosh x.

math.sin_f

fn sin_f(x: Float) -> Float

sin x. NaN for an infinite argument.

math.cos_f

fn cos_f(x: Float) -> Float

cos x. NaN for an infinite argument.

math.tan_f

fn tan_f(x: Float) -> Float

tan x.

math.atan_f

fn atan_f(x: Float) -> Float

atan x.

math.atan2_f

fn atan2_f(y: Float, x: Float) -> Float

atan2(y, x), the angle of the point (x, y), with C's signed-zero and infinity conventions.

math.pow_f

fn pow_f(x: Float, y: Float) -> Float

x^y with C's pow special cases: x^0 = 1 and 1^y = 1 for every operand (NaN included), a negative base only to an integer power, zero and infinity by the sign and parity rules.

tileinterp/mma

The fused operations on plain lists: fma and the three matrix products, each with the accumulation order and the rounding the interpreter models, and the reason for it written beside it.

mma.fma_in

fn fma_in(fmt: String, a: Float, b: Float, c: Float) -> Float

a * b + c in format fmt, rounded once.

Binary64: Boldo and Melquiond's emulation (math.fma_f64). A narrower format: the product of two operands of at most 24 significand bits is exact in binary64, c is added by rounding to odd, and that value is rounded into fmt; rounding to odd first loses nothing for a format with at most 51 significand bits, so the second step is the only rounding.

mma.mma_float

fn mma_float(out: String, m: Int, k: Int, n: Int, a: List[Float], b: List[Float], acc: List[Float]) -> List[Float]

acc + a * b for an m by k matrix a, a k by n matrix b and an m by n accumulator, all row major, accumulating into format out. See the header for the order and the roundings.

mma.mma_scaled

fn mma_scaled(m: Int, k: Int, n: Int, v: Int, a: List[Float], b: List[Float], acc: List[Float], sa: List[Float], sb: List[Float]) -> List[Float]

The block-scaled product of Ops.td's MmaFScaledOp: acc plus the sum over p of (a[i][p] * sa[i][p / v]) * (b[p][j] * sb[p / v][j]), where sa is m by k / v and sb is k / v by n. The scales are powers of two (f8E8M0FNU), so every factor and product is exact in binary64, and the sum is accumulated wide and rounded into f32 once, like the unscaled product into a narrower accumulator.

mma.mma_int

fn mma_int(m: Int, k: Int, n: Int, a: List[Int], b: List[Int], acc: List[Int]) -> List[Int]

The integer product of mmai: acc plus the sum of the products of signed i8 operands, wrapped to the i32 accumulator. Integer addition is associative modulo 2^32, so the order is not a model, it is exact.

tileinterp/selftest

Kernels recorded with packages/tileir and run through the interpreter, each answer written down by hand. No function here is public.