A linear system is a question about data, and some data answers no: a
singular matrix has no inverse, an indefinite one no Cholesky factor. Those are
not program errors — no assertion was violated, no index was wrong — so the
linalg unit is born on the failure channel:
imp
Solve, Invert, Cholesky, Qr : linalg;
LinalgError, QrDecomposition : linalg;
Singular, NotPositiveDefinite : linalg;Every routine fails LinalgError, and the members say what the data was:
Singular from an exactly zero pivot — the reference implementations’ own
rule, no epsilon to tune — NotPositiveDefinite from Cholesky’s first
non-positive pivot, and ShapeMismatch where a routine meets sizes that fit
no question at all.
The example is examples/Linalg
in the tutorial repository.
Solve
x := Solve(a, b) on fail e do
begin
WriteLn("unexpected failure");
leave;
end;Solve answers the x with a·x = b through one partial-pivot LU
factorization and one substitution pass. The pivoting is part of the answer’s
quality, not an option — the example’s system deliberately starts with a zero
on the diagonal, so the row exchange must happen for the answer to exist.
Invert, and when not to
Invert answers the whole inverse: the factorization once, then one
substitution per identity column. The example verifies it the strongest way
available — with the shipped operators themselves:
inv := Invert(a) on fail e do
begin
WriteLn("unexpected failure");
leave;
end;
idm := a * inv; { multiplies back to the identity }Where a program only ever multiplies the inverse onto vectors, Solve per
vector is the cheaper and better-conditioned road; Invert exists for the
answers that are themselves the matrix.
Cholesky as a question
Cholesky answers the lower factor with a = l·lᵀ, reading only the lower
triangle — the symmetric contract’s own convention. Because the first
non-positive pivot refuses as NotPositiveDefinite, the call doubles as the
cheapest definiteness test there is:
l := Cholesky(a) on fail e do
begin
if e = NotPositiveDefinite then
WriteLn("not positive definite — and that is the answer");
end;Qr
Qr answers the reduced orthogonal-triangular pair by Householder
reflections — the numerically sound road — as one record:
parts := Qr(a) on fail e do
begin
WriteLn("unexpected failure");
leave;
end;
back := parts.q * parts.r; { restores the argument }(Inside a function that itself declares fails LinalgError, the shorter
on fail leave forwards the failure instead.)
For a 3×2 argument, parts.q is 3×2 with orthonormal columns and parts.r is
2×2 upper-triangular, with exact zeros below the diagonal — the triangle is a
stated fact of the answer, not a rounding accident.
The value model is the working storage
Every classic factorization overwrites its matrix in place. In Mica that costs no caution at all: a parameter arrives as the callee’s own deep copy, so the elimination works in place on that copy while the caller’s matrix never moves. What LAPACK documents as “on exit, A has been overwritten” is here simply not observable.
Every fold runs ascending in a fixed order, so results are bit-identical across optimization tiers and architectures — the same promise the operators keep. No BLAS is linked; the compiler’s own code generation carries the unit.