mtl5-python

Python bindings for MTL5 — NumPy/SciPy/JAX/PyTorch interop with hardware accelerator dispatch.

Built with nanobind for minimal overhead and zero-copy array interop.

Install

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:

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 — needs Python 3.10+ and a C++20 compiler (GCC 12+, Clang 15+, MSVC 2022):

pip install .

See Development for an editable install.

Quick start

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 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=

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:

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

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:

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:

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:

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

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:

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:

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:

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:

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:

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:

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.

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_asserts 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:

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:

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:

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:

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:

# 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), 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:

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:

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:

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:

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:

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

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:

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:

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:

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:

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‖₂

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:

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:

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:

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:

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 by default. Everything else is opt-in, since it needs a library on the build machine or makes the binary non-portable:

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:

>>> 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.

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 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

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