README.md

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 only. 2d@2d→{M,N}, 2d@1d→{M}, 1d@2d→{N}, 1d@1d→scalar
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; TL_BATCH_MATMUL_BIAS overrides the batched-matmul one at runtime.

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(). #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/cuda/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-Instruct GGUF model 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 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?"

Decode throughput on an RTX 3090 has gone from 3.5 to ~330 tok/s over the course of the CUDA decode-path work — about 1.35x off llama.cpp on the same shape. See bench/cuda/speed/bench_qwen_decode.cpp and docs/performance-notes.md for the measurement methodology and the full tuning history.

License

MIT license (c) 2026 yhirose