cpp-tensorlib

A C++17 header-only, cross-platform tensor library with zero third-party dependencies — CPU (own SIMD kernels) and GPU (own Metal / CUDA / WGSL kernels) on macOS, Linux, Windows and the browser, from a single #include <tensorlib.h>. No BLAS, no cuBLAS/cuDNN, no CUTLASS: every accelerated kernel — GEMM, GEMV, attention, quantized (int4) matmul — is hand-written in this repo, so the only thing a consumer links against is the OS itself (and the CUDA driver, which is dlopen'd at runtime, not linked).

As a full end-to-end proof of that promise, the library also drives a real local-LLM chat pipeline — GGUF loader, BPE tokenizer, and a Qwen2.5 decoder, all on the same own-kernel, zero-dependency foundation (see LLM inference below).

• Header-only, C++17 minimum (builds cleanly through C++23) — #include <tensorlib.h>, nothing to link or build separately.
• Cross-platform: macOS (Accelerate + Metal), Linux/Windows (own BLIS-style CPU microkernels — scalar/AVX2/NEON — and own CUDA kernels), and WebAssembly (own WGSL kernels over WebGPU). One dispatch seam, no platform #ifdefs in consumer code.
• Zero third-party dependencies: no OpenBLAS, cuBLAS, cuDNN, or CUTLASS anywhere in the accelerated paths (see Dependency policy).
• MLX-style lazy evaluation: ops build a graph; tl::eval() evaluates multiple arrays in one topological pass, with build-time peephole fusion (every node carries an affine epilogue, so scalar chains and dot-then-scale collapse to a single dispatch).
• Switchable backend via tl::use_cpu() / tl::use_gpu() / tl::use_auto() (auto picks CPU or GPU per op from measured per-kernel-class thresholds).
• Zero-copy views (transpose / reshape / slice), numpy broadcast rules.
• Data types: float (primary), plus bf16/int4 storage dtypes for the bandwidth-bound LLM decode path (see the milestones).

Design informed by silarray (the macOS-only experiment this succeeds), rebuilt for three platforms with a single dispatch seam the backends plug into. Serves as the matrix foundation for culebra's Tensor type (culebra keeps autograd/VJP; this library owns the graph, fusion, and execution).

Example

#include <tensorlib.h>

auto a = tl::array::ones({1000, 1000});
auto b = tl::array::ones({1000, 1000});

auto c = a.dot(b) * 0.5f + 1.0f;   // lazy: one fused GEMM+epilogue dispatch
auto d = (a + b).relu();           // also lazy

tl::eval(c, d);                    // evaluate both in one pass

tl::use_gpu();                     // Metal (macOS), CUDA (Linux/Windows), WebGPU (wasm)
auto e = a.dot(b);
tl::use_auto();                    // CPU or GPU per op, by measured size

API reference

#include <tensorlib.h> pulls in the whole library; everything lives in namespace tl. The CMake target is INTERFACE, so it deliberately does not impose a standard on you — set C++17 or newer yourself.

tokenizer.h and gguf.h are not in the umbrella header; include them directly if you want them.

Creating arrays

Shapes are tl::shape_t (= std::vector<int64_t>), so braced literals work. {} is a valid shape — it makes a rank-0 scalar.

array::zeros(shape) / ones(shape)	filled with 0 / 1
array::full(shape, v)	filled with v
array::empty(shape)	allocated, uninitialized
array::from(vector<float> v)	1-d from host data (copies)
array::from(vector<float> v, shape)	same, reshaped; throws if the count mismatches
array()	undefined array (defined() is false)
tl::concat(vector<array>)	join along axis 0; other dims must match

Inspecting

shape(), strides(), rank(), size(), contiguous(), defined(), dt(). Note it is rank(), not ndim().

Reading elements: at({i, j}) (no bounds checking), item() for a single-element array, raw() for a const float*, data() for a mutable float*. All of these evaluate the graph first, and flush any pending GPU work — mixing them with lazy ops is safe.

There is no printing. No operator<<, no to_string — format arrays yourself from raw().

For tests: allclose(a, b, rtol = 1e-5f, atol = 1e-6f) and array_equal(a, b) (the same thing with zero tolerance). Both return false on a shape mismatch.

Operations

Everything below is lazy unless noted.

Group	API
Arithmetic	+ - * / in array⊗array, array⊗float and float⊗array forms; pow(a, b), pow(a, s)
Comparison	> < >= <= == !=, same three forms — result is an f32 mask of 1.0/0.0; where(cond, a, b)
Unary	.exp() .log() .sqrt() .sigmoid() .relu()
Softmax	.softmax() — last axis only, numerically stable, throws on rank 0
Matmul	a.dot(b) — rank 1 or 2, or batched rank ≥ 3 (same rank, identical leading batch dims, no broadcast). 2d@2d→{M,N}, 2d@1d→{M}, 1d@2d→{N}, 1d@1d→scalar, [...,M,K]@[...,K,N]→[...,M,N]
Axis reductions	.sum(axis, keepdims=false), .mean(...), .max(...), .argmax(...) — any axis, negatives count from the end; argmax returns indices as f32
Scalar reductions	float sum(), float max(), float mean(), int64_t argmax() — eager, return host values
Views	.transpose(), .transpose({axes}), .reshape(shape), .slice(start, count) (axis 0 only), .clone()
Broadcast VJP	.sum_to(shape) — reduce back to a shape that broadcasts to this one
In-place	.add_(b) — eager, mutates shared storage
Model ops	array::attn_decode(q, K, V, scale), array::rope(x, pos, base=10000.0f), array::rmsnorm(x, weight, eps=1e-5f), array::silu(x), array::swiglu(gate, up)

Broadcasting follows numpy's trailing-dimension rules; an incompatible pair throws. attn_decode takes q [H,D], K/V [H,ctx,D] and returns [H,D].

Lazy evaluation

Ops build a graph. Nothing runs until something forces it:

auto c = (a.dot(b) * 0.5f + 1.0f).relu();   // nothing computed yet
c.eval();                                    // now it runs
tl::eval(c, d, e);                           // or: several roots, one pass

.eval() (member) is the idiomatic form; the free tl::eval(...) is variadic and evaluates multiple roots in a single topological pass, which is what you want when their graphs share subexpressions.

Evaluation also happens implicitly on data(), raw(), item(), at(), the scalar reductions, add_(), and when a view is taken of a lazy source.

Choosing a backend

tl::use_cpu();    // default
tl::use_gpu();    // Metal / CUDA / WebGPU, per platform
tl::use_auto();   // per op, by measured size thresholds
bool ok = tl::gpu_available();   // compiled in AND a device is present

Selection is a process-wide switch, and which GPU backend you get is decided at compile time. When no device is present, use_gpu() and use_auto() silently run on the CPU. The auto thresholds are measured per kernel class and per backend; the matmul one is re-derived on the host at the first auto-mode graph eval (a four-point CPU-vs-GPU census, a few ms; TL_AUTO_MATMUL pins it, TL_AUTO_TRACE=1 prints the census) and TL_BATCH_MATMUL_BIAS overrides the batched-matmul one. The CPU gemm splits across threads only when each would get a floor of multiply-adds that is derived from the pool's measured wake-up on first use (TL_CPU_MIN_WORK pins it; misc/census_pool_latency.cpp shows the arithmetic); below that it runs on the calling thread.

The GPU ops are written once, in gpu_ops.h, over views (gpu::span: a device handle and a byte offset) and a kernel ABI every backend shares (gpu_abi.h). A backend is a device core — memory, a kernel table, one dispatch — plus kernel source, and the ops it runs its own way as members of its own struct; an op a backend has no kernel for answers false and the evaluator falls back to the CPU. gpu::census(kernel) counts launches, so a test can tell the two apart. No machine here runs CUDA kernels, so tools/cuda_trace records what the CUDA backend asks of the driver (kernel, grid, every argument) against a stand-in libcuda in a Linux container, and tools/cuda_trace/compare.sh <ref> diffs that across a change to the backend's host side. -DTENSORLIB_HOST_GPU=ON selects a reference backend whose device is the CPU (gpu_host.h), which runs that whole layer and its conformance tests with no GPU. docs/backends.md covers the layers, the kernel ABI, and how to add an op or a backend.

Profiling

tl::profile (profile.h) says where evaluation goes: a profile::scope labels a block, the evaluator opens one per op under it, and every kernel launch, transfer and blocking wait a backend performs lands under the innermost open scope. profile::start() / stop() bracket a session, rows() returns it, report(stderr) prints it as a table; TL_PROFILE=1 does all three for a whole process. CUDA stamps each launch with events, so its rows carry a device time; Metal knows the GPU time per command buffer (summarize().batch_device_us); WebGPU counts launches. A row's device_timed says how many of its launches the device time covers.

Storage dtypes

Compute and results are always f32. bf16 and q4 are weight-container formats for the bandwidth-bound decode path, not general dtypes:

auto w = weights.to_bf16();   // f32 -> bf16 (round-to-nearest-even)
auto q = weights.to_q4();     // f32 [K,N] -> grouped int4; needs K % 32 == 0
auto f = q.to_f32();          // back to f32 (a no-op if already f32)

Ops other than the native decode GEMV widen transparently. Direct element access requires f32 — raw(), data(), at() and item() throw on a bf16/q4 array, so call to_f32() first.

Things that will bite you

• data() requires a contiguous array; raw() does not. Calling data() on a transposed view throws — clone() it first.
• Views alias their base. transpose/reshape/slice share storage, so a write through data() is visible through every view of it — including from unevaluated graphs still holding that array.
• add_() mutates shared storage and is not safe on an array that is still feeding a lazy graph.
• Evaluation is single-threaded by design. There is no documented support for sharing an array across threads.

Opt-in headers

#include <tokenizer.h> gives tl::tokenizer(gguf_path) with encode(), decode(), bos_id(), eos_id() — GPT-2 byte-level BPE, read from the model's own GGUF metadata, so it serves any model of that family (Qwen2 is the one the tests pin against). #include <gguf.h> gives tl::gguf::model(path) — an mmap'd reader with tensor(name), tensors(), kv(key), metadata(). Both are used by the LLM pipeline below.

Backends

Platform	CPU	GPU
macOS	Accelerate (vDSP / vForce / CBLAS) ✅	Metal / MSL (#embed JIT), STEEL SGEMM ✅
Linux / Windows	own BLIS-style microkernels (scalar / AVX2 / NEON, runtime-dispatched) ✅	own CUDA kernels, PTX JIT'd by a dlopen'd driver ✅
WebAssembly	own scalar microkernels ✅	own WGSL kernels over WebGPU ✅
any	reference strided implementation (oracle + fallback) ✅	—

Everything reduces to strides through one index walker, so views, broadcast and transposed operands share a single code path. Accelerated backends replace it per-op at the graph::eval_one dispatch seam; the seam carries no platform #ifdefs (non-Apple builds get inline stubs).

All three GPU backends cover the array surface — GEMM, elementwise unary/binary, rank-2 broadcast, row reductions (softmax / row sum / row max) — each with the same fused affine epilogue, and each validated against the reference implementation as an oracle. The LLM decode kernels (GEMV, attention, RoPE, int4 dequant-matmul) are CUDA-only; on Metal and WebGPU those ops decline at the seam and run on the CPU.

The WebGPU backend is built with emcc --use-port=emdawnwebgpu -DTENSORLIB_WEBGPU, and differs from the native two in two ways worth knowing before you target it:

• Chrome in practice. It needs JSPI to keep flush() and sync_to_host() synchronous. Firefox has it behind a flag and Safari has not shipped it; there, device acquisition fails, available() stays false, and every op routes to CPU — correct, just not accelerated. That is also the fallback if the page hands in no device.
• No STEEL-equivalent tiling. The auto-mode thresholds are calibrated against measurements, but the WGSL SGEMM is a straightforward tiled kernel; macOS gets a considerably more tuned one.

CI runs the wasm suite on every push — test/wasm/ builds the same oracle suite as the native backends and runs it headless via Deno (a browser is not required, locally either), asserting a per-kernel dispatch census so a silent all-CPU fallback fails the build rather than passing green.

Dependency policy

Zero third-party dependencies, on every platform. No OpenBLAS, no cuBLAS/cuDNN, no CUTLASS — every accelerated kernel (CPU SIMD, CUDA GEMM/ GEMV/attention, int4 dequant-matmul) is hand-written in this repo. What each platform links against:

• macOS — only OS frameworks (Accelerate, Metal, Foundation).
• Linux / Windows (CPU) — nothing; own microkernels compiled straight into the consumer binary.
• Linux / Windows (CUDA) — nothing linked. Kernels are compiled to PTX at build time (nvcc, build-time only — not a runtime dependency) and the CUDA driver is dlopen'd (LoadLibrary on Windows) at run time, so a binary built with -DTENSORLIB_CUDA=ON still runs — falling back to CPU — on a machine with no NVIDIA driver at all.
• WebAssembly — nothing linked. The WGSL kernels are committed as a C string and the WebGPU entry points come from Emscripten's bundled emdawnwebgpu port, so the toolchain is the only requirement and it is a build-time one.
• Tests only — doctest is vendored (not a runtime dependency of the library itself).

Build and test

cmake -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build
ctest --test-dir build --output-on-failure   # cpu / gpu / auto modes
./build/tensorlib_bench                       # micro-benchmarks

On Linux/Windows with an NVIDIA GPU, add -DTENSORLIB_CUDA=ON (needs nvcc on PATH at build time only — see Dependency policy):

cmake -B build -DCMAKE_BUILD_TYPE=Release -DTENSORLIB_CUDA=ON
cmake --build build
ctest --test-dir build --output-on-failure

Requires C++17 or newer — the headers use inline variables, std::optional, and structured bindings, but nothing from C++20/23, so any g++ >= 11 / clang++ >= 13 works (Apple Clang too). The bundled CMake build compiles the tests at C++23, but consuming the headers only needs C++17. Exception: the macOS/Metal backend #embeds its shader source, so building on macOS additionally needs a #embed-capable compiler (Clang 19+); non-Apple builds never reach that #embed. The WebGPU backend has the same problem and solves it the other way — its WGSL is committed as a generated C string (kernels/tensorlib_webgpu_wgsl.inc, regenerated by kernels/gen_wgsl_inc.sh when the .wgsl changes), so it needs no #embed and no build step.

For a WebAssembly build, see test/wasm/build.sh — a flat emcc line with --use-port=emdawnwebgpu -DTENSORLIB_WEBGPU, plus -sJSPI=1, which is what lets flush() keep a synchronous signature.

LLM inference

bench/models/chat_qwen.cpp is the end-to-end proof of the header-only, zero-dependency, cross-platform claims above: it loads a real Qwen2.5, Qwen3 or Qwen3.5 GGUF model (F16, BF16, or quantized as Q4_K / Q6_K / Q8_0 blocks, which stay quantized on the device: Qwen3.8-27B Q4_K_M, 17.7 GiB, runs on a 24 GB RTX 3090) and chats, using nothing but this library's own kernels — own GGUF v3 reader (include/gguf.h), own GPT-2-byte-BPE tokenizer (include/tokenizer.h), and own CUDA prefill/decode kernels. No GGML, no llama.cpp, no third-party inference runtime anywhere in the path.

cmake -B build -DCMAKE_BUILD_TYPE=Release -DTENSORLIB_CUDA=ON
cmake --build build --target tensorlib_chat_qwen
./build/tensorlib_chat_qwen path/to/qwen2.5-0.5b-instruct.gguf "What is the capital of France?"

On an RTX 3090, Qwen2.5-0.5B (F16 GGUF, bf16 weight storage):

Path	tok/s
decode, CUDA-graph replay	~460
decode, imperative (no graph capture)	~305
prefill, 512-token prompt	~44000
prefill, 2048-token prompt	~43000

Decode started at 3.5 tok/s and reached ~460 by removing host-side graph construction, then per-launch overhead, then GEMV inefficiency. Against llama.cpp on the same files and GPU, decode is at parity or ahead for Qwen2.5, Qwen3 and Qwen3.5, with or without 2048 tokens cached, and prefill is at parity at 512 tokens and up to 1.3x behind with 2048 cached (both sum the prefill's products in f16 by default): docs/llama-cpp-comparison.md has the numbers, what differs between the two runs, and how to reproduce them. Quantized files prefill their Q4_K and Q8_0 projections on the int8 tensor cores, the activation rounded to 8 bits as llama.cpp's MMQ does (TL_NO_MMQ=1 takes an f16 path instead), and are level with or ahead of llama.cpp there, Qwen3.8-27B included.

Prefill is a different regime and was slower than decode for a long time, because the prompt was fed through the decode path one token at a time. Its tokens share one set of weights, so batching turns the projections into GEMMs rather than GEMVs and the attention into one tiled causal pass — worth ~50x over replaying the decode graph per prompt token, and ~170x over the array-graph path chat used originally.

bench/cuda/speed/bench_qwen_decode.cpp and bench_qwen_prefill.cpp produce these numbers and also guard them: every prefill path is checked against the array-graph reference for the same first generated token, and the decode comparison is greedy-exact. Measure with min-of-rounds — a single run on WSL2 swings ±20 tok/s with clock ramp.

License

MIT license (c) 2026 yhirose