Metadata-Version: 2.2
Name: mtl5
Version: 5.9.8
Summary: Python bindings for MTL5 — NumPy/SciPy/JAX/PyTorch interop with hardware accelerator dispatch
License: MIT
Requires-Python: >=3.10
Requires-Dist: numpy>=1.24
Provides-Extra: scipy
Requires-Dist: scipy>=1.10; extra == "scipy"
Provides-Extra: torch
Requires-Dist: torch>=2.0; extra == "torch"
Provides-Extra: jax
Requires-Dist: jax[cpu]>=0.4; extra == "jax"
Provides-Extra: pandas
Requires-Dist: pandas>=2.0; extra == "pandas"
Provides-Extra: sklearn
Requires-Dist: scikit-learn>=1.3; extra == "sklearn"
Provides-Extra: notebooks
Requires-Dist: matplotlib>=3.7; extra == "notebooks"
Requires-Dist: jupyterlab>=4.0; extra == "notebooks"
Requires-Dist: scipy>=1.10; extra == "notebooks"
Provides-Extra: all
Requires-Dist: scipy>=1.10; extra == "all"
Requires-Dist: torch>=2.0; extra == "all"
Requires-Dist: jax[cpu]>=0.4; extra == "all"
Requires-Dist: pandas>=2.0; extra == "all"
Requires-Dist: scikit-learn>=1.3; extra == "all"
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: numpy>=1.24; extra == "dev"
Requires-Dist: scipy>=1.10; extra == "dev"
Description-Content-Type: text/markdown

# mtl5-python

Python bindings for [MTL5](https://github.com/stillwater-sc/mtl5) — NumPy/SciPy/JAX/PyTorch interop with hardware accelerator dispatch.

Built with [nanobind](https://github.com/wjakob/nanobind) for minimal overhead and zero-copy array interop.

## Install

```bash
pip install mtl5
```

Prebuilt wheels are published for Linux, macOS, and Windows on CPython
3.10–3.12, so no compiler is required. Optional ecosystem integrations are
available as extras:

```bash
pip install "mtl5[scipy]"   # SciPy sparse interop
pip install "mtl5[all]"     # scipy, torch, jax, pandas, scikit-learn
```

### From source

Building the extension yourself — for contributing, or to enable a
[build option](#build-options) — needs Python 3.10+ and a C++20 compiler
(GCC 12+, Clang 15+, MSVC 2022):

```bash
pip install .
```

See [Development](#development) for an editable install.

## Quick start

```python
import numpy as np
import mtl5

# Vectors and norms
v = mtl5.vector(np.array([3.0, 4.0]))
print(mtl5.norm(np.array([3.0, 4.0])))  # 5.0

# Dot product
a = np.array([1.0, 2.0, 3.0])
b = np.array([4.0, 5.0, 6.0])
print(mtl5.dot(a, b))  # 32.0

# Solve Ax = b
A = np.array([[2.0, 1.0], [1.0, 3.0]])
b = np.array([5.0, 7.0])
x = mtl5.solve(A, b)
print(x)  # [1.6, 1.8]
```

## Mixed precision

A mixed-precision operation has three independent precisions. MTL5 supplies the
precision-generic kernels, [Universal](https://github.com/stillwater-sc/universal)
supplies the number systems, and this package composes them:

| | chosen by |
|---|---|
| **element** (storage) | the container you pass in — `mtl5.convert(x, "posit16")` |
| **accumulator** (compute) | `accumulator=` |
| **result** (delivery) | `result=` |

```python
import numpy as np, mtl5

x = np.random.default_rng(0).standard_normal(4000)
v = mtl5.convert(x, "posit16")  # store narrow

mtl5.mixed.dot(v, v)  # accumulate in posit16 too
mtl5.mixed.dot(v, v, accumulator="f64")  # ...or in double
mtl5.mixed.dot(v, v, accumulator="quire")  # ...or exactly, in Universal's quire
```

Accumulating 4000 posit16 products, against the exact value of the same
posit16 data:

| accumulator | relative error |
|---|---|
| none (posit16) | 1.8 × 10⁻¹ |
| `"f32"` | 3.8 × 10⁻⁷ |
| `"f64"` | 1.1 × 10⁻¹⁶ |
| `"fma"` | 1.1 × 10⁻¹⁶ |
| `"quire"` | **0** — bit-exact |

`mtl5.mixed.accumulators(dtype)` lists what a given element type supports. The
quire is available for the posit, cfloat, lns and fixpnt families; `f32`/`f64`
have none (Universal defines no quire for the native types). Exactness varies by
family — it is genuinely exact for posit and fixpnt, while the cfloat and lns
quires have known upstream limitations documented in
`python/include/mtl/math/quire_accumulator.hpp`.

`accumulator=` is available on `dot`, `norm` (ord=2), `frobenius_norm`,
`matvec` and `matmul`.

### Iterative refinement

Factor cheaply in a low precision, then recover accuracy with a residual formed
in float64:

```python
x, info = mtl5.mixed.lu_iterative_refine(A, b, working="posit16", rel_tol=1e-14)
# info -> {'iters': 7, 'rel_residual': ..., 'converged': True}
mtl5.mixed.backward_error(A, x, b)
```

On a 100×100 system, forward error of the refined solution:

| working precision | iterations | forward error |
|---|---|---|
| `fp16` | 5 | 1.9 × 10⁻⁸ |
| `posit16` | 7 | 1.9 × 10⁻¹⁵ |
| `f32` | 2 | 3.8 × 10⁻¹⁶ |
| `f64` | 0 | 3.1 × 10⁻¹⁶ |

The result is always the *best* iterate found, so an over-long `max_iter` never
degrades the answer.

`mtl5.mixed.iterative_refine(A, M, b)` is the sparse counterpart, refining
through any factorization exposing `solve()` — the sparse direct factorizations
below, or the ILU(0)/IC(0) preconditioners.

## Matrix Market I/O and sparsity pictures

```python
A = mtl5.io.mm_read("circuit.mtx")  # -> CSR
D = mtl5.io.mm_read_dense("small.mtx")  # -> dense
mtl5.io.mm_write_sparse("out.mtx", A, "comment")

mtl5.io.spy(A, "pattern.png")  # binary non-zero pattern
mtl5.io.spy_magnitude(A, "mag.png", log_scale=True)
mtl5.io.spy_density(A, "dens.png", max_pixels=512)
```

Two things differ from what a SciPy user expects.

**The function picks the container, not the file.** `scipy.io.mmread` returns an
ndarray for an `array` file and COO for a `coordinate` one. Here `mm_read`
always gives CSR and `mm_read_dense` always gives dense, and *both* accept
either file format — so reading a dense `.mtx` into CSR is a deliberate call
rather than an error.

**The PNGs are uncompressed.** MTL5 writes them with a from-scratch encoder
emitting DEFLATE *stored* blocks, which is what lets `spy` work with no image
library and no plotting stack in the process. File size is therefore about
`width × height × channels` — measured at the default `max_pixels=1024`, that
is ~1.0 MB for `spy` (grayscale) and ~3.1 MB for `spy_magnitude` /
`spy_density` (RGB). Pipe the output through any PNG optimizer if that matters.

`.gz` inputs are read transparently only when MTL5 is built with zlib
(`-C cmake.define.MTL5_WITH_ZLIB=ON`); `mtl5.build_info()["zlib"]` reports it,
and without it a `.gz` path raises rather than silently mis-parsing.

## Test matrices

A named catalog of matrices with known pathologies — the inputs a
mixed-precision experiment actually wants, and tedious to hand-roll correctly:

```python
g = mtl5.generators

g.hilbert(8)  # cond ~1.5e10 — the canonical ill-conditioned matrix
g.clement(7)  # eigenvalues exactly -6,-4,-2,0,2,4,6
g.randspd(4, [1, 10, 100, 1000])  # SPD with exactly that spectrum
g.randsym(4, [-5, -1, 2, 7])  # controlled *indefinite* matrix
g.randsvd(20, 20, kappa=1e6)  # condition number exactly 1e6
g.laplacian_2d(64, 64)  # sparse, the usual solver benchmark
```

Also `frank`, `pascal`, `wilkinson`, `rosser`, `magic`, `lehmer`, `lotkin`,
`minij`, `ones`, `forsythe`, `kahan`, `moler`, `companion`, `vandermonde`,
`randorth`, `laplacian_1d`, `poisson2d`, and the published catalog via
`testsuite_names()` / `testsuite_kappa(name)`.

Every dense generator takes `dtype=`, which is what makes them useful here:

```python
H = g.hilbert(8, dtype="posit16")  # correctly rounded posit16 Hilbert
mtl5.cholesky(H)  # ...then watch it fail
```

Generation happens in float64 and `dtype=` rounds. That is the right semantics
for a test matrix — the definitions are over the reals, so you want the
correctly rounded representation of the exact entry, not the result of
evaluating the formula in low-precision arithmetic.

Range vectors follow NumPy, and take `dtype=` too:

```python
mtl5.linspace(0, 1, 5)
mtl5.arange(0, 10, 3)
mtl5.logspace(0, 3, 4)
mtl5.geomspace(1, 1000, 4)  # 1, 10, 100, 1000 — a true geometric
# progression, not logspace's exponents
```

## Dense factorizations

```python
qr = mtl5.qr(A)  # Householder QR; tall or square
x = qr.solve(b)  # least squares
qr.Q, qr.R

lq = mtl5.lq(A)  # the row-space counterpart
lq.L, lq.Q

ld = mtl5.ldlt(A)  # symmetric, indefinite allowed
ld.solve(b)
ld.diagonal()  # D — its signs are the inertia
```

All accept an MTL5 matrix or a float32/float64 NumPy array, alongside the
existing `mtl5.lu` and `mtl5.cholesky`.

**Every** dense factorization — `lu`, `qr`, `lq`, `cholesky`, `ldlt` and
`bunch_kaufman` — is available for float32, float64 and all fifteen Universal
dtypes, which is what makes a like-for-like comparison across number systems
possible. The integer element types are not supported.

Availability is not usability, though, and the two diverge sharply at narrow
widths. The orthogonalizing pair (`qr`, `lq`) fails first: below 16 bits every
Householder reflector rounds to zero, so `Q` comes back as the exact identity
and the factorization looks perfect while reconstructing nothing. The solving
three (`cholesky`, `ldlt`, `bunch_kaufman`) have no orthogonality to lose and
stay usable to 8 bits. `tests/test_universal_factorizations.py` records where
each format stops.

### Cholesky vs LDLᵀ across number systems

`ldlt` and `cholesky` are both available in every bound precision, which is what
makes the interesting comparison possible. Cholesky takes square roots, so
it refuses a matrix that has drifted out of positive-definiteness — the failure
mode of a Kalman covariance update in low precision. LDLᵀ has no square roots,
survives, and records what happened in `D`:

```python
P = np.eye(6)
P[3, 3] = -1e-3  # covariance went indefinite

mtl5.cholesky(mtl5.convert(P, "posit16"))  # RuntimeError: not SPD
d = mtl5.ldlt(mtl5.convert(P, "posit16")).diagonal()
(d < 0).any()  # True — D names the bad direction
```

`mtl5.bunch_kaufman(A)` is the pivoting variant, for a symmetric matrix that
plain `ldlt` rejects on a zero pivot.

## Eigenvalues, BLAS 2/3, and matrix properties

The eigen entry points mirror `numpy.linalg`, and return NumPy arrays:

```python
mtl5.eigvalsh(A)  # symmetric eigenvalues, real, ascending
w, Q = mtl5.eigh(A)  # ...with eigenvectors:  A = Q diag(w) Qᵀ
mtl5.eigvals(A)  # general spectrum, complex
w, V = mtl5.eig(A)  # ...with right eigenvectors
mtl5.spectral_radius(A)
mtl5.inertia(A)  # {'positive': …, 'negative': …, 'zero': …}
```

BLAS levels 2 and 3 write into a caller-supplied output, as BLAS does — that
in-place accumulation is the point:

```python
mtl5.ger(alpha, x, y, A)  # A += alpha x yᵀ
mtl5.symv(alpha, A, x, beta, y)  # y = alpha A x + beta y
mtl5.trsv(A, x, upper=True)  # x = A⁻¹x
mtl5.trmm(alpha, A, B, upper=True)  # B = alpha A B
mtl5.trsm(alpha, A, B, upper=True)  # solve A X = alpha B
mtl5.symm(alpha, A, B, beta, C)  # C = alpha A B + beta C
mtl5.syrk(alpha, A, beta, C)  # C = alpha A Aᵀ + beta C
mtl5.syr2k(alpha, A, B, beta, C)
```

Property predicates come in two cost classes. The **docstrings are
authoritative**; the split is:

- **O(n²) or cheaper** — the structural checks (`is_square`, `is_empty`,
  `is_symmetric`, `is_hermitian`, `is_triangular` and the upper/lower variants,
  `is_diagonal`, `is_banded`, `is_diagonally_dominant`) and every vector
  predicate.
- **O(n³)** — anything that forms a product, factorizes, or runs an
  eigensolve: `is_orthogonal`, `is_unitary`, `is_normal`, `is_spd`,
  `is_positive_definite`, `is_singular`, `is_nonsingular`, `is_invertible`,
  `spectral_radius`, `inertia`, `is_indefinite`. Don't put these inside a loop.

SVD and the queries built on it:

```python
U, s, V = mtl5.svd(A)  # s is the vector of singular values (NumPy's convention)
mtl5.svdvals(A)  # singular values only — cheaper
mtl5.condition_number(A)  # σ_max / σ_min
mtl5.rcond(A)  # σ_min / σ_max, safer near singular
mtl5.numerical_rank(A)
mtl5.nullity(A)
```

`svd` takes a `tol`. `V` and the reconstruction are accurate to machine
precision regardless, but `U`'s orthogonality is bounded by the iteration's
tolerance rather than by eps. Over 160 matrices (n = 3..20, four tolerances)
`‖UᵀU − I‖` had a median of about 1.4×tol and a worst case of 4.5×tol; the
regression test asserts 20×tol to leave headroom. Treat `tol` as the knob:
tighten it if you need an orthonormal `U` specifically, rather than relying on
a particular multiple.

## Complex numbers

`mtl5.vector()` and `mtl5.matrix()` accept `complex64` and `complex128` arrays
and give the same zero-copy views as the real types, with dtypes `c64`/`c128`:

```python
A = mtl5.matrix(np.array([[2 + 0j, 1 - 1j], [1 + 1j, 3 + 0j]]))
b = mtl5.vector(np.array([3 + 1j, 1 + 4j]))
x = mtl5.solve(A, b)  # complex LU with partial pivoting
mtl5.norm(x, 2)  # a real float, not a complex
A.real.to_numpy()  # real part, as a real matrix
```

Three things differ from the real case, and getting them wrong is quiet rather
than loud, so they are worth stating.

**`dot` is Hermitian.** It computes `sum(conj(a[i]) * b[i])`, conjugating the
*first* argument — that is NumPy's `vdot`, not NumPy's `dot`. The unconjugated
bilinear product is `mtl5.dot_real`, which is what `np.dot` does for 1-D
complex. Both exist because both are wanted; the names say which is which.

**`.T` does not conjugate.** MTL5's transpose is the plain one. `.H` (or
`mtl5.adjoint`) is the conjugate transpose. For real elements the two coincide,
which is exactly why the distinction has to be explicit here.

**Hermitian and symmetric are different properties.** `mtl5.is_hermitian(A)`
tests `A == Aᴴ` and `mtl5.is_symmetric(A)` tests `A == Aᵀ`; for complex those
are not the same matrix, and which one you have decides which solver is right.

What is available: containers and factories, all four norms, `dot`/`dot_real`,
`matmul`/`matvec`, `solve`/`lu`/`inv`, `transpose`/`adjoint`/`conj`,
`ldlt_solve`, `cholesky`, and `qr`/`lq`. What is not: `bunch_kaufman`, the eigen
and SVD family, and the Krylov solvers — MTL5 has no complex implementation of
any of them, and complex input raises a `TypeError` naming the alternative
rather than silently taking a real part.

**Complex least squares** works through `qr`, which uses MTL5's complex
Householder:

```python
f = mtl5.qr(A)  # A complex, num_rows >= num_cols
x = f.solve(mtl5.vector(b))  # min ||Ax - b||_2
f.Q, f.R  # Q is unitary: Q^H Q = I
```

`lq` is the row-space counterpart, for the underdetermined case. Both were
checked for the answer rather than the compile: the least-squares residual is
orthogonal to range(A) to 3.5e-15, which is what distinguishes a solve applying
Qᴴ from one applying Qᵀ — the latter would reconstruct `A = QR` perfectly and
still solve the wrong problem.

**Symmetric and Hermitian need different factorizations**, and `mtl5.ldlt_solve`
picks from the matrix: LDLᵀ for a complex symmetric one, LDLᴴ for a Hermitian
one. They are not interchangeable — sending either input to the other's kernel
gives a wrong answer — so a matrix that is neither is refused rather than
guessed at.

```python
S = mtl5.matrix(np.array([[2 + 1j, 1 - 1j], [1 - 1j, 3 + 2j]]))  # A == A^T
mtl5.ldlt_solve(S, b)  # LDL^T

H = mtl5.matrix(np.array([[2 + 0j, 1 - 1j], [1 + 1j, 3 + 0j]]))  # A == A^H
mtl5.ldlt_solve(H, b)  # LDL^H
```

`mtl5.cholesky` follows the same principle. For complex it computes A = L·Lᴴ,
since the plain L·Lᵀ form and its diagonal ordering test mean nothing for
complex elements — MTL5 `static_assert`s against that and this routes to its
Hermitian variant. A non-real diagonal is reported as *not Hermitian* rather
than as a definiteness failure, since those are different mistakes.

Complex is not in `mtl5.dtypes()`, which lists what `mtl5.convert()` accepts —
the Universal number systems are real-only, so there is no complex target.

## N-dimensional arrays

`mtl5.array` exposes MTL5's `mtl/array` layer — ranks 1 through 4, float32 and
float64:

```python
import mtl5

x = mtl5.array.asarray(a)  # zero-copy view of any strided NumPy array
x[1, 2]  # element access
x[:, 1:3]  # a view, never a copy
x.T.sum_axis(0)  # transpose is a view too
x.to_numpy()  # back to NumPy, strides included
```

Rank is a C++ template parameter, so the ranks are fixed at build time. Rank 5+
raises rather than silently flattening.

The point of the layer here is that it makes the dense containers sliceable
without a copy:

```python
M = mtl5.matrix(a)
col = mtl5.array.as_ndarray(M)[:, 2]  # a column of a DenseMatrix, no copy
```

`asarray` accepts any strided layout, so a NumPy transpose or slice comes
through without materialising. Three things raise rather than quietly doing
something else, all because the result aliases your buffer:

* **negative strides** — MTL5's strides are unsigned
* **any dtype but float32/float64** — converting would return a view of a
  temporary, which is neither zero-copy nor the dtype you asked for
* **read-only arrays** — the view can write through

Two details worth knowing:

**`.strides` is in elements, not bytes.** That is how MTL5 holds them. NumPy's
`.strides` is in bytes, so the two differ by itemsize — `[3, 1]` here is
`(24, 8)` there for float64.

**Broadcasting is same-rank only.** An extent of 1 stretches, so `(2,3) + (2,1)`
works, but NumPy's rank promotion — `(2,3) + (3,)` — does not, because MTL5's
`broadcast_shape` takes two shapes of equal rank. Reshape the operand first.

`reshape` and `ravel` return a view when the layout allows it and a copy
otherwise, and never error; `flatten` always copies. All three give NumPy's
element order even for a transposed or sliced source. MTL5's own `reshape`
throws rather than copying in that case, which is right for a C++ caller but is
not NumPy's contract — see `docs/gap-analysis-2026-08.md` §3.9.

## Accumulator policy on the sparse factorizations

The four sparse direct factorizations that use a dense numeric workspace —
`splu`, `klu`, `supernodal_lu`, `supernodal_ldlt` — take an `accumulator=`
argument that types that workspace. A float32 factor can accumulate its updates
in float64:

```python
f = ms.splu(A32, accumulator="f64")  # narrow factor, wide arithmetic
f.accumulator  # 'f64'
```

The factor itself stays in the element type. Only the arithmetic that produced
it widens: the accumulator removes the intermediate roundings inside each
column's update chain, and each L/U entry is still rounded once on the way out.
This is the mixed-precision knob a fixed-precision library structurally cannot
offer.

On a badly scaled random sparse matrix (n = 400, cond ≈ 2.8e10) factored in
float32, forward error of the direct solve — measured on x86-64 Linux:

| ordering | nnz(L+U) | `f32` | `f64` | gain |
|---|---|---|---|---|
| `colamd` | 81 712 | 8.19 × 10⁻³ | 6.39 × 10⁻³ | 1.28× |
| `amd` | 92 026 | 5.35 × 10⁻² | 3.17 × 10⁻² | 1.69× |
| `natural` | 98 729 | 3.52 × 10⁻² | 1.39 × 10⁻² | 2.54× |

The *direction* is robust — a wider accumulator gives a smaller forward error in
every ordering, matrix and platform tried, and the test suite asserts that. The
*magnitude* is not: it ranges from roughly 1.3× to 3.3×, and which ordering
benefits most changes between platforms, so read the column above as one
measurement rather than a trend. (An earlier draft of this section claimed the
gain grows with fill; that held on Linux and reversed on macOS.)

A 1.3–3× improvement is worth having but is not the headline. The headline is
what it does to iterative refinement, which is what you would actually pair a
narrow factor with — same matrix, `natural` ordering, refined against a float64
residual:

| accumulator | iterations | forward error |
|---|---|---|
| `f32` | 6 | 2.0 × 10⁻⁹ |
| `f64` | **4** | **1.4 × 10⁻¹⁰** |

A third fewer iterations and 14× closer, from a factor occupying the same
memory.

Valid accumulators are `"f32"`, `"f64"`, `"fma32"`, `"fma64"`, or `None` for the
element type. `fma64` is measurably identical to `f64` here. Two are refused
rather than accepted quietly:

* **narrower than the element type** — `f32` on a float64 matrix loses precision
  rather than gaining it
* **`quire`** — `SparseMatrix` is float32/float64 only, and Universal defines a
  quire only for its own number systems

One asymmetry worth knowing: `refactor()` replays the stored pivot sequence
through an upstream entry point that takes no accumulator, so on `splu`, `klu`
and `supernodal_lu` it **refuses** when a non-default accumulator is in effect
rather than silently reverting to element precision — construct a new
factorization instead. `supernodal_ldlt.refactor()` re-runs the full numeric
factorization, so it carries the policy through.

## Sparse storage formats

CSR is the default and what every solver here takes. Two other formats are
available when the shape of the work suits them:

```python
import mtl5.sparse as ms

C = ms.coo_matrix(row, col, data, shape=(n, n))  # triplets
C.insert(3, 7, 1.5)  # append; no pattern needed up front
A = C.tocsr()  # duplicates sum here

E = ms.ell_matrix(A)  # fixed-width, vectorisable
E.padding_ratio  # ...if this is small
```

**COO** is the format you build in when you don't know the pattern yet. It
accumulates: a repeated `(row, col)` sums rather than overwrites, exactly as
`scipy.sparse.coo_matrix` does. The duplicates stay separate — `nnz` counts them
individually and a round trip through `to_scipy` preserves them — until
`tocsr()` folds them into one entry.

**ELL** stores `nrows × max_width` slots and pads the short rows, which is what
makes it regular enough to vectorise over. The cost is that one long row forces
every other row to carry empty slots, so it is only a good idea when the row
widths are near-uniform. `padding_ratio` is that decision as a number — 0.0 when
every slot is occupied, and rising toward 1.0 as the matrix gets more ragged:

```python
# a 2-D Laplacian: every row has 5 entries bar the boundary
ms.ell_matrix(mtl5.generators.laplacian_2d(32, 32)).padding_ratio  # 0.025

# an arrow matrix: one dense row and column, the rest diagonal
ms.ell_matrix(ms.from_scipy(arrow)).padding_ratio  # 0.985 — use CSR
```

ELL has no incremental build path: MTL5's `ell_matrix` has no element setter, so
CSR is the only way in.

**There is no CSC.** MTL5's `compressed2D` used to accept a `col_major`
orientation tag and ignore it — the inserter, element access and `mult` all
treated the storage as row-major regardless, so genuine CSC arrays came back
transposed with no error. That is now a compile-time rejection upstream
([mtl5#355](https://github.com/stillwater-sc/mtl5/issues/355)), but there is
still no column-major container, so `from_scipy` converts a `csc_matrix` to CSR
at the boundary — which is honest about the cost.

The thing CSC is usually wanted for here is a transpose product, and that is
available directly, using MTL5's own transpose view rather than a scipy round
trip:

```python
A.rmatvec(x)  # A.T @ x, no second copy of the matrix
```

`mtl5.sparse.as_linear_operator(A)` uses it, so both directions of a SciPy
`LinearOperator` now stay inside MTL5.

## Tensor algebra

`mtl5.tensor` is index-notation tensor algebra — distinct from `mtl5.array`,
which is NumPy-shaped N-D data. This one has fixed dimensions, Einstein
summation and a metric:

```python
import numpy as np, mtl5.tensor as mtt

A = mtt.asarray(np.arange(9.0).reshape(3, 3))
x = mtt.asarray(np.array([1.0, 2.0, 3.0]))

mtt.contract(A, "ij", x, "j")  # A @ x
mtt.contract(A, "ji", x, "j")  # A.T @ x
mtt.outer(x, x)  # rank 2
```

Rank and dimension are both compile-time in MTL5, so each pair is a separate
type. **Ranks 1, 2 and 4 over dimensions 2, 3 and 4** cover the module: nothing
upstream uses rank 3, and rank 4 arises only as `outer(rank2, rank2)`. Anything
else raises rather than silently reshaping — use `mtl5.array` for general N-D
data.

**Index strings are read at runtime**, which MTL5's own `contract` cannot do —
it takes index names as compile-time template parameters. The space is small
enough to enumerate, though: the repeated index sits in one of two positions on
each side, so four instantiations cover every rank-2 contraction up to
relabelling. Python therefore reaches MTL5's real contraction rather than a
reimplementation. The letters are arbitrary; only which positions they share
matters, so `contract(A, "pq", B, "qr")` and `contract(A, "ij", B, "jk")` agree.

A trace (`"ii"`), two shared indices, or none at all are each refused with a
message naming which it was.

The metric operations cover raising and lowering. `raise` is a Python keyword,
so the rank-1 pair is `raise_index` / `lower_index`, alongside upstream's
`raise_first` / `lower_second` for rank 2:

```python
g = mtt.minkowski_metric()  # signature (-, +, +, +)
v = mtt.asarray(np.array([1.0, 2.0, 3.0, 4.0]))
mtt.lower_index(v, g)  # [-1, 2, 3, 4]
mtt.raise_index(mtt.lower_index(v, g), g)  # back to v
```

`is_symmetric` and `is_antisymmetric` take a tolerance, and
`SymmetricTensor_d{dim}_{dtype}` stores D(D+1)/2 components instead of D².

## Matrix views

Eight ways to read part of a matrix:

```python
mtl5.view.lower(A)  # and upper, strict_lower, strict_upper
mtl5.view.transposed(A)
mtl5.view.banded(A, 1, 1)  # tridiagonal
mtl5.view.map(A, [2, 0], [1, 2])  # out[i, j] = A[rows[i], cols[j]]
mtl5.view.hermitian(A)
```

These **materialise** — each returns an ordinary `DenseMatrix` that every other
binding accepts. Upstream they are lazy accessors holding a `const Matrix&`,
which buys nothing here (Python has no expression template to feed) and is a
lifetime hazard: a caller passing a NumPy array or a converted scipy matrix
hands over a temporary the view would outlive.

Two carry a trap worth knowing.

**`banded` takes bandwidths, not signed diagonal offsets.** `banded(A, 1, 1)` is
tridiagonal and `banded(A, 0, 0)` the diagonal. Writing `-1` for "one
subdiagonal" reads plausibly but asks upstream for the wrong region entirely —
negative values are refused rather than silently returning a band you did not
mean.

**`hermitian` is not the adjoint.** It reads the **upper triangle** and mirrors
it conjugated into the lower, discarding whatever was stored there — the usual
"this matrix is Hermitian, only one triangle is meaningful" convention. For Aᴴ
use `mtl5.adjoint(A)`. The two agree only when the input is already Hermitian,
which is exactly the case you would reach for first when checking:

```python
Z = np.array([[1 + 0j, 2 - 1j], [9 + 9j, 4 + 0j]])  # lower triangle inconsistent on purpose
mtl5.view.hermitian(mtl5.matrix(Z))  # [[1, 2-1j], [2+1j, 4]] — lower discarded
mtl5.adjoint(mtl5.matrix(Z))  # [[1, 9-9j], [2+1j, 4]] — a real transpose
```

A non-real diagonal is refused: the view leaves the diagonal alone, so it would
produce a matrix that is not in fact Hermitian.

## Multigrid and smoothers

```python
import mtl5, mtl5.sparse as ms

M = mtl5.mg.multigrid_1d(A, n_levels=5)  # builds the whole hierarchy
x = M.vcycle(mtl5.vector(x0), mtl5.vector(b), cycles=10)
M.level_sizes  # [127, 63, 31, 15, 7]
```

On 1-D Poisson the residual drops by a factor of about **0.03 per V-cycle**, and
that factor holds as the problem grows — mesh independence is what separates
multigrid from a lone smoother, and the test suite asserts it rather than just
checking the residual went down.

`multigrid_1d` builds the hierarchy for you: standard 1-D coarsening, Galerkin
coarse operators, the chosen smoother at every level. `n_levels` is an upper
bound — coarsening stops before a level would fall below 4 rows, so check
`level_sizes`. `wcycle` is also available.

The seven smoothers are reachable directly, for use outside a hierarchy:

```python
x = mtl5.mg.smooth(A, x, b, kind="symmetric_gauss_seidel", sweeps=5)
mtl5.mg.smoothers()
# ['jacobi', 'gauss_seidel', 'backward_gauss_seidel',
#  'symmetric_gauss_seidel', 'sor', 'backward_sor', 'symmetric_sor']
```

`omega` applies to the SOR variants; at `omega=1.0` SOR is exactly Gauss-Seidel,
which the tests use to confirm the parameter reaches the kernel.

**Jacobi is a poor multigrid smoother**, and that is a property of the method
rather than of this binding. Multigrid needs a smoother that damps
*high-frequency* error — the part the coarse grid cannot represent — and
undamped Jacobi leaves the highest modes almost untouched. The usual remedy is
damped Jacobi at `omega ≈ 2/3`, which MTL5's `jacobi` has no parameter for. Use
a Gauss-Seidel or SOR variant in a hierarchy.

Grid transfer is exposed too, along with a sparse Galerkin product:

```python
R = mtl5.mg.make_restriction_1d(127)  # 63 x 127, full weighting
P = mtl5.mg.make_prolongation_1d(63)  # 127 x 63, linear interpolation
Ac = mtl5.mg.galerkin(R, A, P)  # the coarse operator, still sparse
```

`galerkin` exists because spelling it `R @ A @ P` would go through a
sparse-times-sparse product that returns a **dense** matrix — a fine-sized
intermediate, which is exactly what you cannot afford in a hierarchy.

## Krylov solvers and preconditioners

Ten solvers and eight preconditioners, in any combination:

```python
import mtl5.sparse as ms

x, info = ms.gmres(A, b, M=ms.ilu0(A))
x, info = ms.cg(A, b, M=ms.ic0(A), rtol=1e-12)
x, info = ms.iterative_solve(A, b, solver="idr_s", M=ms.ssor(A, omega=1.4), s=8)
```

`ms.solvers()` and `ms.preconditioners()` list them:

| | |
|---|---|
| **symmetric** | `cg`, `minres` |
| **general** | `bicgstab`, `bicgstab_ell`, `cgs`, `gmres`, `idr_s`, `tfqmr`, `bicg`, `qmr` |
| **preconditioners** | `identity`, `diagonal`, `ic0`, `ildl`, `ilu0`, `ssor`, `ilut`, `block_diagonal` |

`M=` defaults to `identity`, so omitting it runs unpreconditioned. Solver-specific
knobs are `restart=` (gmres), `ell=` (bicgstab_ell) and `s=` (idr_s).

The preconditioner is type-erased rather than being a template parameter of each
solver, so the pairing is chosen at runtime. Binding the cross product directly
would have meant 10 × 8 × 2 dtypes = **160** instantiations of a full iterative
solver; this is 20.

**`bicg` and `qmr` require a symmetric preconditioner.** They are the only two
that apply Mᵀ, and MTL5 implements a preconditioner's adjoint as its forward
solve — exact when M is symmetric, wrong otherwise. `identity`, `diagonal`,
`ic0` and `ildl` are always symmetric; `ilu0`, `ssor`, `ilut` and
`block_diagonal` are symmetric only when A is, which each one determines at
construction and reports as `.is_symmetric`. A pairing that would break down is
refused rather than run:

```python
ms.bicg(nonsymmetric_A, b, M=ms.ilu0(nonsymmetric_A))  # ValueError, by design
ms.bicg(nonsymmetric_A, b, M=ms.diagonal(nonsymmetric_A))  # fine
```

One behaviour worth knowing: **TFQMR can stagnate** unpreconditioned — it
plateaus rather than converging slowly, and more iterations do not help. Giving
it any real preconditioner fixes it. That is a property of the method, not of
this binding.

## Sparse direct solvers

Seven factorizations, one interface — construct, `.solve(b)`, `.refactor(A2)`:

| | for | notes |
|---|---|---|
| `ms.splu` | general square | Gilbert–Peierls, threshold pivoting |
| `ms.klu` | circuit matrices | block triangular form + per-block LU |
| `ms.supernodal_lu` | general square | dense block updates; `.nsuper` |
| `ms.cholesky` | symmetric positive definite | cheapest when it applies |
| `ms.ldlt` | symmetric, possibly indefinite | `.diagonal()` gives the inertia |
| `ms.supernodal_ldlt` | symmetric | dense block updates |
| `ms.qr` | least squares, rectangular | `min ‖Ax − b‖₂` |

```python
import mtl5.sparse as ms

lu = ms.splu(A, ordering="amd")  # analyze (ordering + symbolic) then factor
x = lu.solve(b)

lu.refactor(A2)  # numeric only — same pattern, new values
x2 = lu.solve(b2)

k = ms.klu(A)  # block triangular form + per-block LU
k.nblocks  # how reducible the matrix turned out to be
```

Two things `scipy.sparse.linalg.splu` cannot do.

**Refactorization.** A sequence of matrices sharing one sparsity pattern — the
circuit-transient case — pays for the ordering and symbolic analysis once. On a
2-D Laplacian, n=3600, nnz(A)=17,760:

| | nnz(factor) | factor | refactor | speedup |
|---|---|---|---|---|
| `splu` | 205,636 | 12.5 ms | 4.5 ms | 2.8× |
| `klu` | 119,530 | 6.7 ms | 2.2 ms | 3.0× |
| `supernodal_lu` | 205,636 | 19.6 ms | 4.0 ms | **5.0×** |
| `cholesky` | 59,765 | 22.3 ms | 20.8 ms | 1.1× |
| `ldlt` | 56,165 | 8.4 ms | 7.0 ms | 1.2× |
| `supernodal_ldlt` | 56,165 | 4.7 ms | 2.2 ms | 2.1× |

Two things to read off that table. Exploiting symmetry cuts the fill to about a
quarter, and `supernodal_ldlt` is both the sparsest and the fastest option for a
symmetric matrix. And the refactor win is small for `cholesky`/`ldlt` — they do
not pivot, so their analysis is a cheap symbolic pass and there is little to
skip; the LU-family factorizations, which must otherwise redo ordering *and* the
pivot search, gain the most.

The saving is the analysis, so it scales with how much of the runtime that is.
On a random pattern with catastrophic fill, numeric work dominates and refactor
is no faster (we measured a slight loss). Structured sparsity is where it pays.

**A factor narrower than the residual.** The factor's precision is chosen by
the dtype of `A`, independent of the precision you refine in:

```python
lu32 = ms.splu(A.astype(np.float32), ordering="amd")
x, info = mtl5.mixed.iterative_refine(ms.from_scipy(A), lu32, b, rel_tol=1e-14)
```

On a 2-D Laplacian, n=1600:

| | forward error |
|---|---|
| float64 factor, direct solve | 2.0 × 10⁻¹⁵ |
| float32 factor, direct solve | 8.9 × 10⁻⁷ |
| **float32 factor + float64 refinement** (3 iters) | **7.0 × 10⁻¹⁶** |

Half the factorization memory and traffic, and the refined answer is *better*
than the float64 direct solve. Every square factorization above works this way;
`qr` is excluded because least squares is not the square system refinement
corrects.

### Orderings

`ms.orderings()` lists `amd`, `colamd`, `rcm` and `natural`. Each is also
available standalone as a permutation, for inspection or for use on a scipy
matrix directly:

```python
p = ms.amd(A)  # or ms.colamd(A), ms.rcm(A), ms.ordering(A, name)
A[p][:, p]
```

The choice matters. Nonzeros in L+U for the 2-D Laplacian above (n=1600,
7840 nonzeros in A):

| ordering | `amd` | `colamd` | `rcm` | `natural` |
|---|---|---|---|---|
| nnz(L+U) | 41,542 | 62,944 | 90,040 | 128,078 |

`colamd` is the default because it suits unsymmetric matrices; `amd` is the
better choice when the pattern is symmetric, as here.

## Performance

### Threading

MTL5's kernels are parallel but **serial by default**. The thread pool is sized
once, on the first MTL5 call, and cannot be resized afterwards — so configure it
before you do any work:

```python
import mtl5

mtl5.set_num_threads(8)  # must precede the first MTL5 operation
print(mtl5.get_num_threads())
```

Equivalently, set `MTL5_NUM_THREADS` in the environment before importing:

```bash
MTL5_NUM_THREADS=8 python my_script.py
```

`set_num_threads()` raises `RuntimeError` rather than silently doing nothing if
the pool has already been sized. The count is clamped to hardware concurrency.

All non-trivial kernels release the GIL, so MTL5 work overlaps with other Python
threads instead of blocking the interpreter.

Double-precision `matmul`, 1000×1000, on a 20-core box:

| threads | 1 | 2 | 4 | 8 |
|---|---|---|---|---|
| GFLOP/s | 14.4 | 27.2 | 48.0 | 72.7 |

### Build options

The wheel builds MTL5's blocked GEMM vectorised through
[Google Highway](https://github.com/google/highway) by default. Everything else
is opt-in, since it needs a library on the build machine or makes the binary
non-portable:

```bash
pip install . -C cmake.define.MTL5_WITH_BLAS=ON      # external BLAS
pip install . -C cmake.define.MTL5_WITH_LAPACK=ON    # geev/syev/gesdd dispatch
pip install . -C cmake.define.MTL5_NATIVE_ARCH=ON    # -march=native, local builds only
```

`MTL5_NATIVE_ARCH=ON` lets Highway target AVX2/AVX-512 rather than the x86-64
baseline — a large additional speedup, but the resulting binary only runs on
machines like the one that built it.

Check what a given install actually has:

```python
>>> mtl5.build_info()
{'blas': False, 'lapack': False, 'native_fast_gemm': True,
 'highway_simd': True, 'kpu': False}
>>> mtl5.get_backend()
'native'
```

`set_backend()` validates against this build rather than silently accepting a
backend that was never compiled in; backend selection itself is compile-time.

### Benchmarking the number systems

The emulated formats are much slower than hardware float, and *how much* is the
number that decides whether a mixed-precision experiment fits in a night.
`benchmarks/bench_blas.py` measures it — `dot`, `gemv`, `gemm`, `lu` and `qr`
across the bound number systems — and reports **slowdown against float64**
rather than raw GFLOP/s, because that is the figure experiment planning is
built from.

```bash
python benchmarks/bench_blas.py            # default sweep
python benchmarks/bench_blas.py --quick    # smallest sizes, seconds
python benchmarks/bench_blas.py --json out.json
```

For scale, one recorded run (x86_64, single thread, no `-march=native`), GEMM at
n=64: `takum32` 434x float64, `cfloat32` 570x, `posit32` 4024x. Those are an
example from one host, not portable constants — every report embeds
`build_info()`, thread count and architecture for that reason.

`cfloat32` is the one to compare against rather than `float`: it is IEEE
binary32 through the same emulation layer, so `posit32 / cfloat32` separates
what a *format* costs from what *emulation* costs.

**[`benchmarks/README.md`](benchmarks/README.md)** documents the methodology and
how to read the output — including two traps: the small-`n` rows understate
every ratio, and a narrow format can time perfectly normally while computing
nothing useful.

## Development

```bash
pip install -e ".[dev]"
pytest -v
```
