| File | Lines | Bytes | Role |
|---|---|---|---|
include/tensor.h |
751 | 25 KB | Public API — single include for all consumers |
src/tensor.c |
1,001 | 32 KB | All matrix + parameter implementations |
src/random.c |
128 | 3.7 KB | xorshift32 PRNG, uniform & Xavier fill |
tests/test_tensor.c |
1,271 | — | 56-test suite, zero external dependencies |
examples/basic_usage.c |
203 | — | Annotated walkthrough of every API call |
Makefile |
72 | — | build / run / demo / valgrind / clean |
README.md |
345 | — | Full API reference + 15-part roadmap |
| Total | 3,771 | ~62 KB | Zero external dependencies |
Test result: 56 / 56 PASS — zero compiler warnings under -Wall -Wextra -Werror -pedantic
typedef struct {
size_t rows;
size_t cols;
float *data; /* contiguous row-major storage */
} Matrix;Why row-major?
Three concrete reasons were evaluated and documented:
- C language native —
arr[r][c]in C is already row-major; our layout matches the language model, eliminating mental translation. - Cache locality — the dominant neural network operation is iterating over a complete row (dot products, activation outputs). Row-major puts consecutive row elements in adjacent cache lines, maximising L1 hit rates.
- SIMD readiness — AVX2/NEON intrinsics (
_mm256_fmadd_ps) operate on contiguous 256-bit/128-bit chunks. Inner loops in row-major layout are mechanically vectorisable.
Element indexing formula: data[r * cols + c]
Matrix* matrix_create(size_t rows, size_t cols); /* calloc — OS-guaranteed zero init */
void matrix_free(Matrix *m); /* poisons pointer after free */Key decisions:
callocnotmalloc+memset—callocgets OS-guaranteed zero pages. On Linux this is often free (zero pages from the kernel).memsetwould re-touch pages we already paid for.- Pointer poisoning — after
free(m->data), the implementation setsm->data = NULL,m->rows = 0,m->cols = 0before freeing the struct. This converts use-after-free into an immediate NULL-dereference crash, which is far easier to diagnose than silent memory corruption. matrix_free(NULL)is a no-op — matches the semantics of the standardfree(), making cleanup code at error paths unconditionally safe.
The library enforces a strict no-silent-failure rule via three mechanisms working together:
| Mechanism | When used | Example |
|---|---|---|
NULL return |
Any function that allocates, on failure | matrix_create OOM |
stderr message |
Every detected error with context | "A is (2×3), B is (3×2)" |
| Early return / no-op | Mutating functions with invalid input | matrix_set OOB write |
Two internal guard macros keep this DRY:
#define REQUIRE_NON_NULL(p, name, action) /* null check + stderr + action */
#define REQUIRE_SAME_SHAPE(A, B, action) /* shape check + stderr + action */These macros embed __func__ so every error message automatically names the function that detected it. Shape mismatch messages always print both actual shapes, e.g.:
[mini_tensor] ERROR: matrix_matmul — inner dimensions do not match:
A is (2×3), B is (4×2). Required: A->cols == B->rows.
Algorithm: xorshift32 (Marsaglia 2003)
static uint32_t rng_state = 2463534242u;
static uint32_t rng_next_u32(void) {
uint32_t x = rng_state;
x ^= x << 13;
x ^= x >> 17;
x ^= x << 5;
rng_state = x;
return x;
}Why xorshift32?
| Property | Value |
|---|---|
| Period | 2³² − 1 ≈ 4 billion samples |
| Operations | 3 XOR-shifts — no multiplication, no table |
| Statistical quality | Passes Diehard tests; adequate for weight init |
| Reproducibility | rng_seed(42) gives identical sequences across runs |
| Dependencies | None — pure integer arithmetic |
Two fill functions built on top:
rng_fill_uniform(data, n, min, max)— linear interpolation from [0,1) into [min,max]rng_fill_xavier(data, n, fan_in, fan_out)— limit = √(6/(fan_in+fan_out)), samples Uniform(−limit, +limit)
Xavier initialisation (Glorot & Bengio, 2010) is the standard weight initialisation for Transformer linear layers. It keeps activation variance stable across layers by scaling the init range to the layer's input/output sizes.
Matrix* matrix_add(const Matrix *A, const Matrix *B);
Matrix* matrix_sub(const Matrix *A, const Matrix *B);
Matrix* matrix_elementwise_mul(const Matrix *A, const Matrix *B); /* Hadamard */
Matrix* matrix_scalar_mul(const Matrix *A, float scalar);
Matrix* matrix_scalar_add(const Matrix *A, float scalar);Architecture note: all three element-wise binary ops share a single static helper:
static Matrix* elementwise_binary(const Matrix *A, const Matrix *B, BinaryOp op);This eliminates triplicated allocation/guard code. The compiler inlines the trivial op_add, op_sub, op_mul callbacks. After the _inplace refactor, the allocating versions became:
Matrix* matrix_relu(const Matrix *A) {
Matrix *out = matrix_copy(A);
if (!out) return NULL;
matrix_relu_inplace(out);
return out;
}The kernel logic lives in exactly one place — the _inplace function — so correctness improvements automatically apply to both calling patterns.
Matrix* matrix_matmul(const Matrix *A, const Matrix *B);
// A: (m×k) B: (k×n) → C: (m×n)Algorithm: classical triple-loop with register-local accumulator.
for (size_t i = 0; i < m; ++i) {
for (size_t j = 0; j < n; ++j) {
float acc = 0.0f; /* stays in register across inner loop */
for (size_t p = 0; p < k; ++p)
acc += A->data[i*k + p] * B->data[p*n + j];
C->data[i*n + j] = acc;
}
}The acc variable prevents the compiler from generating load/store to C->data on every multiply-add. This is the most important scalar optimisation for matmul.
Access pattern analysis (documented in comments):
| Array | Access as p advances | Cache behaviour |
|---|---|---|
A->data[i*k + p] |
Sequential stride-1 | ✅ cache friendly |
B->data[p*n + j] |
Strided (stride = n) | |
C->data[i*n + j] |
Single write per (i,j) | ✅ optimal |
Three FUTURE: hooks are explicitly marked at the right loop levels:
FUTURE (tiling): split i/j/p into TILE_SIZE blocks to make B-tiles fit in L1FUTURE (SIMD): replace inner p-loop with_mm256_fmadd_psover 8 floatsFUTURE (OpenMP):#pragma omp parallel forover the i-loop
Matrix* matrix_softmax(const Matrix *A);
void matrix_softmax_inplace(Matrix *A);Algorithm: three-pass row-wise log-sum-exp stabilisation.
Pass 1: max_val = max_j A[i][j]
Pass 2: e[j] = exp(A[i][j] - max_val), sum = Σ e[j]
Pass 3: out[j] = e[j] / sum
Subtracting max_val before expf prevents IEEE 754 overflow (expf(800) → ∞) which would corrupt the entire row. This is the standard implementation used in PyTorch, JAX, and every production inference engine.
Direct Transformer relevance:
Attention(Q,K,V) = softmax(Q·Kᵀ / √d_k) · V
The Q·Kᵀ scores can be large (proportional to d_k). Without the max subtraction, large models would produce NaN attention weights on the very first forward pass.
typedef struct {
size_t rows;
size_t cols;
size_t stride; /* = parent->cols, NOT view->cols */
float *data; /* non-owning pointer into parent */
} MatrixView;
MatrixView matrix_view(const Matrix *parent,
size_t row_start, size_t col_start,
size_t rows, size_t cols);
Matrix* matrix_from_view(const MatrixView *view);
Matrix* matrix_slice(const Matrix *A, size_t r0, size_t r1,
size_t c0, size_t c1);The stride field is the critical insight. Given a parent of shape (M × N) and a view starting at column c0:
view.data = parent->data + row_start * N + c0
view.stride = N ← parent's full width
view.cols = c1 - c0 ← only the view's width
To read view[r][c]: view.data[r * view.stride + c]
The stride skips over the columns that belong to other views without any copy. This is how numpy views and torch.Tensor slices work internally.
Primary Transformer use case — Q/K/V projection split:
/* Project once: X·W_QKV shape (seq_len × 3*d_model) */
Matrix *proj = matrix_matmul(X, W_QKV);
/* Zero-copy views — no allocation */
MatrixView vQ = matrix_view(proj, 0, 0, seq_len, d_model);
MatrixView vK = matrix_view(proj, 0, d_model, seq_len, d_model);
MatrixView vV = matrix_view(proj, 0, 2*d_model,seq_len, d_model);
/* Materialise only when a contiguous buffer is required */
Matrix *Q = matrix_from_view(&vQ);matrix_from_view copies row-by-row using memcpy, which is safe because each source row is view.cols wide starting at offset i * view.stride.
Matrix* matrix_add_rowwise(const Matrix *A, const Matrix *bias); /* bias: 1×n */
Matrix* matrix_mul_rowwise(const Matrix *A, const Matrix *scale); /* scale: 1×n */
Matrix* matrix_add_colwise(const Matrix *A, const Matrix *bias); /* bias: m×1 */Why these three specifically?
Every normalisation and projection layer in a Transformer reduces to one of:
| Operation | Formula | Function |
|---|---|---|
| Linear bias add | y = xW + b |
matrix_add_rowwise(y, b) |
| LayerNorm scale | y * γ |
matrix_mul_rowwise(y, gamma) |
| LayerNorm shift | y + β |
matrix_add_rowwise(y, beta) |
| Per-token scaling | y + pos_bias |
matrix_add_colwise(y, pos) |
Shape contract is strictly enforced:
_rowwiserequiresbias->rows == 1 && bias->cols == A->cols_colwiserequiresbias->cols == 1 && bias->rows == A->rows
This precision matters because silent shape broadcasting errors (numpy-style) are the most common source of incorrect gradient flows in custom autograd engines.
size_t matrix_argmax_row(const Matrix *A, size_t row);
size_t* matrix_argmax_rows(const Matrix *A); /* caller free()s result */Transformer inference loop:
/* Forward pass produces logit matrix: (1 × vocab_size) */
size_t next_token = matrix_argmax_row(logits, 0);
/* Or for a batched decode of an entire sequence: */
size_t *tokens = matrix_argmax_rows(output_logits);
/* tokens[i] = predicted token at position i */
free(tokens);matrix_argmax_rows returns a heap-allocated size_t[] — documented clearly so the caller knows to free() it. This avoids the ambiguity of returning a pointer to a static buffer (not thread-safe) or requiring the caller to pre-allocate.
Tie-breaking policy: first (leftmost) maximum — deterministic, consistent with numpy and PyTorch defaults.
void matrix_relu_inplace(Matrix *A);
void matrix_softmax_inplace(Matrix *A);
void matrix_scalar_mul_inplace(Matrix *A, float scalar);
void matrix_scalar_add_inplace(Matrix *A, float scalar);
void matrix_add_inplace(Matrix *dst, const Matrix *src);
void matrix_sub_inplace(Matrix *dst, const Matrix *src);
void matrix_elementwise_mul_inplace(Matrix *dst, const Matrix *src);The design pattern applied:
Before (duplicated logic):
matrix_relu → allocate + loop
(no inplace)
After (single source of truth):
matrix_relu_inplace → actual kernel loop
matrix_relu → matrix_copy → matrix_relu_inplace
This means:
- Fixing a bug in
relu_inplaceautomatically fixesrelu— no drift. - Adding SIMD to
relu_inplaceautomatically acceleratesrelu. - Training code that needs to avoid allocation can call
_inplacedirectly. - The existing API is 100% backwards-compatible — nothing broke.
When to use which:
| Scenario | Use |
|---|---|
| Forward pass, clarity matters | matrix_relu(A) — allocates, pure |
| Training inner loop, perf matters | matrix_relu_inplace(buf) — zero alloc |
Residual connection: x = x + sublayer(x) |
matrix_add_inplace(x, sublayer_out) |
| Layer Norm γ scaling | matrix_elementwise_mul_inplace(x, gamma) |
typedef struct {
Matrix *data; /* forward-pass weights */
Matrix *grad; /* accumulated gradient — NULL until first use */
} Parameter;
Parameter* parameter_create(size_t rows, size_t cols);
Parameter* parameter_wrap(Matrix *data); /* takes ownership */
void parameter_free(Parameter *p);
int parameter_zero_grad(Parameter *p); /* lazy alloc on first call */
int parameter_accumulate_grad(Parameter *p, const Matrix *grad_delta);
void parameter_print(const Parameter *p, const char *label);Lazy gradient allocation is a deliberate design decision. In a Transformer:
- Embedding tables may have millions of parameters.
- During forward-only inference, gradients are never needed.
- Allocating
gradupfront would double peak memory for no benefit.
The pattern parameter_zero_grad → parameter_accumulate_grad × N → optimizer step maps exactly to PyTorch's:
optimizer.zero_grad()
loss.backward() # accumulates gradients
optimizer.step()Ownership semantics of parameter_wrap:
Matrix *W = matrix_random_xavier(...);
Parameter *p = parameter_wrap(W);
/* W is now owned by p — do NOT call matrix_free(W) separately */
parameter_free(p); /* frees both p and W */This is the same move-semantics pattern used by std::unique_ptr in C++ and Rust's ownership system — made explicit through documentation and naming.
Matrix* matrix_row_mean(const Matrix *A); /* (m×1) */
Matrix* matrix_row_variance(const Matrix *A, const Matrix *mean); /* (m×1) */Layer Normalisation formula:
LayerNorm(x) = (x − μ) / √(σ² + ε) × γ + β
Implementation in Part 3 using existing primitives:
Matrix *mu = matrix_row_mean(x);
Matrix *sigma2 = matrix_row_variance(x, mu); /* pass pre-computed mu */
Matrix *x_shift = matrix_sub_rowwise(x, mu); /* (x - μ) */
/* add ε to sigma2, take sqrt, divide x_shift ... */
Matrix *normed = matrix_mul_rowwise(x_hat, gamma);
Matrix *out = matrix_add_rowwise(normed, beta);The optional mean parameter in matrix_row_variance avoids recomputing it when the caller already has it — a concrete micro-optimisation that matters during training when LN runs at every sub-layer.
Total: 52 tests across 11 groups
| Group | Tests | What is validated |
|---|---|---|
| Creation | 5 | dimensions, zero-init, invalid args, identity, fill |
| Element access | 2 | get/set correctness, OOB safety |
| Addition | 2 | known values, shape mismatch |
| Subtraction | 2 | known values, shape mismatch |
| Scalar mul | 2 | correctness, multiply-by-zero |
| Elementwise mul | 1 | Hadamard product |
| Transpose | 3 | shape, values, symmetric matrix |
| Matrix mul | 4 | 2×3 × 3×2 hand-verified, identity, mismatch, 1×1 |
| Softmax | 2 | row sums = 1.0, uniform input → equal probs |
| ReLU | 1 | negative/zero/positive elements |
| Row statistics | 2 | mean correctness, variance formula |
| Copy | 1 | deep copy independence |
| Random init | 3 | range bounds, deterministic seeds, Xavier range |
| NULL safety | 1 | 12 NULL inputs, no crash, no UB |
| matrix_mean_all | 4 | known values, MSE pattern, single element, NULL |
| MatrixView / slice | 4 | stride, Q/K/V split, exact values, invalid bounds |
| Broadcasting | 4 | rowwise add, rowwise mul, colwise add, mismatch |
| Argmax | 3 | row search, all rows, single element |
| In-place ops | 5 | relu, softmax, scalar mul, add, inplace==allocating |
| Parameter | 5 | create, lazy grad, accumulate, wrap, shape mismatch |
Key test: test_inplace_vs_allocating_relu
This test is not just a correctness check — it's a contract test. It verifies that the refactor did not introduce any behavioural divergence between the allocating and inplace paths. Any future change to the kernel that breaks this equivalence will be caught immediately.
Key test: test_view_materialise
Simulates the exact Q/K/V split that will appear in Part 4:
Matrix *X = ...; /* shape (2 × 6), d_k = 2 */
MatrixView vQ = matrix_view(X, 0, 0, 2, 2);
MatrixView vK = matrix_view(X, 0, 2, 2, 2);
MatrixView vV = matrix_view(X, 0, 4, 2, 2);
Matrix *Q = matrix_from_view(&vQ);
/* Verifies Q[1][1] == 7.0, V[1][1] == 11.0, etc. */Decision: Matrix holds one float* buffer — no array-of-pointers.
Rationale: Array-of-pointers (float **) requires N+1 allocations, destroys cache locality, and complicates memcpy-based operations. A single flat buffer is optimal for everything from SIMD to serialisation.
Decision: Use calloc for all matrix data allocations.
Rationale: calloc requests zero pages from the OS, which are often already zeroed (recycled from other processes) — potentially free. memset touches all pages unconditionally. Zero-initialisation by default also eliminates an entire class of bugs where operations on newly created matrices produce garbage.
Decision: Allocating variants are wrappers over _inplace, not vice versa.
Rationale: If _inplace were a wrapper over the allocating version (copy, then copy back), it would (a) allocate unnecessarily and (b) require the caller to not be aliased. The chosen direction gives zero overhead for the performance path and no logic duplication.
Decision: p->grad is NULL at creation; allocated on first parameter_zero_grad or parameter_accumulate_grad call.
Rationale: Inference-only models never need gradients. In large models (embedding tables), pre-allocating grad would double memory. Lazy allocation is safe because parameter_zero_grad guarantees the buffer exists before any accumulation.
Decision: matrix_view() returns a MatrixView struct on the stack.
Rationale: The struct is 32 bytes (4 × 8 bytes on 64-bit). Stack allocation is free, avoids heap fragmentation, and makes the non-owning semantics obvious — you can't matrix_free() a struct on the stack. PyTorch's TensorView and NumPy's ndarray slice views follow the same principle.
Decision: Returns size_t*, caller calls free().
Rationale: Static buffers are not thread-safe. Stack buffers of variable length (VLAs) are a C99 pitfall. Heap allocation with documented ownership is the safest and most flexible contract, consistent with standard C string functions like strdup.
Decision: The summation loop uses a double accumulator; the result is cast back to float.
Rationale: Summing millions of small float values with float arithmetic introduces catastrophic cancellation (the Kahan summation problem). A double accumulator costs one extra register and avoids meaningful precision loss across batch sizes up to 2²⁶ elements — well beyond any realistic loss computation. This matches the strategy used by ATLAS cblas_sdot and Python's math.fsum.
Before closing Part 0, all three readiness checkpoints were formally verified against the live API.
Matrix *y = matrix_matmul(x, W); // x:(batch×in) W:(in×out) → y:(batch×out)
Matrix *out = matrix_add_rowwise(y, b); // b:(1×out) broadcast to every row| Requirement | Function | Status |
|---|---|---|
| Matrix multiplication | matrix_matmul |
✅ dimension-checked |
| Bias broadcast | matrix_add_rowwise |
✅ shape contract enforced |
Result: ✅ Linear layer implementable with zero new code.
Parameter *W_param = parameter_create(in, out); // W_param->data = W
// W_param->grad = dW (lazy)
Parameter *b_param = parameter_create(1, out); // b_param->data = b
// b_param->grad = db (lazy)| Operation | Function | Status |
|---|---|---|
| Allocate W + dW together | parameter_create |
✅ |
| Zero gradients before batch | parameter_zero_grad |
✅ lazy alloc on first call |
| Accumulate dW | parameter_accumulate_grad |
✅ shape-validated |
| Free both matrices | parameter_free |
✅ no double-free risk |
Result: ✅ Training foundation ready.
MSE = mean( (y_pred − y_true)² )
Matrix *diff = matrix_sub(y_pred, y_true); // element-wise difference
Matrix *sq = matrix_elementwise_mul(diff, diff); // element-wise square
float loss = matrix_mean_all(sq); // global scalar reduction| Step | Function | Status |
|---|---|---|
y_pred − y_true |
matrix_sub |
✅ |
diff² (Hadamard) |
matrix_elementwise_mul |
✅ |
| Global mean | matrix_mean_all |
✅ added as final Part 0 item |
Result: ✅ MSE loss computable in 3 lines. Part 1 can start immediately.
matrix_mean_allwas the one function identified as missing during the audit. It was added, documented, and covered by 4 tests (including the exact MSE pattern above) before Part 0 was closed.
Every performance insertion point is marked in source with FUTURE: comments.
matrix_matmul inner p-loop → AVX2 _mm256_fmadd_ps (8× float throughput)
matrix_matmul i/j/p loops → cache tiling (TILE_SIZE = 32–64)
matrix_matmul i-loop → #pragma omp parallel for
matrix_relu element loop → _mm256_max_ps(x, zero)
matrix_scalar_mul element loop → _mm256_mul_ps
elementwise_binary loop → _mm256_add/sub/mul_ps
Tiling sketch for matmul (the most impactful optimisation):
#define TILE 32
for (size_t ii = 0; ii < m; ii += TILE)
for (size_t jj = 0; jj < n; jj += TILE)
for (size_t pp = 0; pp < k; pp += TILE)
/* 32×32 tile of B fits in ~4 KB — stays in L1 cache */
for (i in [ii, ii+TILE]) for (j in [jj, jj+TILE]) for (p in [pp, pp+TILE])
acc += A[i][p] * B[p][j];This shows exactly which mini_tensor functions each future Part will call.
Part 1 — Linear Layer y = xW + b
matrix_matmul(x, W)
matrix_add_rowwise(y, b)
parameter_create / parameter_accumulate_grad
Part 2 — Embeddings lookup: token_id → vector
matrix_get / matrix_set (embedding table lookup)
matrix_create(vocab_size, d_model)
Part 3 — Layer Normalisation (x−μ)/√(σ²+ε) × γ + β
matrix_row_mean
matrix_row_variance
matrix_sub_rowwise [to be added]
matrix_scalar_add_inplace (add ε)
matrix_mul_rowwise(x_hat, γ)
matrix_add_rowwise(scaled, β)
Part 4 — Multi-Head Attention softmax(QKᵀ/√d_k) V
matrix_view / matrix_from_view (Q,K,V split)
matrix_matmul(Q, K_T) (scores)
matrix_scalar_mul_inplace(scores, 1/√d_k)
matrix_softmax_inplace(scores) (attention weights)
matrix_matmul(weights, V) (context)
Part 5 — FFN sub-layer ReLU(xW₁+b₁)W₂+b₂
matrix_matmul × 2
matrix_add_rowwise × 2
matrix_relu_inplace
Part 6 — Positional Encoding x + PE
matrix_add_inplace(x, pe)
Part 7 — Transformer Block compose Parts 3–6
All of the above + matrix_add_inplace (residual connections)
Part 8 — Backpropagation dL/dW, dL/db, ...
parameter_zero_grad
parameter_accumulate_grad
matrix_elementwise_mul (gradient of ReLU)
matrix_matmul (gradient of Linear: dL/dX = dL/dY · Wᵀ)
matrix_transpose (gradient of Linear: dL/dW = Xᵀ · dL/dY)
Part 9–15 — Adam, BPE, Encoder/Decoder, Training Loop
All prior + matrix_scalar_mul_inplace (weight update: W -= lr * grad)
matrix_elementwise_mul_inplace (Adam momentum update)
- A correct, safe, zero-dependency foundation for a full Transformer implementation in pure C
- A learning-oriented codebase where every decision is documented and explained
- A forward-compatible architecture — every API is designed to grow without breaking changes
- A reference implementation — readable enough that you can verify correctness against the paper's equations by inspection
- It is not a production inference engine (no BLAS, no SIMD)
- It is not a training framework (no optimizer, no computation graph)
- It does not handle batching across the leading dimension (all ops are 2D — batching is via explicit loops)
These are intentional choices for Part 0. The architecture does not prevent any of them from being added later.
mini_tensor is a 3,771-line, zero-dependency, 56-test pure-C matrix library that forms the mathematical bedrock for a complete Transformer implementation from scratch.
Every major design decision — row-major layout, calloc zero-init, _inplace refactor pattern, lazy gradient allocation, MatrixView stride semantics, double accumulation in matrix_mean_all — was made with the full 15-part roadmap in mind.
The library compiles cleanly under the strictest GCC flags (-Wall -Wextra -Werror -pedantic) and passes all 56 tests including valgrind-verified memory safety. A formal three-checkpoint audit confirmed that the linear layer forward pass, parameter storage for training, and MSE loss computation are all implementable without adding any further API surface. Part 1: Linear Layer is unblocked.