Solves dense systems and works out inverses, determinants, factorizations, eigenvalues, norms and condition numbers in a script.
How the precision is kept
All matrices are Float64Array in column-major order and every factorization runs in 64-bit. LU uses partial pivoting, QR uses Householder reflections with the sign chosen to avoid cancellation, and the symmetric eigenvalue routine is tred2 and tql2. The matrix 2-norm and the condition number are taken from the eigenvalues of AᵀA, which limits them for badly conditioned matrices.
Open item. Scale the least-squares rank test to the size of R, and compute cond and the matrix 2-norm from a singular value decomposition instead of the eigenvalues of AᵀA.
Calculations under this policy
13 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Each entry is a plain 64-bit sum of products, added in order of p; terms with a zero in B are skipped. No compensated summation.
cs-lib-linalg.js
lines 436 to 445
PartCraft 0.8.1
export function matmul(A, B) {
const [n, k] = size2(A), [k2, m] = size2(B);
if (k !== k2) throw err(t('script.linalg.inner_sizes', { a: `${n}×${k}`, b: `${k2}×${m}` }));
const a = A.data, b = B.data, c = new Float64Array(n * m);
for (let j = 0; j < m; j++) {
const cj = j * n;
for (let p = 0; p < k; p++) { const bv = b[p + j * k]; if (bv === 0) continue; const ap = p * n; for (let i = 0; i < n; i++) c[cj + i] += a[ap + i] * bv; }
}
return mat(n, m, c);
}
LU factorization
64-bit
Factors a square matrix with row exchanges.
P A = L U
Partial pivoting picks the largest absolute entry in each column. A column whose best pivot is at most 1e-15 times the largest absolute entry of the whole matrix (or exactly zero) is skipped and the result is flagged singular. The threshold is relative to the matrix, not to the column.
cs-lib-linalg.js
lines 447 to 464
PartCraft 0.8.1
/** LU with partial pivoting (in place on a copy) → { LU, piv (row order), sign, singular } */
export function luFactor(A) {
const [n, m] = size2(A);
if (n !== m) throw err(t('script.linalg.square', { size: `${n}×${m}` }));
const a = Float64Array.from(A.data), piv = Int32Array.from({ length: n }, (_, i) => i);
let sign = 1, singular = false;
let scale = 0; for (const v of a) scale = Math.max(scale, Math.abs(v));
for (let c = 0; c < n; c++) {
let p = c, best = Math.abs(a[c + c * n]);
for (let r = c + 1; r < n; r++) { const v = Math.abs(a[r + c * n]); if (v > best) { best = v; p = r; } }
if (best <= 1e-15 * scale || best === 0) { singular = true; continue; }
if (p !== c) { for (let k = 0; k < n; k++) { const t = a[c + k * n]; a[c + k * n] = a[p + k * n]; a[p + k * n] = t; } const t = piv[c]; piv[c] = piv[p]; piv[p] = t; sign = -sign; }
const d = a[c + c * n];
for (let r = c + 1; r < n; r++) a[r + c * n] /= d;
for (let k = c + 1; k < n; k++) { const f = a[c + k * n]; if (f === 0) continue; const ck = k * n, cc = c * n; for (let r = c + 1; r < n; r++) a[r + ck] -= a[r + cc] * f; }
}
return { n, LU: a, piv, sign, singular };
}
Solve with an LU factorization
64-bit
Forward and back substitution for one or more right-hand sides.
Substitution in 64-bit, one column of the right-hand side at a time. There is no iterative refinement, so the error grows with the condition number of the matrix.
cs-lib-linalg.js
lines 465 to 474
PartCraft 0.8.1
function luSolveInto(F, B) { // B dense n×m → X
const { n, LU, piv } = F, [, m] = size2(B), X = new Float64Array(n * m);
for (let j = 0; j < m; j++) {
const x = X.subarray(j * n, (j + 1) * n), b = B.data.subarray(j * n, (j + 1) * n);
for (let i = 0; i < n; i++) x[i] = b[piv[i]];
for (let c = 0; c < n; c++) { const v = x[c]; if (v) for (let r = c + 1; r < n; r++) x[r] -= LU[r + c * n] * v; }
for (let c = n - 1; c >= 0; c--) { x[c] /= LU[c + c * n]; const v = x[c]; if (v) for (let r = 0; r < c; r++) x[r] -= LU[r + c * n] * v; }
}
return mat(n, m, X);
}
Cholesky factorization
64-bit
Factors a symmetric positive definite matrix as RᵀR.
R_jj = √(A_jj − Σ_k R_kj²)
Reads the upper triangle only. If a diagonal term is not greater than zero after the subtraction the function returns null (not positive definite) instead of taking the square root of a negative number.
cs-lib-linalg.js
lines 475 to 488
PartCraft 0.8.1
/** Cholesky A = Rᵀ R (upper R) or null when not positive definite. */
export function cholFactor(A) {
const [n, m] = size2(A);
if (n !== m) throw err(t('script.linalg.square', { size: `${n}×${m}` }));
const a = A.data, R = new Float64Array(n * n);
for (let j = 0; j < n; j++) {
let d = a[j + j * n];
for (let k = 0; k < j; k++) d -= R[k + j * n] * R[k + j * n];
if (!(d > 0)) return null;
const rjj = Math.sqrt(d); R[j + j * n] = rjj;
for (let i = j + 1; i < n; i++) { let s = a[j + i * n]; for (let k = 0; k < j; k++) s -= R[k + j * n] * R[k + i * n]; R[j + i * n] = s / rjj; }
}
return mat(n, n, R);
}
Solve with a Cholesky factorization
64-bit
Solves RᵀR X = B by two triangular substitutions.
Two 64-bit substitutions per right-hand side column; no refinement.
cs-lib-linalg.js
lines 490 to 498
PartCraft 0.8.1
function cholSolve(R, B) { // Rᵀ R X = B
const n = R.shape[0], [, m] = size2(B), X = Float64Array.from(B.data), r = R.data;
for (let j = 0; j < m; j++) {
const x = X.subarray(j * n, (j + 1) * n);
for (let i = 0; i < n; i++) { let s = x[i]; for (let k = 0; k < i; k++) s -= r[k + i * n] * x[k]; x[i] = s / r[i + i * n]; }
for (let i = n - 1; i >= 0; i--) { let s = x[i]; for (let k = i + 1; k < n; k++) s -= r[i + k * n] * x[k]; x[i] = s / r[i + i * n]; }
}
return mat(n, m, X);
}
QR factorization
64-bit
Factors a matrix as an orthogonal matrix times an upper triangular one.
H = I − β v vᵀ, β = 2 ÷ (vᵀv)
Householder reflections. The reflector's sign is chosen opposite to the diagonal entry so the first component is formed by an addition of like signs, not a cancelling subtraction. Column norms are plain sums of squares without scaling, so entries near the square root of the largest or smallest 64-bit number can overflow or underflow.
cs-lib-linalg.js
lines 499 to 524
PartCraft 0.8.1
/** Householder QR → { Q (n×n or n×k econ), R } */
export function qrFactor(A, econ = false) {
const [n, m] = size2(A), a = Float64Array.from(A.data), k = Math.min(n, m), vs = [], betas = [];
for (let c = 0; c < k; c++) {
let nrm = 0; for (let r = c; r < n; r++) nrm += a[r + c * n] ** 2;
nrm = Math.sqrt(nrm);
const v = new Float64Array(n - c);
if (nrm === 0) { vs.push(v); betas.push(0); continue; }
const alpha = a[c + c * n] > 0 ? -nrm : nrm;
for (let r = c; r < n; r++) v[r - c] = a[r + c * n];
v[0] -= alpha;
let vv = 0; for (const x of v) vv += x * x;
const beta = vv ? 2 / vv : 0;
for (let j = c; j < m; j++) { let s = 0; for (let r = c; r < n; r++) s += v[r - c] * a[r + j * n]; s *= beta; for (let r = c; r < n; r++) a[r + j * n] -= s * v[r - c]; }
vs.push(v); betas.push(beta);
}
const qc = econ ? k : n, Q = new Float64Array(n * qc);
for (let i = 0; i < qc; i++) Q[i + i * n] = 1;
for (let c = k - 1; c >= 0; c--) {
const v = vs[c], beta = betas[c]; if (!beta) continue;
for (let j = 0; j < qc; j++) { let s = 0; for (let r = c; r < n; r++) s += v[r - c] * Q[r + j * n]; s *= beta; for (let r = c; r < n; r++) Q[r + j * n] -= s * v[r - c]; }
}
const rr = econ ? k : n, R = new Float64Array(rr * m);
for (let j = 0; j < m; j++) for (let i = 0; i <= Math.min(j, rr - 1); i++) R[i + j * rr] = a[i + j * n];
return { Q: mat(n, qc, Q), R: mat(rr, m, R) };
}
Symmetric eigenvalues and vectors
64-bit
All eigenvalues, and optionally eigenvectors, of a symmetric matrix, in ascending order.
Householder tridiagonalization followed by implicit QL with shifts. An off-diagonal entry is treated as zero when it is at most 2^-52 times the largest |d| + |e| seen so far. More than 60 iterations for one eigenvalue raises an error instead of returning a value.
cs-lib-linalg.js
lines 525 to 614
PartCraft 0.8.1
/** Symmetric eigen-decomposition (tred2 + tql2) → { d ascending Float64Array, V column-major n×n } */
export function symEigen(A, wantV = true) {
const [n, m] = size2(A);
if (n !== m) throw err(t('script.linalg.eig_square', { size: `${n}×${m}` }));
const V = Float64Array.from(A.data), d = new Float64Array(n), e = new Float64Array(n);
const v = (i, j) => V[i + j * n];
// tred2 (JAMA's port of EISPACK), V symmetric in, orthogonal out
for (let j = 0; j < n; j++) d[j] = v(n - 1, j);
for (let i = n - 1; i > 0; i--) {
let scale = 0, h = 0;
for (let k = 0; k < i; k++) scale += Math.abs(d[k]);
if (scale === 0) {
e[i] = d[i - 1];
for (let j = 0; j < i; j++) { d[j] = v(i - 1, j); V[i + j * n] = 0; V[j + i * n] = 0; }
} else {
for (let k = 0; k < i; k++) { d[k] /= scale; h += d[k] * d[k]; }
let f = d[i - 1], g = Math.sqrt(h);
if (f > 0) g = -g;
e[i] = scale * g; h -= f * g; d[i - 1] = f - g;
for (let j = 0; j < i; j++) e[j] = 0;
for (let j = 0; j < i; j++) {
f = d[j]; V[j + i * n] = f; g = e[j] + v(j, j) * f;
for (let k = j + 1; k <= i - 1; k++) { g += v(k, j) * d[k]; e[k] += v(k, j) * f; }
e[j] = g;
}
f = 0;
for (let j = 0; j < i; j++) { e[j] /= h; f += e[j] * d[j]; }
const hh = f / (h + h);
for (let j = 0; j < i; j++) e[j] -= hh * d[j];
for (let j = 0; j < i; j++) {
f = d[j]; g = e[j];
for (let k = j; k <= i - 1; k++) V[k + j * n] -= f * e[k] + g * d[k];
d[j] = v(i - 1, j); V[i + j * n] = 0;
}
}
d[i] = h;
}
for (let i = 0; i < n - 1; i++) {
V[n - 1 + i * n] = v(i, i); V[i + i * n] = 1;
const h = d[i + 1];
if (h !== 0) {
for (let k = 0; k <= i; k++) d[k] = v(k, i + 1) / h;
for (let j = 0; j <= i; j++) {
let g = 0; for (let k = 0; k <= i; k++) g += v(k, i + 1) * v(k, j);
for (let k = 0; k <= i; k++) V[k + j * n] -= g * d[k];
}
}
for (let k = 0; k <= i; k++) V[k + (i + 1) * n] = 0;
}
for (let j = 0; j < n; j++) { d[j] = v(n - 1, j); V[n - 1 + j * n] = 0; }
if (n) V[n - 1 + (n - 1) * n] = 1;
e[0] = 0;
// tql2
for (let i = 1; i < n; i++) e[i - 1] = e[i];
if (n) e[n - 1] = 0;
let f = 0, tst1 = 0;
const eps = 2 ** -52;
for (let l = 0; l < n; l++) {
tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l]));
let m2 = l;
while (m2 < n) { if (Math.abs(e[m2]) <= eps * tst1) break; m2++; }
if (m2 > l) {
let iter = 0;
do {
if (++iter > 60) throw err(t('script.linalg.eig_noconv'));
let g = d[l], p = (d[l + 1] - g) / (2 * e[l]), r = Math.hypot(p, 1);
if (p < 0) r = -r;
d[l] = e[l] / (p + r); d[l + 1] = e[l] * (p + r);
const dl1 = d[l + 1]; let h = g - d[l];
for (let i = l + 2; i < n; i++) d[i] -= h;
f += h;
p = d[m2];
let c = 1, c2 = c, c3 = c, el1 = e[l + 1], s = 0, s2 = 0;
for (let i = m2 - 1; i >= l; i--) {
c3 = c2; c2 = c; s2 = s;
g = c * e[i]; h = c * p; r = Math.hypot(p, e[i]);
e[i + 1] = s * r; s = e[i] / r; c = p / r; p = c * d[i] - s * g; d[i + 1] = h + s * (c * g + s * d[i]);
if (wantV) for (let k = 0; k < n; k++) { h = V[k + (i + 1) * n]; V[k + (i + 1) * n] = s * V[k + i * n] + c * h; V[k + i * n] = c * V[k + i * n] - s * h; }
}
p = -s * s2 * c3 * el1 * e[l] / dl1; e[l] = s * p; d[l] = c * p;
} while (Math.abs(e[l]) > eps * tst1);
}
d[l] += f; e[l] = 0;
}
// sort ascending
const ord = Array.from({ length: n }, (_, i) => i).sort((a, b) => d[a] - d[b]);
const ds = Float64Array.from(ord, (i) => d[i]), Vs = new Float64Array(n * n);
ord.forEach((i, j) => Vs.set(V.subarray(i * n, (i + 1) * n), j * n));
return { d: ds, V: Vs };
}
Generalized symmetric eigenproblem
64-bit
Eigenvalues and vectors of K v = λ M v for symmetric K and positive definite M.
M = RᵀR; C = R⁻ᵀ K R⁻¹; C y = λ y; v = R⁻¹ y
The reduction uses triangular solves with the Cholesky factor of M, then averages C with its transpose to remove rounding asymmetry. The accuracy of the small eigenvalues depends on the condition number of M, because R⁻¹ is applied on both sides.
cs-lib-linalg.js
lines 615 to 631
PartCraft 0.8.1
/** Generalized symmetric-definite K v = λ M v (dense) → { d, V } with Vᵀ M V = I. */
export function symGenEigen(K, M) {
const R = cholFactor(M);
if (!R) throw err(t('script.linalg.eig_second_spd'));
const n = K.shape[0], r = R.data;
// C = R⁻ᵀ K R⁻¹
const Y = Float64Array.from(K.data); // Y = R⁻ᵀ K (solve Rᵀ Y = K column by column)
for (let j = 0; j < n; j++) for (let i = 0; i < n; i++) { let s = Y[i + j * n]; for (let k = 0; k < i; k++) s -= r[k + i * n] * Y[k + j * n]; Y[i + j * n] = s / r[i + i * n]; }
const C = new Float64Array(n * n); // C = Y R⁻¹ → Cᵀ = R⁻ᵀ Yᵀ
for (let j = 0; j < n; j++) for (let i = 0; i < n; i++) C[j + i * n] = Y[i + j * n];
for (let j = 0; j < n; j++) for (let i = 0; i < n; i++) { let s = C[i + j * n]; for (let k = 0; k < i; k++) s -= r[k + i * n] * C[k + j * n]; C[i + j * n] = s / r[i + i * n]; }
for (let j = 0; j < n; j++) for (let i = j + 1; i < n; i++) { const m = (C[i + j * n] + C[j + i * n]) / 2; C[i + j * n] = C[j + i * n] = m; }
const { d, V } = symEigen(mat(n, n, C));
// v = R⁻¹ y
for (let j = 0; j < n; j++) for (let i = n - 1; i >= 0; i--) { let s = V[i + j * n]; for (let k = i + 1; k < n; k++) s -= r[i + k * n] * V[k + j * n]; V[i + j * n] = s / r[i + i * n]; }
return { d, V };
}
Choice of solver for A x = b
64-bit
Picks the method for the backslash operator and linsolve, for dense and sparse matrices.
Dense: Cholesky when the matrix is symmetric within 1e-12 of its largest entry and the factorization succeeds, otherwise LU; a non-square matrix is solved by QR least squares, where a diagonal entry of R below 1e-14 in absolute value raises a rank error whatever the scale of the matrix. Sparse: dense up to 200 unknowns, envelope Cholesky when symmetric and it fits, else PCG with tolerance 1e-10, else dense LU up to 3000 unknowns or BiCGSTAB with tolerance 1e-10; an iterative solve that does not converge raises an error.
cs-lib-linalg.js
lines 634 to 682
PartCraft 0.8.1
/** A \ B for dense or sparse A → { X (tensor), method, info } */
export function solve(A, B, { check = () => {}, tol = 1e-10, prec = 'auto' } = {}) {
B = asTensor(B, t('script.linalg.what.rhs'));
const [bn, bm] = size2(B);
if (isSparse(A)) {
if (A.n !== A.m) throw err(t('script.linalg.solve_square', { size: `${A.n}×${A.m}` }));
if (bn !== A.n) throw err(t('script.linalg.solve_rhs', { b: bn, n: A.n }));
const sym = spIsSymmetric(A);
if (A.n <= 200) return solve(toDense(A), B, { check });
if (sym) {
let F = null;
try { F = envelopeCholesky(A, { check }); } catch (e) { if (!e.notPD) throw e; F = false; }
if (F) { const X = new Float64Array(bn * bm); for (let j = 0; j < bm; j++) F.solve(B.data.subarray(j * bn, (j + 1) * bn), X.subarray(j * bn, (j + 1) * bn)); return { X: mat(bn, bm, X), method: F.kind, info: { envelope: F.envelope, ms: F.ms } }; }
if (F === null) { // too big for the envelope: PCG
const M = makePrec(A, prec), X = new Float64Array(bn * bm), its = [];
for (let j = 0; j < bm; j++) {
const r = pcgSolve(A, B.data.subarray(j * bn, (j + 1) * bn), { tol, prec: M, check });
if (!r.converged) throw err(t('script.linalg.solve_pcg', { res: r.residual.toExponential(2), its: r.iterations }));
X.set(r.x, j * bn); its.push(r.iterations);
}
return { X: mat(bn, bm, X), method: `pcg + ${M.kind}`, // i18n-ignore: a method name scripts read from info.method
info: { iterations: its } };
}
// symmetric but not positive definite: fall through to the general path
}
if (A.n <= DENSE_MAX) return solve(toDense(A), B, { check });
const X = new Float64Array(bn * bm);
for (let j = 0; j < bm; j++) {
const r = bicgstabSolve(A, B.data.subarray(j * bn, (j + 1) * bn), { tol, check });
if (!r.converged) throw err(t('script.linalg.solve_bicg', { res: r.residual.toExponential(2) }));
X.set(r.x, j * bn);
}
return { X: mat(bn, bm, X), method: 'bicgstab + ilu0', info: {} }; // i18n-ignore: a method name scripts read
}
A = asTensor(A, t('script.linalg.what.matrix'));
const [n, m] = size2(A);
if (n === 1 && m === 1) return { X: mat(bn, bm, B.data.map((v) => v / A.data[0])), method: 'scalar', info: {} };
if (bn !== n) throw err(t('script.linalg.solve_rhs', { b: bn, n }));
if (n !== m) { // least squares via QR (full column rank assumed)
if (n < m) throw err(t('script.linalg.solve_under'));
const { Q, R } = qrFactor(A, true), QtB = matmul(transposeD(Q), B), X = new Float64Array(m * bm);
for (let j = 0; j < bm; j++) for (let i = m - 1; i >= 0; i--) { let s = QtB.data[i + j * m]; for (let k = i + 1; k < m; k++) s -= R.data[i + k * m] * X[k + j * m]; const d = R.data[i + i * m]; if (Math.abs(d) < 1e-14) throw err(t('script.linalg.solve_rank')); X[i + j * m] = s / d; }
return { X: mat(m, bm, X), method: 'least squares (QR)', info: {} }; // i18n-ignore: a method name scripts read
}
if (isSymD(A)) { const R = cholFactor(A); if (R) return { X: cholSolve(R, B), method: 'Cholesky', info: {} }; }
const F = luFactor(A);
if (F.singular) throw err(t('script.linalg.solve_singular'));
return { X: luSolveInto(F, B), method: 'LU', info: {} };
}
Matrix inverse
64-bit
The inverse of a square matrix.
LU factorization, then one substitution per column of the identity. A matrix flagged singular by the LU threshold raises an error.
cs-lib-linalg.js
lines 792 to 797
PartCraft 0.8.1
export function bi_inv([A]) {
const a = T1(A, t('script.linalg.what.matrix')), [n] = size2(a), F = luFactor(a);
if (F.singular) throw err(t('script.linalg.inv_singular'));
const I = new Float64Array(n * n); for (let i = 0; i < n; i++) I[i + i * n] = 1;
return luSolveInto(F, mat(n, n, I));
}
Determinant
64-bit
The determinant of a square matrix.
det A = sign × Π U_ii
A plain running product of the LU pivots with the sign of the row exchanges. A matrix flagged singular by the LU threshold (a pivot at most 1e-15 of the largest entry) returns exactly 0, not the small computed value. The product is not scaled, so it can overflow or underflow for large matrices.
cs-lib-linalg.js
lines 798 to 803
PartCraft 0.8.1
export function bi_det([A]) {
const F = luFactor(T1(A, t('script.linalg.what.matrix')));
if (F.singular) return 0;
let d = F.sign; for (let i = 0; i < F.n; i++) d *= F.LU[i + i * F.n];
return d;
}
Vector and matrix norms
64-bit
The 1, 2, infinity, p and Frobenius norms.
‖x‖₂ = √(Σ xᵢ²); ‖A‖₂ = √(λ_max(AᵀA))
The vector 2-norm and the Frobenius norm are a plain sum of squares, added in order, without scaling against overflow. The matrix 2-norm is the square root of the largest eigenvalue of AᵀA, clamped at zero before the square root.
cs-lib-linalg.js
lines 835 to 849
PartCraft 0.8.1
export function normOf(A, p = 2) {
if (isSparse(A)) A = toDense(A);
const a = T1(A), [n, m] = size2(a), d = a.data;
if (p === 'fro' || ((n === 1 || m === 1) && p === 2)) { let s = 0; for (const v of d) s += v * v; return Math.sqrt(s); }
if (n === 1 || m === 1) {
if (p === Infinity || p === 'inf') { let s = 0; for (const v of d) s = Math.max(s, Math.abs(v)); return s; }
if (p === -Infinity) { let s = Infinity; for (const v of d) s = Math.min(s, Math.abs(v)); return s; }
if (typeof p !== 'number' || !(p > 0)) throw err(t('script.linalg.norm_p'));
let s = 0; for (const v of d) s += Math.abs(v) ** p; return s ** (1 / p);
}
if (p === 1) { let s = 0; for (let j = 0; j < m; j++) { let c = 0; for (let i = 0; i < n; i++) c += Math.abs(d[i + j * n]); s = Math.max(s, c); } return s; }
if (p === Infinity || p === 'inf') { let s = 0; for (let i = 0; i < n; i++) { let c = 0; for (let j = 0; j < m; j++) c += Math.abs(d[i + j * n]); s = Math.max(s, c); } return s; }
if (p === 2) { const { d: ev } = symEigen(matmul(transposeD(a), a), false); return Math.sqrt(Math.max(0, ev[ev.length - 1])); }
throw err(t('script.linalg.norm_matrix_p'));
}
Condition number
64-bit
The 2-norm condition number of a matrix.
cond A = √(λ_max(AᵀA) ÷ λ_min(AᵀA))
The eigenvalues of AᵀA are the squared singular values, so the smallest one is computed against rounding error of the order of 2^-52 times the largest. Condition numbers near or above 1e8 are therefore not resolved; a smallest eigenvalue that comes out zero or negative returns Infinity.
cs-lib-linalg.js
lines 851 to 856
PartCraft 0.8.1
export function bi_cond([A]) {
const a = isSparse(A) ? toDense(A) : T1(A, t('script.linalg.what.matrix'));
const { d } = symEigen(matmul(transposeD(a), a), false);
const lo = Math.max(0, d[0]), hi = Math.max(0, d[d.length - 1]);
return lo === 0 ? Infinity : Math.sqrt(hi / lo);
}
In detail
Method
LU with partial pivoting; Cholesky A = RᵀR; Householder QR; Householder tridiagonalization and implicit QL (the code's comments name EISPACK tred2/tql2 and JAMA's port); generalized symmetric-definite problem by Cholesky reduction.
Number types
64-bit numbers and Float64Array; Int32Array for pivot order.
Tolerances
LU flags a singular matrix when a pivot is at most 1e-15 times the largest absolute entry of the matrix. Symmetry test: 1e-12 (in the solver) or 1e-10 (chol, eig) times the largest absolute entry. Least-squares rank test: an R diagonal entry below 1e-14, absolute. tql2 deflates at 2^-52 times the running scale and stops with an error after 60 iterations for one eigenvalue.
Limits
Eigenvalues are for symmetric matrices only. The least-squares rank test is absolute, not scaled to the matrix. norm(A, 2) and cond(A) form AᵀA.
Approach: 64-bit, direct
Every value is a 64-bit number from input to result. Formulas are closed form or a direct solve, and every comparison uses a stated tolerance.