//! NumPy-parity tests for the slice-level linear-algebra kernels.
//!
//! Every `data::*` constant below was produced once by `numpy 2.3.2` (float64 unless noted) from the
//! small inputs listed next to it (inputs are rounded to 3 decimals or come from the integer formulas
//! visible in the tests); nothing here is recomputed by the code under test. The one-off generator
//! scripts are not part of the repository.
//! Where an exact answer exists (Hilbert determinant and inverse) it was computed with
//! Python `fractions.Fraction`. Tolerances are normwise: `max|actual - expected| <= tol * max(1, max|expected|)`.
use super::kernels::*;
use data::*;
// ---------------------------------------------------------------------------
// helpers
// ---------------------------------------------------------------------------
fn max_abs(values: &[f64]) -> f64 {
values.iter().fold(0.0, |best, v| f64::max(best, v.abs()))
}
#[track_caller]
fn assert_close(actual: &[f64], expected: &[f64], tol: f64) {
assert_eq!(actual.len(), expected.len(), "length mismatch");
let scale = max_abs(expected).max(1.0);
let worst = actual
.iter()
.zip(expected)
.fold(0.0_f64, |w, (a, e)| w.max((a - e).abs()));
assert!(
worst <= tol * scale,
"max abs error {worst:e} > {:e}\nactual = {actual:?}\nexpected = {expected:?}",
tol * scale
);
}
#[track_caller]
fn assert_rel(actual: f64, expected: f64, tol: f64) {
assert!(
(actual - expected).abs() <= tol * expected.abs(),
"{actual:e} vs {expected:e} (rel tol {tol:e})"
);
}
#[track_caller]
fn assert_close_f32(actual: &[f32], expected: &[f64], tol: f32) {
assert_eq!(actual.len(), expected.len());
let scale = (max_abs(expected) as f32).max(1.0);
for (a, e) in actual.iter().zip(expected) {
assert!(
(a - *e as f32).abs() <= tol * scale,
"{a} vs {e}\nactual = {actual:?}"
);
}
}
fn mm(a: &[f64], b: &[f64], m: usize, k: usize, n: usize) -> Vec<f64> {
gemm(m, n, k, a, k as isize, 1, b, n as isize, 1)
}
fn transpose<T: Copy>(a: &[T], rows: usize, cols: usize) -> Vec<T> {
let mut out = Vec::with_capacity(a.len());
for j in 0..cols {
for i in 0..rows {
out.push(a[i * cols + j]);
}
}
out
}
fn eye(n: usize) -> Vec<f64> {
let mut out = vec![0.0; n * n];
for i in 0..n {
out[i * n + i] = 1.0;
}
out
}
fn hilbert(n: usize) -> Vec<f64> {
(0..n * n)
.map(|index| 1.0 / ((index / n + index % n + 1) as f64))
.collect()
}
fn sub(a: &[f64], b: &[f64]) -> Vec<f64> {
a.iter().zip(b).map(|(x, y)| x - y).collect()
}
fn scaled(a: &[f64], factor: f64) -> Vec<f64> {
a.iter().map(|x| x * factor).collect()
}
fn to_f32(a: &[f64]) -> Vec<f32> {
a.iter().map(|&x| x as f32).collect()
}
/// Columns of the `rows x cols` matrix `q` are orthonormal: `qT q = I`.
#[track_caller]
fn assert_orthonormal_columns(q: &[f64], rows: usize, cols: usize, tol: f64) {
let gram = mm(&transpose(q, rows, cols), q, cols, rows, cols);
assert_close(&gram, &eye(cols), tol);
}
/// Rows of the `rows x cols` matrix `q` are orthonormal: `q qT = I`.
#[track_caller]
fn assert_orthonormal_rows(q: &[f64], rows: usize, cols: usize, tol: f64) {
let gram = mm(q, &transpose(q, rows, cols), rows, cols, rows);
assert_close(&gram, &eye(rows), tol);
}
/// The component of largest magnitude in column `col` of a row-major `rows x cols` matrix is positive.
#[track_caller]
fn assert_sign_rule(m: &[f64], rows: usize, cols: usize, col: usize) {
let mut best = 0.0_f64;
let mut value = 0.0_f64;
for row in 0..rows {
let v = m[row * cols + col];
if v.abs() > best {
best = v.abs();
value = v;
}
}
assert!(value >= 0.0, "column {col} violates the sign rule: {value}");
}
/// Full SVD contract: values match NumPy, `U S Vt = A`, orthonormal factors, descending, sign rule.
#[track_caller]
fn check_svd(a: &[f64], rows: usize, cols: usize, s_numpy: &[f64], tol: f64) {
let (u, s, vt, k) = svd(a, rows, cols).unwrap();
assert_eq!(k, rows.min(cols));
assert_eq!(s.len(), k);
assert_close(&s, s_numpy, tol);
assert!(s.windows(2).all(|w| w[0] >= w[1]), "not descending: {s:?}");
let mut us = u.clone();
for r in 0..rows {
for c in 0..k {
us[r * k + c] *= s[c];
}
}
assert_close(&mm(&us, &vt, rows, k, cols), a, tol);
assert_orthonormal_columns(&u, rows, k, tol);
assert_orthonormal_rows(&vt, k, cols, tol);
for c in 0..k {
assert_sign_rule(&u, rows, k, c);
}
}
/// Full eigh contract: values match NumPy, `A v = lambda v`, orthonormal V, ascending, sign rule.
#[track_caller]
fn check_eigh(a: &[f64], n: usize, w_numpy: &[f64], tol: f64) -> Vec<f64> {
let (w, v) = eigh(a, n).unwrap();
assert_close(&w, w_numpy, tol);
assert!(w.windows(2).all(|p| p[0] <= p[1]), "not ascending: {w:?}");
let av = mm(a, &v, n, n, n);
let mut lv = v.clone();
for r in 0..n {
for c in 0..n {
lv[r * n + c] *= w[c];
}
}
assert_close(&av, &lv, tol);
assert_orthonormal_columns(&v, n, n, tol);
for c in 0..n {
assert_sign_rule(&v, n, n, c);
}
v
}
/// Two eigenvectors/singular vectors agree up to sign.
#[track_caller]
fn assert_columns_match_up_to_sign(
actual: &[f64],
expected: &[f64],
rows: usize,
cols: usize,
tol: f64,
) {
for c in 0..cols {
let dot: f64 = (0..rows)
.map(|r| actual[r * cols + c] * expected[r * cols + c])
.sum();
let sign = if dot < 0.0 { -1.0 } else { 1.0 };
for r in 0..rows {
assert!(
(sign * actual[r * cols + c] - expected[r * cols + c]).abs() <= tol,
"column {c} differs: {actual:?} vs {expected:?}"
);
}
}
}
/// Constants generated once with NumPy 2.3.2 (see the module docs).
#[allow(dead_code, clippy::approx_constant)]
#[rustfmt::skip]
mod data {
// ---- GEMM ----
pub const G33_A: [f64; 9] = [0.346, 0.822, 0.33, -1.303, 0.905, 0.446, -0.537, 0.581, 0.365];
pub const G33_B: [f64; 9] = [0.189, -0.523, -0.413, -2.441, 1.8, 1.144, -0.325, 0.774, 0.281];
pub const G33_C: [f64; 9] = [-2.048358, 1.554062, 0.8901999999999999, -2.6003220000000002, 2.6556729999999997, 1.698785, -1.6383389999999998, 1.609161, 0.98901];
pub const G34_A: [f64; 12] = [2.041, -2.556, 0.418, -0.568, -0.453, -0.216, -2.02, -0.232, -0.865, 3.323, 0.226, -0.353];
pub const G42_B: [f64; 8] = [-0.652, -0.175, 1.664, 0.659, -1.641, -0.005, -0.623, 0.149];
pub const G32_C: [f64; 6] = [-5.915990000000001, -2.128301, 3.395288, -0.087537, 5.942505, 2.2875050000000003];
pub const G15_A: [f64; 5] = [-0.802, -1.324, -0.248, 0.42, 1.136];
pub const G53_B: [f64; 15] = [1.053, 1.776, -2.553, -0.138, 1.014, 1.352, 0.654, 1.497, 0.29, 0.551, 0.179, -1.074, -0.847, 0.38, -0.58];
pub const G13_C: [f64; 3] = [-1.5547579999999999, -2.631284, -0.924422];
pub const G51_A: [f64; 5] = [0.001, 0.299, -0.274, -0.891, -0.455];
pub const G14_B: [f64; 4] = [-1.738, -1.337, -1.361, -0.352];
pub const G54_C: [f64; 20] = [-0.001738, -0.001337, -0.001361, -0.000352, -0.519662, -0.399763, -0.406939, -0.105248, 0.476212, 0.366338, 0.372914, 0.096448, 1.548558, 1.191267, 1.212651, 0.31363199999999997, 0.79079, 0.608335, 0.619255, 0.16016];
pub const G_AT: [f64; 12] = [-0.803, 0.243, -1.656, 0.656, 1.143, -0.453, 0.43, 0.251, -0.394, -0.862, -2.033, 1.41];
pub const G_BT: [f64; 8] = [-1.103, -0.725, -0.782, 0.267, -0.249, 0.126, 0.843, 0.858];
pub const GTA_C: [f64; 6] = [1.446536, 0.44224100000000005, 2.598184, 0.40654, 0.09404400000000011, 0.2033329999999999];
pub const GTB_C: [f64; 6] = [-0.876655, -0.9652350000000001, 2.1739550000000003, -1.8163349999999998, -1.7260630000000001, 0.5217269999999999];
pub const GTAB_C: [f64; 6] = [-0.15630499999999994, -0.09450299999999995, -1.835797, -1.4492099999999999, 2.839571, 1.2329039999999998];
pub const GS_BIGA: [f64; 48] = [0.034, 1.36, 1.225, -0.51, -0.298, -0.527, 0.57, -0.056, 0.747, -1.847, 1.567, -0.096, 0.68, -0.137, -0.379, 0.463, 0.825, -0.203, -0.153, 0.686, -0.87, -1.514, 0.395, -0.671, -1.92, -0.814, -0.468, -1.193, -1.492, 0.037, 0.897, -0.233, -0.744, 0.385, 0.717, -0.3, 0.545, 1.043, -0.207, -0.814, 0.348, 0.248, 1.099, -1.285, -0.662, -0.838, -1.734, 0.126];
pub const GS_BIGB: [f64; 48] = [-0.007, 1.046, 0.742, 0.724, 1.619, -1.206, -0.627, -1.321, -0.108, 0.999, -0.022, 0.496, -1.911, 0.147, -0.907, 1.775, 0.887, 0.949, -0.058, 0.613, 0.658, -0.344, -0.497, -0.115, -0.605, -0.594, -0.283, -0.728, 0.766, -1.596, 0.824, -0.626, -0.546, -1.351, -0.144, -0.248, 0.191, -0.534, 0.094, 1.82, 0.409, -0.574, 0.953, -0.129, 0.594, 0.613, -0.391, -1.93];
pub const GS_C: [f64; 6] = [0.8190290000000001, -0.7895310000000001, -0.6899489999999997, 0.028725999999999988, 1.107878, -0.337524];
pub const GAB_C0: [f64; 6] = [1.827, -3.078, 0.958, 0.07, 1.318, 0.386];
pub const GAB_C: [f64; 6] = [-10.918480000000002, -5.795602, 7.269576, -0.140074, 12.54401, 4.76801];
pub const GF32_C: [f64; 6] = [-5.915990034619568, -2.128300923400537, 3.3952880297179213, -0.08753699926707137, 5.942505051801799, 2.287504901648388];
pub const GF32_TC: [f64; 6] = [1.446536032690048, 0.4422409897506383, 2.5981841080341344, 0.406539974793829, 0.09404398485994392, 0.2033330020519048];
pub const GI_A: [i64; 12] = [-4, 4, 2, -2, -5, 2, -2, 4, -2, 2, -3, 1];
pub const GI_B: [i64; 8] = [-5, 3, -4, 2, 3, 0, 3, 1];
pub const GI_C: [i64; 6] = [4, -6, 23, -7, -4, -1];
pub const GI_AT: [i64; 12] = [4, 3, 4, -5, -5, 2, 2, 3, 3, -1, 4, 0];
pub const GI_TC: [i64; 6] = [3, 1, 26, 3, -19, 16];
pub const GBIG_SUM: f64 = 6.0;
pub const GBIG_SUMSQ: f64 = 36602140.0;
pub const GBIG_00: f64 = 4.0;
pub const GBIG_LAST: f64 = -18.0;
pub const GBIG_MID: f64 = -27.0;
// ---- det / slogdet / inv / solve ----
pub const DET_A3: [f64; 9] = [-0.36, 1.204, 1.397, 0.317, 0.414, -0.49, -0.914, -0.9, -0.998];
pub const DET_A3_V: f64 = 1.3576851359999997;
pub const DET_A4: [f64; 16] = [0.359, 1.511, -1.786, 1.687, -0.047, -0.8, -0.803, -1.083, -0.224, 0.834, 0.584, 0.638, -1.695, -1.571, 1.554, 0.969];
pub const DET_A4_V: f64 = 3.7877804443600014;
pub const DET_PIVOT: [f64; 9] = [0.0, 2.0, 1.0, 3.0, 0.0, 0.0, 0.0, 1.0, 4.0];
pub const DET_PIVOT_V: f64 = -21.0;
pub const DET_P4: [f64; 16] = [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0];
pub const DET_P4_V: f64 = -1.0;
pub const DET_P5: [f64; 25] = [0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0];
pub const DET_P5_V: f64 = -1.0;
pub const SING3: [f64; 9] = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0];
pub const SING3_DET: f64 = 0.0;
pub const DET_H8_EXACT: f64 = 2.737050113791513e-33;
pub const INV_H8_EXACT: [f64; 64] = [64.0, -2016.0, 20160.0, -92400.0, 221760.0, -288288.0, 192192.0, -51480.0, -2016.0, 84672.0, -952560.0, 4656960.0, -11642400.0, 15567552.0, -10594584.0, 2882880.0, 20160.0, -952560.0, 11430720.0, -58212000.0, 149688000.0, -204324120.0, 141261120.0, -38918880.0, -92400.0, 4656960.0, -58212000.0, 304920000.0, -800415000.0, 1109908800.0, -776936160.0, 216216000.0, 221760.0, -11642400.0, 149688000.0, -800415000.0, 2134440000.0, -2996753760.0, 2118916800.0, -594594000.0, -288288.0, 15567552.0, -204324120.0, 1109908800.0, -2996753760.0, 4249941696.0, -3030051024.0, 856215360.0, 192192.0, -10594584.0, 141261120.0, -776936160.0, 2118916800.0, -3030051024.0, 2175421248.0, -618377760.0, -51480.0, 2882880.0, -38918880.0, 216216000.0, -594594000.0, 856215360.0, -618377760.0, 176679360.0];
pub const SLD_A3_SIGN: f64 = 1.0;
pub const SLD_A3_LOG: f64 = 0.30578114363752873;
pub const SLD_A3: [f64; 9] = [-0.36, 1.204, 1.397, 0.317, 0.414, -0.49, -0.914, -0.9, -0.998];
pub const SLD_A4_SIGN: f64 = 1.0;
pub const SLD_A4_LOG: f64 = 1.3317802128098248;
pub const SLD_A4: [f64; 16] = [0.359, 1.511, -1.786, 1.687, -0.047, -0.8, -0.803, -1.083, -0.224, 0.834, 0.584, 0.638, -1.695, -1.571, 1.554, 0.969];
pub const SLD_NEG_SIGN: f64 = -1.0;
pub const SLD_NEG_LOG: f64 = 0.6931471805599455;
pub const SLD_NEG: [f64; 4] = [1.0, 2.0, 3.0, 4.0];
pub const SLD_BIG_LOG: f64 = 909.5211117326481;
pub const SLD_BIG_SIGN: f64 = 1.0;
pub const INV_A4: [f64; 16] = [-0.1591932734374424, -1.6697410161186446, -1.9223341772203204, -0.3233449567605392, -0.005022815413794514, 0.6530608846886208, 1.4399959723435336, -0.20947265493738576, -0.34470873092556276, -0.613120637300401, -0.05532801150395351, -0.04869633620783043, 0.2662063223599468, -0.8787026853565079, -0.9392600907735891, 0.20487395650280865];
pub const INV_PIVOT: [f64; 9] = [0.0, 0.3333333333333333, 0.0, 0.5714285714285714, 0.0, -0.14285714285714285, -0.14285714285714285, 0.0, 0.2857142857142857];
pub const SOL_B1: [f64; 4] = [-1.398, -1.204, -1.302, -0.623];
pub const SOL_X1: [f64; 4] = [4.937243386475064, -2.523636701181903, 1.3224749415792454, 1.7810817577959928];
pub const SOL_B2: [f64; 12] = [0.553, 0.218, -0.058, -2.319, 0.431, -2.126, 0.91, 0.606, 0.83, 0.828, 0.299, -0.535];
pub const SOL_X2: [f64; 12] = [1.7670418107000134, -2.0159771650234135, 2.1365547949016337, -0.38027283197227957, 1.0903795029544912, -0.08085159011737415, 1.1405337728491127, -0.38749047751578264, 1.3036178716172433, 1.4998325769871519, -0.8286221811286417, 0.9634885002999777];
pub const SOL_BP: [f64; 3] = [1.0, 2.0, 3.0];
pub const SOL_XP: [f64; 3] = [0.6666666666666666, 0.14285714285714285, 0.7142857142857143];
// ---- cholesky ----
pub const CH_A3: [f64; 9] = [4.474, -2.838, -0.582, -2.838, 6.206, 0.096, -0.582, 0.096, 3.298];
pub const CH_L3: [f64; 9] = [2.1151832071950647, 0.0, 0.0, -1.3417277474339728, 2.0989918179368297, 0.0, -0.27515347040400706, -0.13014869505891052, 1.790349654368971];
pub const CH_A5: [f64; 25] = [3.129, -1.579, -1.048, 0.272, 1.556, -1.579, 3.982, 1.58, 0.136, -1.586, -1.048, 1.58, 6.593, -2.468, -1.88, 0.272, 0.136, -2.468, 5.17, 1.285, 1.556, -1.586, -1.88, 1.285, 3.701];
pub const CH_L5: [f64; 25] = [1.768897962009115, 0.0, 0.0, 0.0, 0.0, -0.8926461751397866, 1.7847080450337836, 0.0, 0.0, 0.0, -0.5924592726703587, 0.5889722407256852, 2.4279834657347803, 0.0, 0.0, 0.15376805550223052, 0.15311213918889227, -1.0161014373427733, 2.022486075840824, 0.0, 0.8796437292701127, -0.4486948953974024, -0.45081776439196214, 0.3759546036915961, 1.5431530318383238];
// ---- qr ----
pub const QR_TALL: [f64; 15] = [-1.142, -1.069, -0.757, 0.769, -1.038, -1.026, -0.393, 1.388, -0.917, -0.697, 0.193, 0.145, 0.386, -0.571, -0.936];
pub const QR_TALL_Q: [f64; 15] = [-0.696966385707227, -0.6048071212084922, -0.2486740516831626, 0.4693232492196652, -0.42838375371080106, -0.44725697591025415, -0.23984920278716307, 0.6265805296577621, -0.7238291940721313, -0.42538141054110085, 0.032776920969056045, 0.06863798076105843, 0.23557707958230267, -0.23878488016984314, -0.4576926249406118];
pub const QR_TALL_R: [f64; 9] = [1.638529523688847, -0.2916242845134963, -0.016160831780671958, 0.0, 2.103567036412672, 0.551041677745429, 0.0, 0.0, 1.7492360895269854];
pub const QR_WIDE: [f64; 15] = [-1.232, 0.267, -0.007, 0.502, -1.327, 1.108, 0.094, -1.171, -1.358, -1.307, -0.718, 1.186, 0.895, -0.544, -0.547];
pub const QR_WIDE_Q: [f64; 9] = [-0.6822358067403569, 0.13127540634977586, -0.7192503539721312, 0.6135692158022041, -0.43217953939706105, -0.6608734093201304, -0.3976017120451106, -0.8921813792129073, 0.21430180858908193];
pub const QR_WIDE_R: [f64; 15] = [1.8058272342613508, -0.5960370845997693, -1.0695674333375724, -0.9594140386905121, 0.3208800869796482, 0.0, -1.0637014589544416, -0.293339021606042, 1.138146738780618, 0.8786794081952669, 0.0, 0.0, 0.9707176334789059, 0.41982222829026683, 1.700983676404201];
pub const QR_SQ: [f64; 16] = [0.305, -1.04, 0.75, 0.941, -1.951, -1.302, 0.128, -0.316, -0.017, -0.853, 0.879, 0.778, 0.066, 1.127, 0.468, -0.859];
pub const QR_SQ_Q: [f64; 16] = [-0.15436220929447608, -0.6614348081866801, -0.30719054697586046, 0.6665660287694309, 0.9874120338804017, -0.07947901949063041, -0.06954440239822912, 0.11774606959066536, 0.008603795272151117, -0.45677490765808687, -0.5454869833844752, -0.7026568219258718, -0.0334029698801161, 0.5895284248565152, -0.7766862268006625, 0.21931480339351017];
pub const QR_SQ_R: [f64; 16] = [-1.975872212467193, -1.1700579548680636, 0.0025472295061608197, -0.42189013780355566, 0.0, 1.8454014149365892, -0.63185526163342, -1.4590705794543652, 0.0, 0.0, -1.0822668062765326, -0.024305677797796993, 0.0, 0.0, 0.0, -0.14502754849196892];
pub const QR_UP: [f64; 9] = [2.0, 1.0, 3.0, 0.0, 4.0, 5.0, 0.0, 0.0, 6.0];
pub const QR_UP_Q: [f64; 9] = [1.0, 0.0, 0.0, -0.0, 1.0, 0.0, -0.0, -0.0, 1.0];
pub const QR_UP_R: [f64; 9] = [2.0, 1.0, 3.0, 0.0, 4.0, 5.0, 0.0, 0.0, 6.0];
pub const QR_ZC: [f64; 12] = [1.0, 0.0, 2.0, 2.0, 0.0, 1.0, 3.0, 0.0, 5.0, 4.0, 0.0, 2.0];
pub const QR_ZC_Q: [f64; 12] = [-0.18257418583505536, -0.3651483716701107, -0.23593668160633888, -0.3651483716701107, 0.887251611839318, -0.07285115482625663, -0.5477225575051661, -0.16912258224102303, -0.7270218130542391, -0.7302967433402214, -0.225496776321364, 0.6406761076053924];
pub const QR_ZC_R: [f64; 9] = [-5.477225575051661, 0.0, -4.929503017546495, 0.0, 0.0, -1.1396515953487465, 0.0, 0.0, -2.898481368099346];
pub const QR_TINY: [f64; 12] = [2.44e-21, 6.78e-21, -5.859999999999999e-21, -9.09e-21, -1.992e-20, 9.72e-21, 1.7e-22, 2.06e-21, -7.84e-21, 1.226e-20, 9.43e-21, -1.22e-21];
pub const QR_TINY_Q: [f64; 12] = [-0.1578574958930603, 0.3207862758929807, -0.24518713223092567, 0.5880838678966879, -0.7163919782317487, -0.23246849385004714, -0.01099826815648371, 0.16470752841742203, -0.9346703865300681, -0.7931692211675899, -0.5972856341873108, -0.11060261343746797];
pub const QR_TINY_R: [f64; 9] = [-1.5456979006261218e-20, -2.0287146658669698e-20, 7.695112994060431e-21, 0.0, 1.1152353135084391e-20, -9.405756154229533e-21, 0.0, 0.0, 6.639953853440211e-21];
// ---- lstsq ----
pub const LS_A: [f64; 12] = [1.0, 0.0, 1.0, 1.0, 1.0, 2.0, 1.0, 3.0, 1.0, 4.0, 1.0, 5.0];
pub const LS_B: [f64; 6] = [1.646, 2.157, 3.301, 3.484, 5.2, 5.443];
pub const LS_X: [f64; 2] = [1.5172857142857155, 0.808485714285714];
pub const LS_A2: [f64; 15] = [-0.575, -0.622, 0.08, 1.252, -0.324, -1.102, -0.799, 1.777, -0.35, -1.184, -0.302, 0.298, 0.287, 1.863, -0.191];
pub const LS_B2: [f64; 10] = [-0.826, 0.243, -0.064, -0.094, 0.323, 0.666, 0.189, 0.178, -1.137, 0.578];
pub const LS_X2: [f64; 6] = [-0.20860968297904584, -0.28930329397950233, -0.1861839010861907, 0.2346627125982878, -0.2048372086366782, -0.298557673249119];
pub const LS_UA: [f64; 8] = [0.105, -0.807, 0.892, -0.929, -1.136, 0.513, 1.319, 1.276];
pub const LS_UB: [f64; 2] = [-0.971, -0.197];
pub const LS_UX: [f64; 4] = [0.054309941616970225, 0.30893366218734253, -0.50816996984814, 0.2950543222605549];
pub const LS_UB2: [f64; 4] = [0.842, -2.976, -0.305, 1.45];
pub const LS_UX2: [f64; 8] = [0.06320216910429532, -0.3096128511170806, -0.30009300260393734, 1.0858363136697258, 0.28909195156560674, -0.9032864427973852, -0.3609458539786329, 1.3579001491129916];
pub const LS_RA: [f64; 12] = [-0.85, 0.322, -0.528, 1.789, 0.879, 2.668, 0.372, 1.488, 1.8599999999999999, 0.729, 2.081, 2.81];
pub const LS_RB: [f64; 4] = [0.332, 0.461, 1.558, 1.334];
pub const LS_RX: [f64; 3] = [-0.3722610361585875, 0.5972962268397468, 0.22503519068115935];
pub const LS_SA: [f64; 9] = [-0.575, 0.092, -2.112, -1.397, -0.757, -0.222, -1.15, 0.184, -4.224];
pub const LS_SB: [f64; 3] = [1.0, 2.0, 2.0];
pub const LS_SX: [f64; 3] = [-1.0842918228605107, -0.5813001555969532, -0.2036040796260068];
pub const LS_QA: [f64; 9] = [-1.007, 1.602, 1.922, -0.016, -0.937, -0.591, -0.095, -1.089, 2.165];
pub const LS_QB: [f64; 3] = [-1.051, -0.362, 0.812];
pub const LS_QX: [f64; 3] = [2.043637199239659, 0.044272966011213125, 0.4870017523852096];
// ---- eigh ----
pub const EH_A3: [f64; 9] = [-2.102, -1.613, 1.992, -1.613, -2.552, 2.817, 1.992, 2.817, -0.354];
pub const EH_W3: [f64; 3] = [-6.037232241184393, -0.7282103682744065, 1.7574426094587985];
pub const EH_V3: [f64; 9] = [-0.5356162798199194, -0.8081183772535123, 0.24507119197698402, -0.6667523980035144, 0.5827948472963642, 0.46453353562619204, 0.5182243148950394, -0.08540991927146147, 0.8509692739115846];
pub const EH_A5: [f64; 25] = [-1.134, -0.4019999999999999, 0.04500000000000015, -0.31100000000000005, -1.271, -0.4019999999999999, 1.682, -1.455, -0.9729999999999999, 1.5489999999999997, 0.04500000000000015, -1.455, 2.358, 0.863, 2.097, -0.31100000000000005, -0.9729999999999999, 0.863, -1.1, -0.0129999999999999, -1.271, 1.5489999999999997, 2.097, -0.0129999999999999, 0.144];
pub const EH_W5: [f64; 5] = [-2.7035928534429283, -1.626856877682279, -0.7613039522039882, 3.0013784266676864, 4.040375256661511];
pub const EH_V5: [f64; 25] = [0.517773574918882, 0.3288762037821929, -0.7485053227812628, 0.24669796104853664, -0.05129196315356757, -0.2945494734213471, 0.40820919431788144, -0.23849331151837563, -0.7313952218414898, -0.39343084679520257, -0.39668626996409545, 0.07243255371487717, -0.34510279703308205, -0.13717563538972918, 0.836349468176839, 0.14072427872201826, 0.8244157079366943, 0.48683460264203854, 0.12786900250132005, 0.2172026835817012, 0.6840963480556441, -0.2007435227591171, 0.1635752159693776, -0.6074812826224969, 0.30971578040372505];
pub const EH_REP1: [f64; 9] = [3.0, 1.0, 1.0, 1.0, 3.0, 1.0, 1.0, 1.0, 3.0];
pub const EH_REP1_W: [f64; 3] = [1.9999999999999998, 2.0, 4.999999999999998];
pub const EH_ONES4: [f64; 16] = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
pub const EH_ONES4_W: [f64; 4] = [-5.660001591138239e-16, -1.2325951644078312e-32, 1.0888646801245488e-17, 3.999999999999999];
pub const EH_A6: [f64; 36] = [-2.044, -1.721, 0.88, 0.057999999999999996, -2.457, -1.079, -1.721, -4.118, -2.238, -1.008, -0.479, -0.941, 0.88, -2.238, -0.128, 3.561, -0.873, -0.41, 0.057999999999999996, -1.008, 3.561, 0.612, 1.6589999999999998, 0.672, -2.457, -0.479, -0.873, 1.6589999999999998, 1.026, 0.707, -1.079, -0.941, -0.41, 0.672, 0.707, 2.532];
pub const EH_W6: [f64; 6] = [-6.248856143901629, -3.5859979661212673, -2.620612523621403, 1.60641949041138, 3.9604469159737845, 4.76860022725913];
pub const EH_W3_BIG: [f64; 3] = [-6.037232241184392e+20, -7.282103682744082e+19, 1.7574426094587983e+20];
pub const EH_W3_TINY: [f64; 3] = [-6.037232241184392e-20, -7.282103682744069e-21, 1.7574426094587984e-20];
// ---- svd ----
pub const SV_TALL: [f64; 12] = [-0.083, -1.012, 0.906, -0.445, 1.404, -1.464, 0.463, -0.019, 0.722, -0.833, 0.272, 1.517];
pub const SV_TALL_S: [f64; 3] = [2.7092022160516342, 1.4784647494778507, 0.6662937318422774];
pub const SV_TALL_U: [f64; 12] = [-0.479046213362035, -0.16675052331557838, -0.4902922628108834, 0.7332887942105764, 0.40766215565643543, -0.07628275783861033, -0.23426506122785146, 0.0005341541857484255, 0.8663440417075494, -0.421808191997981, 0.8977781154042654, -0.05694242076121958];
pub const SV_TALL_VT: [f64; 9] = [-0.016112188625200342, 0.5182527293884169, -0.8550757310666481, -0.6190010434126909, 0.6664309683037136, 0.41558088591726, 0.7852248758360283, 0.5359886873500515, 0.3100613187759207];
pub const SV_WIDE: [f64; 15] = [0.049, 0.518, -0.11, -0.089, 0.503, 0.511, 0.345, -1.421, -0.57, -0.307, 0.74, 0.045, -0.394, 0.378, 0.326];
pub const SV_WIDE_S: [f64; 3] = [1.7479453070398756, 0.9073358115024143, 0.6616199277247706];
pub const SV_SQ: [f64; 9] = [-1.003, -0.028, 0.586, -1.687, -1.397, -0.958, 1.003, 1.113, 0.321];
pub const SV_SQ_S: [f64; 3] = [2.874996618728019, 1.028459400522748, 0.3000261718231435];
pub const SV_TALL2: [f64; 28] = [-1.093, 0.778, -0.382, 0.296, 0.336, -1.947, 0.501, -0.966, -1.244, 1.116, -1.613, -1.589, -0.841, 1.406, 0.924, -0.167, 1.807, -0.76, -1.795, 0.67, 0.217, 0.345, -1.095, 0.564, 0.307, -1.251, 0.537, 0.279];
pub const SV_TALL2_S: [f64; 4] = [3.8541228229295337, 3.0413725035123864, 2.227820448664921, 0.8728542886203338];
pub const SV_R1: [f64; 9] = [4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 12.0, 15.0, 18.0];
pub const SV_R1_S: [f64; 3] = [32.83291031876401, 4.435120956233329e-15, 0.0];
pub const SV_R2: [f64; 12] = [0.246, 0.207, 0.03900000000000001, -0.229, -0.865, 0.636, 0.042, 0.176, -0.13399999999999998, 0.814, 1.543, -0.729];
pub const SV_R2_S: [f64; 3] = [2.193042089754217, 0.3552413159620629, 1.4085179895077719e-16];
pub const SV_RW: [f64; 18] = [-1.109, 0.396, 0.406, 0.436, 0.51, 0.025, -2.218, 0.792, 0.812, 0.872, 1.02, 0.05, -0.571, 1.751, 0.98, 0.507, 1.073, -1.564];
pub const SV_RW_S: [f64; 3] = [3.8301038427837413, 1.8838764698071941, 7.407969752399326e-17];
pub const SV_ROW_S: [f64; 1] = [5.477225575051661];
pub const SV_H8_S: [f64; 8] = [1.6959389969219494, 0.29812521131693065, 0.026212843578119035, 0.0014676881177418427, 5.436943369749976e-05, 1.2943320918688685e-06, 1.798873745949091e-08, 1.1115389812031695e-10];
pub const SV_SQ_BIG_S: [f64; 3] = [2.8749966187280184e+20, 1.0284594005227477e+20, 3.000261718231431e+19];
pub const SV_SQ_TINY_S: [f64; 3] = [2.874996618728018e-20, 1.0284594005227479e-20, 3.0002617182314347e-21];
// ---- pinv / rank ----
pub const PI_TALL: [f64; 12] = [0.964, -0.411, -0.58, -1.199, -1.177, 2.553, 0.133, 0.112, 0.867, 0.291, -0.927, -0.986];
pub const PI_TALL_P: [f64; 12] = [0.8871691770974355, -0.09518808777187265, 0.6537308590433596, -0.19349762272567003, -0.2495005301375285, -0.31741503767781115, 0.1450128003565304, -0.5475904521324274, 0.27348559290475316, 0.19058838716485343, 0.4197152361108416, -0.3125318273273697];
pub const PI_WIDE: [f64; 15] = [0.566, 0.775, -0.719, 0.501, 1.231, 0.045, 1.34, 1.149, -0.248, 0.668, 0.588, 1.374, 0.321, -1.818, 0.405];
pub const PI_WIDE_P: [f64; 15] = [0.18761236748379106, -0.1689889271189681, 0.16298722538144064, 0.1450821361865459, 0.2643162009825666, 0.08314311566125818, -0.3433460471845425, 0.5545099425682972, -0.18258801843914965, 0.18927400464072766, 0.2708487467457606, -0.47572279921562416, 0.35717356637604414, 0.1249391580624097, -0.04031709752566317];
pub const PI_R1: [f64; 9] = [4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 12.0, 15.0, 18.0];
pub const PI_R1_P: [f64; 9] = [0.003710575139146569, 0.007421150278293135, 0.011131725417439703, 0.004638218923933211, 0.009276437847866418, 0.013914656771799627, 0.005565862708719855, 0.011131725417439705, 0.016697588126159558];
pub const PI_R2: [f64; 12] = [0.246, 0.207, 0.03900000000000001, -0.229, -0.865, 0.636, 0.042, 0.176, -0.13399999999999998, 0.814, 1.543, -0.729];
pub const PI_R2_P: [f64; 12] = [1.2168494022476162, 1.3567994807441282, -0.313106442206705, 0.8491145766300242, 0.08658958755892734, -0.1273875965669482, 0.02457120931784461, 0.3462267130081776, 1.1302598146886886, 1.484187077311076, -0.33767765152454954, 0.5028878636218465];
pub const RANK_H12: usize = 11;
// ---- norms ----
pub const NV_X: [f64; 6] = [3.0, -4.0, 2.5, 1.5, -2.0, 7.0];
pub const NV_Z: [f64; 7] = [3.0, -4.0, 0.0, 1.5, -2.0, 0.0, 7.0];
pub const NV_X_1: f64 = 20.0;
pub const NV_X_2: f64 = 9.300537618869138;
pub const NV_X_INF: f64 = 7.0;
pub const NV_X_NEGINF: f64 = 1.5;
pub const NV_X_0: f64 = 6.0;
pub const NV_X_NEG1: f64 = 0.4361370716510904;
pub const NV_X_NEG2: f64 = 0.9766147891977883;
pub const NV_X_3: f64 = 7.725032379779304;
pub const NV_Z_1: f64 = 17.5;
pub const NV_Z_2: f64 = 8.958236433584458;
pub const NV_Z_INF: f64 = 7.0;
pub const NV_Z_NEGINF: f64 = 0.0;
pub const NV_Z_0: f64 = 5.0;
pub const NV_Z_NEG1: f64 = 0.0;
pub const NM_M: [f64; 12] = [3.629, -0.956, -0.219, 1.046, -0.627, 0.014, 3.514, 2.402, 1.298, 0.007, 0.236, 3.177];
pub const NM_M_FRO: f64 = 6.750604195181347;
pub const NM_M_1: f64 = 6.625;
pub const NM_M_2: f64 = 5.050535905057741;
pub const NM_M_INF: f64 = 6.5569999999999995;
pub const NM_M_NUC: f64 = 10.861429486169532;
pub const NM_M_NEG1: f64 = 0.977;
pub const NM_M_NEG2: f64 = 1.6445935147502861;
pub const NM_M_NEGINF: f64 = 4.718;
pub const NA_M: [f64; 9] = [0.0, 3.0, -4.0, 1.0, 0.0, 2.0, -5.0, 0.0, 0.0];
pub const NA_1_0: [f64; 3] = [6.0, 3.0, 6.0];
pub const NA_1_1: [f64; 3] = [7.0, 3.0, 5.0];
pub const NA_2_0: [f64; 3] = [5.0990195135927845, 3.0, 4.47213595499958];
pub const NA_2_1: [f64; 3] = [5.0, 2.23606797749979, 5.0];
pub const NA_INF_0: [f64; 3] = [5.0, 3.0, 4.0];
pub const NA_INF_1: [f64; 3] = [4.0, 2.0, 5.0];
pub const NA_NEGINF_0: [f64; 3] = [0.0, 0.0, 0.0];
pub const NA_NEGINF_1: [f64; 3] = [0.0, 0.0, 0.0];
pub const NA_0_0: [f64; 3] = [2.0, 1.0, 2.0];
pub const NA_0_1: [f64; 3] = [2.0, 2.0, 1.0];
// ---- hilbert ----
pub const COND_H8: f64 = 15257575538.072489;
// ---- eig / eigvals ----
pub const E3: [f64; 9] = [1.5, 0.5, -0.5, -0.5, 2.5, 0.5, -1.0, 1.0, 2.0];
pub const E3_RE: [f64; 3] = [2.9999999999999987, 1.000000000000001, 2.0000000000000004];
pub const E3_IM: [f64; 3] = [0.0, 0.0, 0.0];
pub const E4: [f64; 16] = [0.343, -1.611, -0.818, 0.562, 1.528, -0.283, 0.778, 1.238, -0.124, 0.69, 1.363, 1.731, 0.227, -0.683, -1.461, -0.01];
pub const E4_RE: [f64; 4] = [-0.15265792167111872, -0.15265792167111872, 0.8591579216711198, 0.8591579216711198];
pub const E4_IM: [f64; 4] = [1.000138950345139, -1.000138950345139, 1.950308271894911, -1.950308271894911];
pub const E6: [f64; 36] = [-1.158, 0.29, 0.781, 0.544, -0.961, 1.071, 0.701, 0.705, 0.745, 1.104, 2.243, -0.611, 0.047, 1.754, -1.338, 0.326, -0.689, -0.02, 0.475, -1.931, -0.992, -1.405, -0.231, -0.689, 1.515, -0.603, 1.714, -0.406, 0.271, 0.04, 0.012, -1.127, 0.335, 0.384, 0.238, 0.621];
pub const E6_RE: [f64; 6] = [-1.1922838834404534, -1.1922838834404534, -1.8217834333519922, 0.870202406612903, 0.870202406612903, 0.16194638700709704];
pub const E6_IM: [f64; 6] = [2.2031351343664127, -2.2031351343664127, 0.0, 0.2699618900090247, -0.2699618900090247, 0.0];
pub const E5: [f64; 25] = [-1.158, 0.29, 0.781, 0.544, -0.961, 1.071, 0.701, 0.705, 0.745, 1.104, 2.243, -0.611, 0.047, 1.754, -1.338, 0.326, -0.689, -0.02, 0.475, -1.931, -0.992, -1.405, -0.231, -0.689, 1.515];
pub const E5_RE: [f64; 5] = [3.1998899911623173, -2.1129331591730427, 0.3667263420463849, 0.3667263420463849, -0.24040951608204128];
pub const E5_IM: [f64; 5] = [0.0, 0.0, 1.742199708127341, -1.742199708127341, 0.0];
pub const E_ROT: [f64; 4] = [0.0, -1.0, 1.0, 0.0];
pub const E_ROT_RE: [f64; 2] = [0.0, 0.0];
pub const E_ROT_IM: [f64; 2] = [1.0, -1.0];
pub const E_PERM4: [f64; 16] = [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0];
pub const E_PERM4_RE: [f64; 4] = [-1.0000000000000004, 8.326672684688674e-17, 8.326672684688674e-17, 0.9999999999999999];
pub const E_PERM4_IM: [f64; 4] = [0.0, 0.9999999999999996, -0.9999999999999996, 0.0];
pub const E_UPPER4: [f64; 16] = [4.0, 1.0, 2.0, 3.0, 0.0, -2.0, 1.0, 5.0, 0.0, 0.0, 7.0, 1.0, 0.0, 0.0, 0.0, 0.5];
pub const E_UPPER4_RE: [f64; 4] = [4.0, -2.0, 7.0, 0.5];
pub const E_UPPER4_IM: [f64; 4] = [0.0, 0.0, 0.0, 0.0];
pub const E_COMP6: [f64; 36] = [21.0, -175.0, 735.0, -1624.0, 1764.0, -720.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0];
pub const E30_RE: [f64; 30] = [-2.690424146697132, -2.690424146697132, -2.57327375865374, -2.57327375865374, -2.1701838880148108, -1.2362802950407918, -1.2362802950407918, -0.7126182008315689, -0.7126182008315689, -0.6638054058772218, -0.6638054058772218, -0.659719903982807, -0.659719903982807, -0.16406018847527348, -0.16406018847527348, 0.48485851954562886, 0.5824094184721136, 0.5824094184721136, 0.9534394309150951, 0.9534394309150951, 1.014218127457717, 1.014218127457717, 1.2115627618802736, 1.2115627618802736, 2.1075567426265254, 2.1075567426265254, 2.398028792895652, 2.398028792895652, 2.543592631137007, 3.3356659879544943];
pub const E30_IM: [f64; 30] = [-1.1265346720733735, 1.1265346720733735, -1.5903669594089764, 1.5903669594089764, 0.0, -1.1089778269757893, 1.1089778269757893, -0.4989061241034241, 0.4989061241034241, -1.1469614623864803, 1.1469614623864803, -2.783658301347655, 2.783658301347655, -2.0357746850655865, 2.0357746850655865, 0.0, -2.6262535615898144, 2.6262535615898144, -2.891666217915717, 2.891666217915717, -0.6461339102410915, 0.6461339102410915, -1.166710496631389, 1.166710496631389, -1.0528186935469983, 1.0528186935469983, -1.3403183707010584, 1.3403183707010584, 0.0, 0.0];
pub const E30_TRACE: f64 = 3.3280000000000003;
pub const E30_DET: f64 = -999906183.1715429;
pub const E_SYM5: [f64; 25] = [-1.134, -0.4019999999999999, 0.04500000000000015, -0.31100000000000005, -1.271, -0.4019999999999999, 1.682, -1.455, -0.9729999999999999, 1.5489999999999997, 0.04500000000000015, -1.455, 2.358, 0.863, 2.097, -0.31100000000000005, -0.9729999999999999, 0.863, -1.1, -0.0129999999999999, -1.271, 1.5489999999999997, 2.097, -0.0129999999999999, 0.144];
pub const E_SYM5_W: [f64; 5] = [-2.7035928534429297, -1.626856877682279, -0.7613039522039883, 3.001378426667687, 4.0403752566615125];
pub const E_C4_SCALED_RE: [f64; 4] = [-1.5265792167111828e-21, -1.5265792167111828e-21, 8.591579216711172e-21, 8.591579216711172e-21];
pub const E_C4_SCALED_IM: [f64; 4] = [1.0001389503451393e-20, -1.0001389503451393e-20, 1.950308271894906e-20, -1.950308271894906e-20];
// ---- slogdet ----
// ---- cond ----
pub const COND_A4_2: f64 = 11.652546104148914;
pub const COND_A4_1: f64 = 20.59515257645428;
pub const COND_A4_INF: f64 = 23.587937108855385;
pub const COND_A4_FRO: f64 = 15.524783679737261;
pub const COND_A4_NEG2: f64 = 0.08581815433829931;
pub const COND_A3_2: f64 = 4.037776860743056;
pub const COND_A3_1: f64 = 6.621956027660305;
pub const COND_A3_INF: f64 = 5.816100569727385;
pub const COND_A3_FRO: f64 = 5.300909673305572;
pub const COND_A3_NEG2: f64 = 0.24766103588398242;
pub const COND_TALL: f64 = 4.06607789714724;
pub const COND_WIDE: f64 = 2.641917562928982;
pub const COND_SING_2: f64 = 2.517588727560788e+16;
pub const COND_SING_1: f64 = f64::INFINITY;
// ---- multi_dot ----
pub const MD3_0: [f64; 48] = [0.252, 0.136, 1.239, -0.032, -1.318, -0.009, 0.52, 1.208, 1.876, -1.355, -0.571, -0.159, 0.895, -0.181, -0.608, -0.223, -1.206, -0.626, 0.39, -0.081, -1.36, -1.776, 1.274, 1.341, -0.109, -0.217, 0.817, 0.876, 1.307, 1.747, 0.257, -2.115, 0.691, -0.318, -1.734, 0.729, 1.707, -0.515, 0.558, 1.009, 0.336, 0.623, -0.27, -0.234, -0.067, 0.441, 0.12, -0.968];
pub const MD3_1: [f64; 24] = [-0.307, -0.839, 0.126, -0.661, 0.158, 2.439, 0.44, 0.77, -0.794, -0.906, 0.79, -0.96, -0.932, 0.415, -0.464, 0.622, -0.898, 0.889, -2.256, -0.893, -1.347, 0.182, -1.505, 1.003];
pub const MD3_2: [f64; 12] = [-1.33, 0.938, -0.577, 0.215, -0.945, 2.22, -0.968, 0.175, -0.325, 2.175, 0.301, -1.014];
pub const MD3_OUT: [f64; 24] = [-10.561035488, 3.6506200989999993, -3.871477621, 16.806282494999998, -0.10749029899999918, -2.0430920940000017, -2.940765676, 1.326180863, -1.135942027, 3.448659064999999, -0.6330074730000002, 1.0350187620000006, -4.051736548000001, 2.7835155740000004, -1.7439388960000002, 0.9492921200000001, -2.7346940040000005, 6.3796142760000025, 0.28516338999999946, 0.5687530660000001, -0.0202712890000001, -3.1069797449999985, -1.2967842649999997, 3.5117801399999995];
pub const MD5_0: [f64; 24] = [-0.631, 0.929, 0.034, -0.716, 1.755, -0.949, -0.419, 1.446, -0.775, 0.152, -0.572, 1.179, 0.355, -0.869, -0.617, 1.286, 0.082, 1.495, -1.352, 1.119, -0.938, -0.276, -0.213, -0.338];
pub const MD5_1: [f64; 16] = [-0.245, -0.578, 0.26, -1.076, 1.469, 0.592, 0.967, -0.256, 2.37, -0.37, -0.588, 0.785, -1.583, -0.869, 0.687, -0.944];
pub const MD5_2: [f64; 18] = [-1.796, -0.47, 0.022, 1.04, -1.018, -0.101, -0.204, -0.172, -0.648, 0.279, 0.214, 0.213, 0.134, -0.17, -2.17, 0.007, 1.291, -2.198];
pub const MD5_3: [f64; 36] = [-0.989, -0.368, 1.288, 0.194, 0.92, 0.577, -0.636, 0.542, -0.317, -0.322, 0.097, -1.526, 1.192, -0.671, 1.0, 0.136, 1.532, -0.66, -0.312, 0.338, -2.207, 0.828, 1.542, 1.127, 0.755, -0.146, 1.282, 1.074, 0.393, 0.005, -0.362, -1.23, 1.226, -2.172, -0.37, 0.164];
pub const MD5_4: [f64; 20] = [-0.301, -0.578, -1.111, -1.185, 0.811, 0.608, -1.028, 1.084, -0.498, 1.119, -1.221, 0.24, -1.917, 0.072, 1.003, -0.085, -0.775, 0.946, 0.986, 0.475];
pub const MD5_OUT: [f64; 15] = [-0.7401630500931674, -3.47501740425973, 9.41362883567508, 11.726874343203226, 3.642232079159037, -1.0105029034203155, -1.943735628619407, 5.6075661822957565, 8.249078523697102, 2.165883098511454, -14.693802437194211, 6.196299411412308, -7.6015298695182265, 24.451752895435327, -3.041604346626336];
pub const MDI_0: [i64; 10] = [2, 0, -1, -3, 2, -1, 1, -1, 3, 2];
pub const MDI_1: [i64; 15] = [3, 0, 2, -3, 3, 0, 3, 2, 3, -3, -2, -1, -3, 2, -2];
pub const MDI_2: [i64; 12] = [1, -1, -2, 0, 3, -2, 0, 0, 0, -2, -1, -3];
pub const MDI_3: [i64; 8] = [2, -1, 2, -3, -3, -2, 0, 0];
pub const MDI_OUT: [i64; 4] = [52, 60, -134, -243];
// ---- tensordot ----
pub const TD1_A: [f64; 24] = [-1.176, 1.8, -1.182, 2.473, -1.669, 0.432, -0.855, -1.39, -0.844, -2.786, 1.193, 0.498, -0.321, -0.52, 1.646, -1.617, 0.776, -0.208, -0.158, -0.059, -0.425, -0.617, -0.139, -1.362];
pub const TD1_B: [f64; 24] = [-1.152, 1.776, -0.022, -0.535, -0.552, -0.678, 0.177, 1.055, -0.524, -0.468, 0.585, 0.049, 1.233, -1.511, 1.094, -0.062, -0.403, -1.944, 0.205, 0.246, 1.109, 1.243, -0.562, 2.114];
pub const TD1_OUT: [f64; 4] = [-3.868196, 0.39006099999999977, 2.524538, -6.7360760000000015];
pub const TD2_A: [f64; 12] = [-0.311, -2.245, -0.215, 0.887, 0.174, -0.196, -1.298, 0.058, 1.243, -0.442, -0.316, -0.584];
pub const TD2_B: [f64; 8] = [-1.12, 1.404, 0.052, 0.495, -1.597, -0.075, -0.136, 1.605];
pub const TD2_OUT: [f64; 6] = [0.45430299999999996, -0.1081590000000001, 1.859946, 0.33771599999999996, -0.8310680000000001, 0.612762];
pub const TD3_A: [f64; 6] = [-0.563, -0.618, -0.218, 0.709, 0.169, 0.282];
pub const TD3_B: [f64; 4] = [-1.424, 0.194, -2.108, -0.745];
pub const TD3_OUT: [f64; 24] = [0.8017119999999999, -0.10922199999999999, 1.186804, 0.41943499999999995, 0.8800319999999999, -0.119892, 1.3027440000000001, 0.46041, 0.310432, -0.042292, 0.459544, 0.16241, -1.0096159999999998, 0.137546, -1.494572, -0.5282049999999999, -0.240656, 0.032786, -0.356252, -0.12590500000000002, -0.4015679999999999, 0.054708, -0.594456, -0.21008999999999997];
pub const TD4_B: [f64; 6] = [-0.514, -0.049, -0.374, -1.218, 0.682, -0.441];
pub const TD4_OUT: f64 = -0.47146999999999994;
pub const TD5_A: [f64; 12] = [-1.566, 0.677, 0.071, -0.427, -0.368, 0.422, 0.835, -0.414, -0.333, -2.6, 0.666, 2.31];
pub const TD5_B: [f64; 15] = [0.189, -0.81, 2.336, 0.559, 0.142, -1.224, -1.841, -2.763, -1.042, -0.339, -1.025, 2.187, 2.187, -0.685, 0.267];
pub const TD5_OUT: [f64; 20] = [0.4957829999999999, 1.2176770000000001, -3.369663, -0.26383300000000004, -0.186531, 2.2764249999999997, -7.0114719999999995, -5.270714, 1.7197190000000002, -0.7411240000000001, -1.691271, -0.13820299999999988, -0.6847069999999997, -1.286591, -0.09516100000000001, -1.941717, 6.160013999999999, 5.198379999999999, -1.389655, 0.696482];
pub const TD6_A: [f64; 24] = [1.673, -0.089, -0.73, 1.008, -0.471, -1.001, 1.066, 0.848, 0.945, 0.848, 0.649, 0.407, -1.141, 0.832, 1.011, 1.175, -0.659, 0.032, -1.216, 0.055, -0.942, 1.039, 0.541, -1.477];
pub const TD6_B: [f64; 40] = [0.068, 0.892, 0.423, 1.616, -0.796, -1.983, 0.756, 1.147, 0.709, 0.914, -1.89, -1.897, 0.549, 0.928, -0.037, 0.227, -0.566, -1.341, -0.325, 0.912, -0.955, 0.294, -0.693, 0.653, 0.695, -0.871, -0.885, -0.518, 0.14, 0.327, 1.265, 1.242, 0.919, 0.734, 0.546, 0.191, -0.356, 0.222, 1.787, -1.23];
pub const TD6_OUT: [f64; 15] = [4.049555, 0.05192200000000025, -0.5962270000000003, 4.046054000000001, -2.68414, 4.298254, 3.385649, -0.864895, -0.9211199999999999, 0.577407, -0.2928789999999999, -1.3226149999999997, -2.2923419999999997, -0.5326129999999998, 1.9698780000000002];
pub const TD7_A: [i64; 12] = [2, 0, 4, -1, 4, 2, -2, 2, 1, 0, -4, 0];
pub const TD7_B: [i64; 12] = [-3, -4, -3, 4, -4, -1, 2, 0, -2, -4, -4, 4];
pub const TD7_OUT: [i64; 9] = [-22, -8, 4, -10, -20, -11, -24, 0, -12];
pub const PA_IN: [f64; 48] = [-0.968, -0.048, -0.783, 0.143, -0.629, 0.415, -0.258, -2.73, 0.265, 0.511, 0.339, 0.157, 0.347, 1.211, -0.465, -0.825, -0.716, -0.389, -0.708, 0.722, -1.267, 0.538, -0.539, 0.62, 0.106, 0.356, 1.298, 0.303, 0.053, -0.764, -0.965, -1.05, -1.915, -1.247, 0.204, -0.484, -0.67, -1.067, 0.345, 1.06, -0.134, 1.985, -0.49, -1.187, 0.634, -1.808, -0.689, -1.63];
pub const PA_OUT: [f64; 48] = [-0.968, 0.265, -0.716, -0.048, 0.511, -0.389, 0.106, -1.915, -0.134, 0.356, -1.247, 1.985, -0.783, 0.339, -0.708, 0.143, 0.157, 0.722, 1.298, 0.204, -0.49, 0.303, -0.484, -1.187, -0.629, 0.347, -1.267, 0.415, 1.211, 0.538, 0.053, -0.67, 0.634, -0.764, -1.067, -1.808, -0.258, -0.465, -0.539, -2.73, -0.825, 0.62, -0.965, 0.345, -0.689, -1.05, 1.06, -1.63];
pub const PA_OUT2: [f64; 24] = [-1.176, -1.669, -0.844, -0.321, 0.776, -0.425, 1.8, 0.432, -2.786, -0.52, -0.208, -0.617, -1.182, -0.855, 1.193, 1.646, -0.158, -0.139, 2.473, -1.39, 0.498, -1.617, -0.059, -1.362];
pub const TD8_A: [f64; 24] = [1.458, 0.444, 0.655, 0.737, 0.063, -1.074, 0.126, 1.12, -0.621, 0.14, 0.576, -0.742, 0.884, -0.726, -0.066, -1.061, 0.117, -1.456, 0.159, -1.48, -0.532, 1.839, -1.374, -0.305];
pub const TD8_B: [f64; 16] = [0.899, -0.401, 0.182, -1.515, 0.46, -0.679, -0.327, 0.007, -0.364, 0.799, 0.567, -0.448, -1.216, 1.211, 0.005, 0.464];
pub const TD8_OUT: [f64; 24] = [0.38037, 0.529718, 0.49523799999999996, -2.1572339999999994, -1.8451870000000001, 2.1609770000000004, 0.43970600000000004, 0.36026900000000006, 0.19872900000000004, -0.28437700000000005, 0.16407999999999995, 0.33945899999999996, 1.774956, -1.199135, 0.355563, -1.807078, 1.1772269999999998, -0.7235320000000001, 0.5801590000000001, -0.945399, 1.238688, -2.50253, -1.47876, 1.292885];
// ---- vdot ----
pub const VD_A: [f64; 6] = [-1.375, 0.395, -0.95, 0.18, -0.147, -2.289];
pub const VD_B: [f64; 6] = [1.924, 1.245, 1.355, -0.27, -1.172, -0.802];
pub const VD_OUT: f64 = -1.4815129999999999;
// ---- matrix_power ----
pub const MP_A: [f64; 9] = [1.637, 1.824, -0.763, 2.802, 3.509, 0.169, 0.958, 0.441, 2.308];
pub const MP_P1: [f64; 9] = [1.637, 1.824, -0.763, 2.802, 3.509, 0.169, 0.958, 0.441, 2.308];
pub const MP_P2: [f64; 9] = [7.0596630000000005, 9.049821000000001, -2.7017789999999997, 14.580993999999999, 17.498457999999996, -1.1548530000000001, 5.0149919999999995, 4.312689, 4.670438999999999];
pub const MP_P5: [f64; 9] = [825.1297028114388, 1027.2700299307783, -190.04717065268397, 1728.177378182794, 2133.9106358255563, -345.6033537982901, 597.4444343704408, 718.1883722741368, -61.45942022713488];
pub const MP_P10: [f64; 9] = [2342601.229098931, 2903243.7890916495, -500161.3441172296, 4907267.774136261, 6080671.118584082, -1044681.3118037707, 1697407.4582621802, 2102147.1271258853, -357973.67415651644];
pub const MP_PM1: [f64; 9] = [2.4635022842082535, -1.3957402395638911, 0.916608467650431, -1.9357168945648127, 1.384342588433282, -0.7412937123042361, -0.6526793924473241, 0.31482845234104434, 0.19445390602991994];
pub const MP_PM3: [f64; 9] = [27.704526885597836, -17.348974489375376, 11.933054296633664, -23.722875833755946, 14.863630535900937, -10.20673926889458, -7.982391897399852, 4.9711774016358135, -3.346842931965031];
pub const MP_PM7: [f64; 9] = [3679.5184756351923, -2305.756658386104, 1587.9044612355299, -3151.6568364545387, 1974.9739216640749, -1360.1039742899306, -1060.4251464729405, 664.5102629758849, -457.62525611660845];
pub const MP_FIB30: [f64; 4] = [1346269.0, 832040.0, 832040.0, 514229.0];
pub const MP_FIBM5: [f64; 4] = [-3.0, 5.0, 5.0, -8.0];
}
// ---------------------------------------------------------------------------
// GEMM
// ---------------------------------------------------------------------------
#[test]
fn gemm_f64_square_matches_numpy() {
let mut c = [f64::NAN; 9];
gemm_f64(3, 3, 3, 1.0, &G33_A, 3, 1, &G33_B, 3, 1, 0.0, &mut c, 3, 1);
assert_close(&c, &G33_C, 1e-14);
assert_close(&gemm(3, 3, 3, &G33_A, 3, 1, &G33_B, 3, 1), &G33_C, 1e-14);
}
#[test]
fn gemm_f64_non_square() {
// 3x4 @ 4x2; note gemm_f64 takes (m, k, n) while the generic gemm takes (m, n, k).
let mut c = [0.0; 6];
gemm_f64(3, 4, 2, 1.0, &G34_A, 4, 1, &G42_B, 2, 1, 0.0, &mut c, 2, 1);
assert_close(&c, &G32_C, 1e-14);
assert_close(&gemm(3, 2, 4, &G34_A, 4, 1, &G42_B, 2, 1), &G32_C, 1e-14);
}
#[test]
fn gemm_f64_row_times_matrix_is_1xn() {
let mut c = [0.0; 3];
gemm_f64(1, 5, 3, 1.0, &G15_A, 5, 1, &G53_B, 3, 1, 0.0, &mut c, 3, 1);
assert_close(&c, &G13_C, 1e-14);
assert_close(&gemm(1, 3, 5, &G15_A, 5, 1, &G53_B, 3, 1), &G13_C, 1e-14);
}
#[test]
fn gemm_f64_column_times_row_is_outer_product() {
let mut c = [0.0; 20];
gemm_f64(5, 1, 4, 1.0, &G51_A, 1, 1, &G14_B, 4, 1, 0.0, &mut c, 4, 1);
assert_close(&c, &G54_C, 1e-14);
assert_close(&gemm(5, 4, 1, &G51_A, 1, 1, &G14_B, 4, 1), &G54_C, 1e-14);
}
#[test]
fn gemm_f64_transposed_a_via_strides() {
// G_AT is a row-major 4x3 buffer; its transpose (3x4) has strides (1, 3).
let mut c = [0.0; 6];
gemm_f64(3, 4, 2, 1.0, &G_AT, 1, 3, &G42_B, 2, 1, 0.0, &mut c, 2, 1);
assert_close(&c, >A_C, 1e-14);
assert_close(&gemm(3, 2, 4, &G_AT, 1, 3, &G42_B, 2, 1), >A_C, 1e-14);
}
#[test]
fn gemm_f64_transposed_b_via_strides() {
// G_BT is 2x4 row-major; its transpose (4x2) has strides (1, 4).
let mut c = [0.0; 6];
gemm_f64(3, 4, 2, 1.0, &G34_A, 4, 1, &G_BT, 1, 4, 0.0, &mut c, 2, 1);
assert_close(&c, >B_C, 1e-14);
assert_close(&gemm(3, 2, 4, &G34_A, 4, 1, &G_BT, 1, 4), >B_C, 1e-14);
}
#[test]
fn gemm_f64_both_operands_transposed() {
let mut c = [0.0; 6];
gemm_f64(3, 4, 2, 1.0, &G_AT, 1, 3, &G_BT, 1, 4, 0.0, &mut c, 2, 1);
assert_close(&c, >AB_C, 1e-14);
}
#[test]
fn gemm_f64_strided_views_with_offset() {
// A = big[::2, ::2] (3x4, strides 16 and 2); B = big2[1::2, 1::3] (4x2, offset 7, strides 12 and 3).
let mut c = [0.0; 6];
gemm_f64(
3,
4,
2,
1.0,
&GS_BIGA,
16,
2,
&GS_BIGB[7..],
12,
3,
0.0,
&mut c,
2,
1,
);
assert_close(&c, &GS_C, 1e-14);
assert_close(
&gemm(3, 2, 4, &GS_BIGA, 16, 2, &GS_BIGB[7..], 12, 3),
&GS_C,
1e-14,
);
}
#[test]
fn gemm_f64_alpha_beta_accumulate_into_c() {
let mut c = GAB_C0;
gemm_f64(3, 4, 2, 2.0, &G34_A, 4, 1, &G42_B, 2, 1, 0.5, &mut c, 2, 1);
assert_close(&c, &GAB_C, 1e-14);
}
#[test]
fn gemm_f64_column_major_output() {
// C stored column-major: element (i, j) at i + 3 j.
let mut c = [0.0; 6];
gemm_f64(3, 4, 2, 1.0, &G34_A, 4, 1, &G42_B, 2, 1, 0.0, &mut c, 1, 3);
assert_close(&c, &transpose(&G32_C, 3, 2), 1e-14);
}
#[test]
fn gemm_k_zero_gives_zeros_and_beta_scaling() {
let empty: [f64; 0] = [];
// beta = 0 must overwrite even NaN garbage (BLAS semantics: C is not read).
let mut c = [f64::NAN; 6];
gemm_f64(2, 0, 3, 1.0, &empty, 0, 1, &empty, 3, 1, 0.0, &mut c, 3, 1);
assert_eq!(c, [0.0; 6]);
// beta = 0.5 scales the existing C.
let mut c = [2.0, 4.0, 6.0, 8.0, 10.0, 12.0];
gemm_f64(2, 0, 3, 1.0, &empty, 0, 1, &empty, 3, 1, 0.5, &mut c, 3, 1);
assert_eq!(c, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]);
// Generic kernel: empty inner dimension gives an m x n zero matrix.
assert_eq!(gemm(2, 3, 0, &empty, 0, 1, &empty, 3, 1), vec![0.0; 6]);
// f32 as well.
let mut c32 = [f32::NAN; 4];
gemm_f32(2, 0, 2, 1.0, &[], 0, 1, &[], 2, 1, 0.0, &mut c32, 2, 1);
assert_eq!(c32, [0.0; 4]);
}
#[test]
fn gemm_with_empty_output_dimension_is_a_noop() {
let mut c: [f64; 0] = [];
gemm_f64(0, 3, 2, 1.0, &[], 3, 1, &G42_B, 2, 1, 0.0, &mut c, 2, 1);
assert!(gemm(0, 2, 3, &[], 3, 1, &G42_B, 2, 1).is_empty());
assert!(gemm(3, 0, 4, &G34_A, 4, 1, &[], 0, 1).is_empty());
}
#[test]
fn gemm_f32_non_square_and_transposed() {
let (a, b, at) = (to_f32(&G34_A), to_f32(&G42_B), to_f32(&G_AT));
let mut c = [0.0_f32; 6];
gemm_f32(3, 4, 2, 1.0, &a, 4, 1, &b, 2, 1, 0.0, &mut c, 2, 1);
assert_close_f32(&c, &GF32_C, 2e-6);
gemm_f32(3, 4, 2, 1.0, &at, 1, 3, &b, 2, 1, 0.0, &mut c, 2, 1);
assert_close_f32(&c, &GF32_TC, 2e-6);
assert_close_f32(&gemm(3, 2, 4, &a, 4, 1, &b, 2, 1), &GF32_C, 2e-6);
}
#[test]
fn gemm_generic_i64_matches_numpy_exactly() {
assert_eq!(gemm(3, 2, 4, &GI_A, 4, 1, &GI_B, 2, 1), GI_C.to_vec());
// Transposed A through strides.
assert_eq!(gemm(3, 2, 4, &GI_AT, 1, 3, &GI_B, 2, 1), GI_TC.to_vec());
}
#[test]
fn gemm_generic_i32_and_u8_use_the_same_loop() {
let a: Vec<i32> = GI_A.iter().map(|&x| x as i32).collect();
let b: Vec<i32> = GI_B.iter().map(|&x| x as i32).collect();
let expected: Vec<i32> = GI_C.iter().map(|&x| x as i32).collect();
assert_eq!(gemm(3, 2, 4, &a, 4, 1, &b, 2, 1), expected);
// Small non-negative values fit u8: [[1,2],[3,4]] @ [[5,6],[7,8]].
let product: Vec<u8> = gemm(2, 2, 2, &[1u8, 2, 3, 4], 2, 1, &[5u8, 6, 7, 8], 2, 1);
assert_eq!(product, vec![19, 22, 43, 50]);
}
#[test]
fn gemm_large_non_square_matches_numpy_checksums() {
// 150x130 @ 130x170 with integer entries (exact in f64); also a regression test for the
// output-aliasing guard, which used to be O((m n)^2) and made this call take minutes.
let (m, k, n) = (150, 130, 170);
let a: Vec<f64> = (0..m * k)
.map(|x| (((x / k) * 7 + (x % k) * 3) % 11) as f64 - 5.0)
.collect();
let b: Vec<f64> = (0..k * n)
.map(|x| (((x / n) * 5 + (x % n) * 2) % 13) as f64 - 6.0)
.collect();
let mut c = vec![0.0; m * n];
gemm_f64(
m, k, n, 1.0, &a, k as isize, 1, &b, n as isize, 1, 0.0, &mut c, n as isize, 1,
);
assert_eq!(c.iter().sum::<f64>(), GBIG_SUM);
assert_eq!(c.iter().map(|x| x * x).sum::<f64>(), GBIG_SUMSQ);
assert_eq!(c[0], GBIG_00);
assert_eq!(c[m * n - 1], GBIG_LAST);
assert_eq!(c[77 * n + 3], GBIG_MID);
// Generic kernel agrees on the same inputs.
let g = gemm(m, n, k, &a, k as isize, 1, &b, n as isize, 1);
assert_eq!(g, c);
}
#[test]
fn gemm_f64_output_guard_is_fast_for_large_outputs_and_accepts_exotic_layouts() {
// 1200 x 1200 output with inner dimension 2. The aliasing guard used to compare every pair of
// output elements (~1e12 comparisons); it must now be effectively free.
let n = 1200;
let a: Vec<f64> = (0..n).flat_map(|i| [(i % 7) as f64, 1.0]).collect();
let b: Vec<f64> = (0..n)
.map(|_| 1.0)
.chain((0..n).map(|j| (j % 5) as f64))
.collect();
let mut c = vec![0.0; n * n];
gemm_f64(
n, 2, n, 1.0, &a, 2, 1, &b, n as isize, 1, 0.0, &mut c, n as isize, 1,
);
for (i, j) in [(0, 0), (6, 4), (1199, 1199), (500, 777)] {
assert_eq!(c[i * n + j], (i % 7 + j % 5) as f64);
}
// Padded rows (leading dimension 5 > 3) and interleaved strides (3, 2) are valid, non-aliasing layouts.
let (a, b) = ([1.0, 2.0, 3.0, 4.0], [1.0, 0.0, 0.0, 1.0]);
let mut padded = [9.0; 10];
gemm_f64(2, 2, 2, 1.0, &a, 2, 1, &b, 2, 1, 0.0, &mut padded, 5, 1);
assert_eq!(padded, [1.0, 2.0, 9.0, 9.0, 9.0, 3.0, 4.0, 9.0, 9.0, 9.0]);
let mut interleaved = [9.0; 6];
gemm_f64(
2,
2,
2,
1.0,
&a,
2,
1,
&b,
2,
1,
0.0,
&mut interleaved,
3,
2,
);
assert_eq!(interleaved, [1.0, 9.0, 2.0, 3.0, 9.0, 4.0]);
}
#[test]
#[should_panic(expected = "must not alias")]
fn gemm_f64_rejects_interleaved_output_strides_that_collide() {
// Offsets i * 2 + j * 4: (0, 1) and (2, 0) both hit element 4.
let mut c = [0.0; 12];
gemm_f64(
3, 1, 2, 1.0, &[1.0; 3], 1, 1, &[1.0; 2], 2, 1, 0.0, &mut c, 2, 4,
);
}
#[test]
#[should_panic(expected = "outside the supplied slice")]
fn gemm_rejects_strides_that_leave_the_slice() {
gemm(2, 2, 2, &[1.0; 4], 3, 1, &[1.0; 4], 2, 1);
}
#[test]
#[should_panic(expected = "outside the supplied slice")]
fn gemm_rejects_negative_strides_without_offset() {
// A negative stride would need a base pointer inside the slice; the slice API cannot express it.
gemm(2, 2, 2, &[1.0; 4], -2, 1, &[1.0; 4], 2, 1);
}
#[test]
#[should_panic(expected = "must not alias")]
fn gemm_f64_rejects_aliasing_output_strides() {
let mut c = [0.0; 6];
gemm_f64(
2, 2, 2, 1.0, &[1.0; 4], 2, 1, &[1.0; 4], 2, 1, 0.0, &mut c, 1, 1,
);
}
// ---------------------------------------------------------------------------
// determinant, slogdet-free LU facts, inverse, solve
// ---------------------------------------------------------------------------
fn det_of(a: &[f64], n: usize) -> f64 {
Lu::factor(a, n).map(|lu| lu.det()).unwrap_or(0.0)
}
#[test]
fn det_of_random_3x3_and_4x4_matches_numpy() {
assert_rel(det_of(&DET_A3, 3), DET_A3_V, 1e-12);
assert_rel(det_of(&DET_A4, 4), DET_A4_V, 1e-12);
}
#[test]
fn det_needing_row_pivots() {
// Leading entry is zero, so partial pivoting must swap rows; sign must follow.
assert_rel(det_of(&DET_PIVOT, 3), DET_PIVOT_V, 1e-14);
assert_eq!(DET_PIVOT_V.round(), DET_PIVOT_V);
}
#[test]
fn det_sign_of_permutation_matrices() {
assert_eq!(det_of(&DET_P4, 4), DET_P4_V); // 4-cycle: odd
assert_eq!(det_of(&DET_P5, 5), DET_P5_V); // (0 2 1)(3 4): even x odd
assert_eq!(DET_P4_V, -1.0);
assert_eq!(DET_P5_V, -1.0);
// A single transposition and a 3-cycle.
assert_eq!(det_of(&[0.0, 1.0, 1.0, 0.0], 2), -1.0);
assert_eq!(
det_of(&[0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0], 3),
1.0
);
}
#[test]
fn det_of_1x1_and_identity() {
assert_eq!(det_of(&[-2.5], 1), -2.5);
assert_eq!(det_of(&eye(6), 6), 1.0);
}
#[test]
fn det_of_singular_matrices_is_zero() {
// Exactly singular: LU meets an exact zero pivot.
assert!(matches!(
Lu::factor(&[1.0, 2.0, 2.0, 4.0], 2),
Err(KernelError::Singular)
));
assert!(matches!(
Lu::factor(&[0.0; 9], 3),
Err(KernelError::Singular)
));
// 1..9: NumPy returns ~0; either exact singularity or a tiny determinant is acceptable.
assert!(det_of(&SING3, 3).abs() < 1e-12);
}
#[test]
fn det_of_hilbert_8_matches_exact_value_within_condition_error() {
// Exact rational value 2.737e-33; cond(H8) ~ 1.5e10 so expect ~1e-6 relative agreement at best.
let d = det_of(&hilbert(8), 8);
assert_rel(d, DET_H8_EXACT, 1e-6);
}
#[test]
fn lu_rejects_wrong_length_input() {
assert!(matches!(
Lu::factor(&[1.0, 2.0, 3.0], 2),
Err(KernelError::Shape)
));
assert!(matches!(inv(&[1.0; 5], 2), Err(KernelError::Shape)));
}
#[test]
fn inv_matches_numpy_4x4() {
let x = inv(&DET_A4, 4).unwrap();
assert_close(&x, &INV_A4, 1e-12);
assert_close(&mm(&DET_A4, &x, 4, 4, 4), &eye(4), 1e-12);
assert_close(&mm(&x, &DET_A4, 4, 4, 4), &eye(4), 1e-12);
}
#[test]
fn inv_of_matrix_that_needs_pivoting() {
assert_close(&inv(&DET_PIVOT, 3).unwrap(), &INV_PIVOT, 1e-14);
}
#[test]
fn inv_of_permutation_is_its_transpose() {
let p = DET_P5;
assert_eq!(inv(&p, 5).unwrap(), transpose(&p, 5, 5));
}
#[test]
fn inv_of_1x1_and_singular_error() {
assert_eq!(inv(&[4.0], 1).unwrap(), vec![0.25]);
assert!(matches!(
inv(&[1.0, 2.0, 2.0, 4.0], 2),
Err(KernelError::Singular)
));
}
#[test]
fn solve_with_1d_rhs_matches_numpy() {
let lu = Lu::factor(&DET_A4, 4).unwrap();
assert_close(&lu.solve(&SOL_B1, 1).unwrap(), &SOL_X1, 1e-12);
}
#[test]
fn solve_with_2d_rhs_matches_numpy() {
let lu = Lu::factor(&DET_A4, 4).unwrap();
let x = lu.solve(&SOL_B2, 3).unwrap();
assert_close(&x, &SOL_X2, 1e-12);
assert_close(&mm(&DET_A4, &x, 4, 4, 3), &SOL_B2, 1e-12);
}
#[test]
fn solve_with_zero_leading_pivot() {
let lu = Lu::factor(&DET_PIVOT, 3).unwrap();
assert_close(&lu.solve(&SOL_BP, 1).unwrap(), &SOL_XP, 1e-14);
}
#[test]
fn solve_rejects_mismatched_rhs() {
let lu = Lu::factor(&DET_A4, 4).unwrap();
assert!(matches!(
lu.solve(&[1.0, 2.0, 3.0], 1),
Err(KernelError::Shape)
));
}
#[test]
fn hilbert_8_inverse_error_is_bounded_by_condition_number() {
// The exact inverse has integer entries up to ~4e10. Forward error is about cond * eps ~ 3e-6
// in a normwise sense (measured: 1e-8), so 1e-6 is sound but 1e-12 would be wrong.
let x = inv(&hilbert(8), 8).unwrap();
let err = max_abs(&sub(&x, &INV_H8_EXACT));
assert!(
err <= 1e-6 * max_abs(&INV_H8_EXACT),
"inverse error {err:e}"
);
assert!(err > 0.0, "Hilbert inverse cannot be exact in f64");
}
#[test]
fn hilbert_8_solve_has_small_residual_but_visible_forward_error() {
let h = hilbert(8);
let ones = vec![1.0; 8];
let b = mm(&h, &ones, 8, 8, 1);
let x = Lu::factor(&h, 8).unwrap().solve(&b, 1).unwrap();
// Backward stable: tiny residual.
assert!(max_abs(&sub(&mm(&h, &x, 8, 8, 1), &b)) < 1e-14);
// Forward error is limited by cond ~ 1.5e10 times eps (measured 4e-7).
let forward = max_abs(&sub(&x, &ones));
assert!(forward < 1e-5, "forward error {forward:e}");
}
// ---------------------------------------------------------------------------
// Cholesky
// ---------------------------------------------------------------------------
#[test]
fn cholesky_3x3_and_5x5_match_numpy() {
assert_close(&cholesky(&CH_A3, 3).unwrap(), &CH_L3, 1e-13);
assert_close(&cholesky(&CH_A5, 5).unwrap(), &CH_L5, 1e-13);
}
#[test]
fn cholesky_reconstructs_and_is_lower_triangular() {
let l = cholesky(&CH_A5, 5).unwrap();
for i in 0..5 {
for j in i + 1..5 {
assert_eq!(l[i * 5 + j], 0.0);
}
assert!(l[i * 5 + i] > 0.0);
}
assert_close(&mm(&l, &transpose(&l, 5, 5), 5, 5, 5), &CH_A5, 1e-13);
}
#[test]
fn cholesky_reads_only_the_lower_triangle_like_numpy() {
let mut a = CH_A3;
a[1] = 99.0;
a[2] = -77.0;
a[5] = 55.0;
assert_close(&cholesky(&a, 3).unwrap(), &CH_L3, 1e-13);
}
#[test]
fn cholesky_of_1x1_and_scalar_multiple_of_identity() {
assert_eq!(cholesky(&[9.0], 1).unwrap(), vec![3.0]);
assert_close(
&cholesky(&scaled(&eye(3), 4.0), 3).unwrap(),
&scaled(&eye(3), 2.0),
1e-15,
);
}
#[test]
fn cholesky_rejects_indefinite_and_semidefinite_input() {
// Indefinite: eigenvalues 3 and -1.
assert!(matches!(
cholesky(&[1.0, 2.0, 2.0, 1.0], 2),
Err(KernelError::NotPositiveDefinite)
));
// Positive semi-definite but singular (pivot is exactly 0): LAPACK also refuses.
assert!(matches!(
cholesky(&[1.0, 1.0, 1.0, 1.0], 2),
Err(KernelError::NotPositiveDefinite)
));
assert!(matches!(
cholesky(&[0.0; 4], 2),
Err(KernelError::NotPositiveDefinite)
));
assert!(matches!(
cholesky(&[-4.0], 1),
Err(KernelError::NotPositiveDefinite)
));
assert!(matches!(cholesky(&[1.0; 3], 2), Err(KernelError::Shape)));
}
#[test]
fn cholesky_hilbert_boundary_matches_numpy() {
// NumPy: float64 succeeds for n = 8, 12, 13 and fails for n = 14; float32 succeeds for n = 8, fails for n = 10.
for n in [8, 12, 13] {
assert!(cholesky(&hilbert(n), n).is_ok(), "f64 Hilbert {n}");
}
assert!(matches!(
cholesky(&hilbert(14), 14),
Err(KernelError::NotPositiveDefinite)
));
let h8: Vec<f32> = to_f32(&hilbert(8));
assert!(cholesky(&h8, 8).is_ok());
let h10: Vec<f32> = to_f32(&hilbert(10));
assert!(matches!(
cholesky(&h10, 10),
Err(KernelError::NotPositiveDefinite)
));
}
#[test]
fn cholesky_of_hilbert_8_reconstructs_to_machine_precision() {
let h = hilbert(8);
let l = cholesky(&h, 8).unwrap();
assert_close(&mm(&l, &transpose(&l, 8, 8), 8, 8, 8), &h, 1e-15);
}
// ---------------------------------------------------------------------------
// QR
// ---------------------------------------------------------------------------
#[test]
fn qr_tall_matches_numpy_including_signs() {
let (q, r, k) = qr(&QR_TALL, 5, 3).unwrap();
assert_eq!(k, 3);
assert_close(&q, &QR_TALL_Q, 1e-13);
assert_close(&r, &QR_TALL_R, 1e-13);
}
#[test]
fn qr_wide_matches_numpy_including_signs() {
let (q, r, k) = qr(&QR_WIDE, 3, 5).unwrap();
assert_eq!(k, 3);
assert_close(&q, &QR_WIDE_Q, 1e-13);
assert_close(&r, &QR_WIDE_R, 1e-13);
}
#[test]
fn qr_square_matches_numpy_including_signs() {
let (q, r, _) = qr(&QR_SQ, 4, 4).unwrap();
assert_close(&q, &QR_SQ_Q, 1e-13);
assert_close(&r, &QR_SQ_R, 1e-13);
}
#[test]
fn qr_reconstruction_orthonormality_and_triangularity() {
for (a, rows, cols) in [
(&QR_TALL[..], 5, 3),
(&QR_WIDE[..], 3, 5),
(&QR_SQ[..], 4, 4),
] {
let (q, r, k) = qr(a, rows, cols).unwrap();
assert_eq!(k, rows.min(cols));
assert_close(&mm(&q, &r, rows, k, cols), a, 1e-13);
assert_orthonormal_columns(&q, rows, k, 1e-13);
for i in 0..k {
for j in 0..i.min(cols) {
assert!(
r[i * cols + j].abs() < 1e-13,
"R[{i}][{j}] = {}",
r[i * cols + j]
);
}
}
}
}
#[test]
fn qr_of_upper_triangular_input_is_identity_times_input() {
// LAPACK leaves already-triangular columns alone (tau = 0), so Q = I and R = A with positive diagonal.
let (q, r, _) = qr(&QR_UP, 3, 3).unwrap();
assert_close(&q, &QR_UP_Q, 1e-15);
assert_close(&r, &QR_UP_R, 1e-15);
assert_eq!(q, eye(3));
}
#[test]
fn qr_with_a_zero_column_matches_numpy() {
let (q, r, _) = qr(&QR_ZC, 4, 3).unwrap();
assert_close(&mm(&q, &r, 4, 3, 3), &QR_ZC, 1e-13);
assert_orthonormal_columns(&q, 4, 3, 1e-13);
assert_close(&q, &QR_ZC_Q, 1e-13);
assert_close(&r, &QR_ZC_R, 1e-13);
}
#[test]
fn qr_of_tiny_scaled_matrix_is_scale_invariant() {
// Entries ~1e-20: an absolute epsilon cut-off would treat every column as zero.
let (q, r, _) = qr(&QR_TINY, 4, 3).unwrap();
assert_close(&q, &QR_TINY_Q, 1e-13);
assert_close(&scaled(&r, 1e20), &scaled(&QR_TINY_R, 1e20), 1e-13);
}
#[test]
fn qr_of_single_column_and_1x1() {
let (q, r, _) = qr(&[3.0, 4.0], 2, 1).unwrap();
assert_close(&q, &[-0.6, -0.8], 1e-15);
assert_close(&r, &[-5.0], 1e-15);
let (q, r, _) = qr(&[-2.0], 1, 1).unwrap();
assert_eq!((q[0], r[0]), (1.0, -2.0));
}
#[test]
fn qr_rejects_wrong_length() {
assert!(matches!(qr(&[1.0; 5], 2, 3), Err(KernelError::Shape)));
}
// ---------------------------------------------------------------------------
// lstsq
// ---------------------------------------------------------------------------
#[test]
fn lstsq_overdetermined_line_fit_matches_numpy() {
let x = lstsq(&LS_A, 6, 2, &LS_B, 1).unwrap();
assert_close(&x, &LS_X, 1e-12);
}
#[test]
fn lstsq_overdetermined_multiple_rhs_matches_numpy() {
let x = lstsq(&LS_A2, 5, 3, &LS_B2, 2).unwrap();
assert_close(&x, &LS_X2, 1e-12);
}
#[test]
fn lstsq_residual_is_orthogonal_to_column_space() {
let x = lstsq(&LS_A2, 5, 3, &LS_B2, 2).unwrap();
let residual = sub(&mm(&LS_A2, &x, 5, 3, 2), &LS_B2);
let normal = mm(&transpose(&LS_A2, 5, 3), &residual, 3, 5, 2);
assert!(max_abs(&normal) < 1e-12, "{normal:?}");
}
#[test]
fn lstsq_underdetermined_returns_numpys_minimum_norm_solution() {
let x = lstsq(&LS_UA, 2, 4, &LS_UB, 1).unwrap();
assert_close(&x, &LS_UX, 1e-12);
// It solves the system exactly...
assert_close(&mm(&LS_UA, &x, 2, 4, 1), &LS_UB, 1e-12);
// ...and lies in the row space (x = Aᵀ y), the defining property of the minimum-norm solution.
let y = Lu::factor(&mm(&LS_UA, &transpose(&LS_UA, 2, 4), 2, 4, 2), 2)
.unwrap()
.solve(&LS_UB, 1)
.unwrap();
assert_close(&x, &mm(&transpose(&LS_UA, 2, 4), &y, 4, 2, 1), 1e-12);
}
#[test]
fn lstsq_underdetermined_with_multiple_rhs() {
let x = lstsq(&LS_UA, 2, 4, &LS_UB2, 2).unwrap();
assert_close(&x, &LS_UX2, 1e-12);
}
#[test]
fn lstsq_rank_deficient_tall_returns_minimum_norm_solution() {
// Third column = first + second, so the system has a one-parameter family of least-squares solutions.
let x = lstsq(&LS_RA, 4, 3, &LS_RB, 1).unwrap();
assert_close(&x, &LS_RX, 1e-12);
}
#[test]
fn lstsq_rank_deficient_square_returns_minimum_norm_solution() {
let x = lstsq(&LS_SA, 3, 3, &LS_SB, 1).unwrap();
assert_close(&x, &LS_SX, 1e-12);
}
#[test]
fn lstsq_square_nonsingular_equals_solve() {
let x = lstsq(&LS_QA, 3, 3, &LS_QB, 1).unwrap();
assert_close(&x, &LS_QX, 1e-12);
}
#[test]
fn lstsq_rejects_mismatched_rhs() {
assert!(matches!(
lstsq(&LS_A, 6, 2, &[1.0; 5], 1),
Err(KernelError::Shape)
));
}
// ---------------------------------------------------------------------------
// eigh
// ---------------------------------------------------------------------------
#[test]
fn eigh_3x3_matches_numpy_values_and_vectors() {
let v = check_eigh(&EH_A3, 3, &EH_W3, 1e-12);
assert_columns_match_up_to_sign(&v, &EH_V3, 3, 3, 1e-10);
}
#[test]
fn eigh_5x5_matches_numpy_values_and_vectors() {
let v = check_eigh(&EH_A5, 5, &EH_W5, 1e-12);
assert_columns_match_up_to_sign(&v, &EH_V5, 5, 5, 1e-9);
}
#[test]
fn eigh_6x6_eigenvalues_match_numpy() {
check_eigh(&EH_A6, 6, &EH_W6, 1e-12);
}
#[test]
fn eigh_repeated_eigenvalues_still_give_an_orthonormal_basis() {
// [[3,1,1],[1,3,1],[1,1,3]] has eigenvalues 2, 2, 5.
check_eigh(&EH_REP1, 3, &EH_REP1_W, 1e-13);
// The all-ones 4x4 has eigenvalues 0, 0, 0, 4.
check_eigh(&EH_ONES4, 4, &EH_ONES4_W, 1e-13);
// 2 * identity: every vector is an eigenvector; the kernel must return an orthonormal basis.
let v = check_eigh(&scaled(&eye(3), 2.0), 3, &[2.0, 2.0, 2.0], 1e-15);
assert_close(&v, &eye(3), 1e-15);
}
#[test]
fn eigh_of_1x1_diagonal_and_already_sorted_inputs() {
let (w, v) = eigh(&[-3.5], 1).unwrap();
assert_eq!((w, v), (vec![-3.5], vec![1.0]));
let (w, _) = eigh(&[3.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 2.0], 3).unwrap();
assert_eq!(w, vec![-1.0, 2.0, 3.0]);
}
#[test]
fn eigh_is_scale_invariant_for_tiny_and_huge_matrices() {
// Absolute thresholds break this: a 1e-20 matrix looks "already diagonal" to an absolute test.
let tiny = scaled(&EH_A3, 1e-20);
let (w, v) = eigh(&tiny, 3).unwrap();
assert_close(&scaled(&w, 1e20), &scaled(&EH_W3_TINY, 1e20), 1e-12);
assert_orthonormal_columns(&v, 3, 3, 1e-12);
let av = mm(&scaled(&tiny, 1e20), &v, 3, 3, 3);
let mut lv = v.clone();
for r in 0..3 {
for c in 0..3 {
lv[r * 3 + c] *= w[c] * 1e20;
}
}
assert_close(&av, &lv, 1e-12);
let (w, _) = eigh(&scaled(&EH_A3, 1e20), 3).unwrap();
assert_close(&scaled(&w, 1e-20), &scaled(&EH_W3_BIG, 1e-20), 1e-12);
}
#[test]
fn eigh_f32_satisfies_the_eigen_equation() {
let a = to_f32(&EH_A5);
let (w, v) = eigh(&a, 5).unwrap();
assert_close_f32(&w, &EH_W5, 2e-5);
let (w64, v64): (Vec<f64>, Vec<f64>) = (
w.iter().map(|&x| x as f64).collect(),
v.iter().map(|&x| x as f64).collect(),
);
let av = mm(&EH_A5, &v64, 5, 5, 5);
for r in 0..5 {
for c in 0..5 {
assert!((av[r * 5 + c] - w64[c] * v64[r * 5 + c]).abs() < 1e-4);
}
}
}
#[test]
fn eigh_rejects_wrong_length() {
assert!(matches!(eigh(&[1.0; 5], 2), Err(KernelError::Shape)));
}
// ---------------------------------------------------------------------------
// SVD
// ---------------------------------------------------------------------------
#[test]
fn svd_tall_matches_numpy_values_and_vectors() {
check_svd(&SV_TALL, 4, 3, &SV_TALL_S, 1e-12);
let (u, _, vt, _) = svd(&SV_TALL, 4, 3).unwrap();
assert_columns_match_up_to_sign(&u, &SV_TALL_U, 4, 3, 1e-10);
assert_columns_match_up_to_sign(
&transpose(&vt, 3, 3),
&transpose(&SV_TALL_VT, 3, 3),
3,
3,
1e-10,
);
}
#[test]
fn svd_wide_matches_numpy_values() {
check_svd(&SV_WIDE, 3, 5, &SV_WIDE_S, 1e-12);
}
#[test]
fn svd_square_and_tall_7x4_match_numpy_values() {
check_svd(&SV_SQ, 3, 3, &SV_SQ_S, 1e-12);
check_svd(&SV_TALL2, 7, 4, &SV_TALL2_S, 1e-12);
}
#[test]
fn svd_rank_one_matrix_gives_orthonormal_factors() {
// [1,2,3] outer [4,5,6]: one non-zero singular value; U and Vt must still be orthonormal.
check_svd(&SV_R1, 3, 3, &SV_R1_S, 1e-12);
}
#[test]
fn svd_rank_two_tall_and_wide_matrices() {
check_svd(&SV_R2, 4, 3, &SV_R2_S, 1e-12);
check_svd(&SV_RW, 3, 6, &SV_RW_S, 1e-12);
}
#[test]
fn svd_zero_matrix_has_zero_singular_values_and_orthonormal_factors() {
let (u, s, vt, k) = svd(&[0.0; 6], 2, 3).unwrap();
assert_eq!((s.clone(), k), (vec![0.0, 0.0], 2));
assert_orthonormal_columns(&u, 2, 2, 1e-15);
assert_orthonormal_rows(&vt, 2, 3, 1e-15);
}
#[test]
fn svd_of_row_and_column_vectors() {
// 1x4: single singular value equal to the 2-norm (sqrt 30).
check_svd(&[1.0, 2.0, 3.0, 4.0], 1, 4, &SV_ROW_S, 1e-14);
check_svd(&[1.0, 2.0, 3.0, 4.0], 4, 1, &SV_ROW_S, 1e-14);
let (u, s, vt, _) = svd(&[-3.0, 4.0], 2, 1).unwrap();
assert_close(&s, &[5.0], 1e-15);
assert_close(&scaled(&u, 5.0 * vt[0]), &[-3.0, 4.0], 1e-15);
}
#[test]
fn svd_is_scale_invariant() {
check_svd(&scaled(&SV_SQ, 1e20), 3, 3, &SV_SQ_BIG_S, 1e-12 * 1e20);
let tiny = scaled(&SV_SQ, 1e-20);
let (_, s, _, _) = svd(&tiny, 3, 3).unwrap();
assert_close(&scaled(&s, 1e20), &scaled(&SV_SQ_TINY_S, 1e20), 1e-12);
}
#[test]
fn svd_of_hilbert_8_matches_numpy_to_absolute_machine_precision() {
let (_, s, _, _) = svd(&hilbert(8), 8, 8).unwrap();
assert_close(&s, &SV_H8_S, 1e-14);
// The smallest value (1.1e-10) is only defined to ~eps * s_max, so compare absolutely.
assert!((s[7] - SV_H8_S[7]).abs() < 1e-15);
}
#[test]
fn svd_f32_reconstructs() {
let a = to_f32(&SV_TALL);
let (u, s, vt, k) = svd(&a, 4, 3).unwrap();
assert_close_f32(&s, &SV_TALL_S, 2e-5);
let mut us = u.clone();
for r in 0..4 {
for c in 0..k {
us[r * k + c] *= s[c];
}
}
let rec = gemm(4, 3, k, &us, k as isize, 1, &vt, 3, 1);
assert_close_f32(&rec, &SV_TALL, 5e-5);
}
#[test]
fn svd_rejects_wrong_length() {
assert!(matches!(svd(&[1.0; 5], 2, 3), Err(KernelError::Shape)));
}
// ---------------------------------------------------------------------------
// pinv, matrix_rank
// ---------------------------------------------------------------------------
#[test]
fn pinv_tall_and_wide_match_numpy() {
assert_close(&pinv(&PI_TALL, 4, 3).unwrap(), &PI_TALL_P, 1e-12);
assert_close(&pinv(&PI_WIDE, 3, 5).unwrap(), &PI_WIDE_P, 1e-12);
}
#[test]
fn pinv_of_rank_deficient_matrices_matches_numpy() {
assert_close(&pinv(&PI_R1, 3, 3).unwrap(), &PI_R1_P, 1e-12);
assert_close(&pinv(&PI_R2, 4, 3).unwrap(), &PI_R2_P, 1e-12);
}
#[test]
fn pinv_satisfies_the_four_penrose_conditions() {
let (a, rows, cols) = (&PI_R2[..], 4, 3);
let p = pinv(a, rows, cols).unwrap();
let ap = mm(a, &p, rows, cols, rows);
let pa = mm(&p, a, cols, rows, cols);
assert_close(&mm(&ap, a, rows, rows, cols), a, 1e-12);
assert_close(&mm(&pa, &p, cols, cols, rows), &p, 1e-12);
assert_close(&ap, &transpose(&ap, rows, rows), 1e-12);
assert_close(&pa, &transpose(&pa, cols, cols), 1e-12);
}
#[test]
fn pinv_of_invertible_matrix_is_inverse_and_of_zero_is_zero() {
assert_close(&pinv(&DET_A4, 4, 4).unwrap(), &INV_A4, 1e-11);
assert_eq!(pinv(&[0.0; 6], 2, 3).unwrap(), vec![0.0; 6]);
}
#[test]
fn matrix_rank_of_assorted_matrices() {
assert_eq!(matrix_rank(&DET_A4, 4, 4).unwrap(), 4);
assert_eq!(matrix_rank(&SING3, 3, 3).unwrap(), 2);
assert_eq!(matrix_rank(&[1.0; 9], 3, 3).unwrap(), 1);
assert_eq!(matrix_rank(&[0.0; 9], 3, 3).unwrap(), 0);
assert_eq!(matrix_rank(&SV_R2, 4, 3).unwrap(), 2);
assert_eq!(matrix_rank(&SV_RW, 3, 6).unwrap(), 2);
assert_eq!(matrix_rank(&SV_WIDE, 3, 5).unwrap(), 3);
assert_eq!(matrix_rank(&SV_TALL, 4, 3).unwrap(), 3);
}
#[test]
fn matrix_rank_of_hilbert_matrices_matches_numpy_tolerance_rule() {
// NumPy: tol = s_max * max(M, N) * eps. H8 is numerically full rank; H12 (cond 1.8e16) has rank 11.
assert_eq!(matrix_rank(&hilbert(8), 8, 8).unwrap(), 8);
assert_eq!(matrix_rank(&hilbert(12), 12, 12).unwrap(), RANK_H12);
assert_eq!(RANK_H12, 11);
}
// ---------------------------------------------------------------------------
// norms
// ---------------------------------------------------------------------------
#[test]
fn vector_norms_match_numpy_for_every_order() {
let x = &NV_X[..];
assert_rel(vector_norm(x, NormOrd::One).unwrap(), NV_X_1, 1e-14);
assert_rel(vector_norm(x, NormOrd::Two).unwrap(), NV_X_2, 1e-14);
assert_rel(vector_norm(x, NormOrd::Fro).unwrap(), NV_X_2, 1e-14);
assert_eq!(vector_norm(x, NormOrd::Inf).unwrap(), NV_X_INF);
assert_eq!(vector_norm(x, NormOrd::NegInf).unwrap(), NV_X_NEGINF);
assert_eq!(vector_norm(x, NormOrd::Zero).unwrap(), NV_X_0);
assert_rel(vector_norm(x, NormOrd::NegOne).unwrap(), NV_X_NEG1, 1e-14);
assert_rel(vector_norm(x, NormOrd::NegTwo).unwrap(), NV_X_NEG2, 1e-14);
}
#[test]
fn vector_norms_with_zero_entries() {
let z = &NV_Z[..];
assert_rel(vector_norm(z, NormOrd::One).unwrap(), NV_Z_1, 1e-14);
assert_rel(vector_norm(z, NormOrd::Two).unwrap(), NV_Z_2, 1e-14);
assert_eq!(vector_norm(z, NormOrd::Inf).unwrap(), NV_Z_INF);
assert_eq!(vector_norm(z, NormOrd::NegInf).unwrap(), NV_Z_NEGINF);
assert_eq!(vector_norm(z, NormOrd::Zero).unwrap(), NV_Z_0);
assert_eq!(vector_norm(z, NormOrd::NegOne).unwrap(), NV_Z_NEG1);
assert_eq!(norm(z), vector_norm(z, NormOrd::Two).unwrap());
}
#[test]
fn vector_norm_rejects_nuclear_order_and_empty_negative_inf() {
assert!(matches!(
vector_norm(&[1.0], NormOrd::Nuc),
Err(KernelError::Shape)
));
assert!(matches!(
vector_norm::<f64>(&[], NormOrd::NegInf),
Err(KernelError::Shape)
));
assert_eq!(vector_norm::<f64>(&[], NormOrd::Two).unwrap(), 0.0);
}
#[test]
fn matrix_norms_match_numpy_for_every_order() {
let m = &NM_M[..];
let t = |ord| matrix_norm(m, 3, 4, ord).unwrap();
assert_rel(t(NormOrd::Fro), NM_M_FRO, 1e-14);
assert_rel(t(NormOrd::One), NM_M_1, 1e-14);
assert_rel(t(NormOrd::Inf), NM_M_INF, 1e-14);
assert_rel(t(NormOrd::Two), NM_M_2, 1e-12);
assert_rel(t(NormOrd::Nuc), NM_M_NUC, 1e-12);
assert_rel(t(NormOrd::NegOne), NM_M_NEG1, 1e-14);
assert_rel(t(NormOrd::NegInf), NM_M_NEGINF, 1e-14);
assert_rel(t(NormOrd::NegTwo), NM_M_NEG2, 1e-9);
}
#[test]
fn matrix_norm_rejects_vector_only_order_and_bad_length() {
assert!(matches!(
matrix_norm(&NM_M, 3, 4, NormOrd::Zero),
Err(KernelError::Shape)
));
assert!(matches!(
matrix_norm(&NM_M, 3, 3, NormOrd::One),
Err(KernelError::Shape)
));
}
#[test]
fn norm_axis_matches_numpy_for_each_order_and_axis() {
let m = &NA_M[..];
let cases: [(NormOrd, &[f64], &[f64]); 5] = [
(NormOrd::One, &NA_1_0, &NA_1_1),
(NormOrd::Two, &NA_2_0, &NA_2_1),
(NormOrd::Inf, &NA_INF_0, &NA_INF_1),
(NormOrd::NegInf, &NA_NEGINF_0, &NA_NEGINF_1),
(NormOrd::Zero, &NA_0_0, &NA_0_1),
];
for (ord, axis0, axis1) in cases {
assert_close(&norm_axis(m, 3, 3, 0, ord).unwrap(), axis0, 1e-14);
assert_close(&norm_axis(m, 3, 3, 1, ord).unwrap(), axis1, 1e-14);
}
assert!(matches!(
norm_axis(m, 3, 3, 2, NormOrd::One),
Err(KernelError::Shape)
));
}
// ---------------------------------------------------------------------------
// ill-conditioned behaviour
// ---------------------------------------------------------------------------
#[test]
fn hilbert_8_two_norm_condition_number_matches_numpy() {
// cond_2 = s_max / s_min with s_min ~ 1.1e-10: ~1e-5 relative agreement is all the data supports.
let (_, s, _, _) = svd(&hilbert(8), 8, 8).unwrap();
assert_rel(s[0] / s[7], COND_H8, 1e-6);
}
#[test]
fn hilbert_tolerances_scale_with_the_problem() {
// An ill-conditioned but invertible system: inversion succeeds, the residual is tiny, and the
// reported inverse is accurate only to about 1/cond relative precision.
let h = hilbert(8);
let x = inv(&h, 8).unwrap();
assert_close(&mm(&h, &x, 8, 8, 8), &eye(8), 1e-5);
// Error against the exact inverse is far above machine epsilon times n...
let rel = max_abs(&sub(&x, &INV_H8_EXACT)) / max_abs(&INV_H8_EXACT);
assert!(rel > 1e-14, "suspiciously exact: {rel:e}");
// ...but pinv through the SVD gives a comparably accurate answer.
let p = pinv(&h, 8, 8).unwrap();
let rel_pinv = max_abs(&sub(&p, &INV_H8_EXACT)) / max_abs(&INV_H8_EXACT);
assert!(rel_pinv < 1e-6, "{rel_pinv:e}");
}
// ===========================================================================
// New kernels: eig, eigvals, eigvalsh, slogdet, cond, multi_dot, tensordot, vdot, matrix_power
// ===========================================================================
/// Deterministic pseudo-random matrix; the same LCG is used in the NumPy generator script.
fn lcg_matrix(n: usize) -> Vec<f64> {
let mut state: u64 = 12345;
(0..n * n)
.map(|_| {
state = (state * 1_103_515_245 + 12_345) % (1 << 31);
(((state >> 16) % 2001) as f64 - 1000.0) / 1000.0
})
.collect()
}
/// Every expected eigenvalue is matched by a distinct computed one (eigenvalue order is not specified).
#[track_caller]
fn assert_eigs_match(re: &[f64], im: &[f64], exp_re: &[f64], exp_im: &[f64], tol: f64) {
assert_eq!(re.len(), exp_re.len());
let scale = max_abs(exp_re).max(max_abs(exp_im)).max(1.0);
let mut taken = vec![false; re.len()];
for (er, ei) in exp_re.iter().zip(exp_im) {
let mut best = None;
let mut best_distance = f64::INFINITY;
for j in 0..re.len() {
let distance = (re[j] - er).hypot(im[j] - ei);
if !taken[j] && distance < best_distance {
best = Some(j);
best_distance = distance;
}
}
let j = best.expect("more expected eigenvalues than computed");
taken[j] = true;
assert!(
best_distance <= tol * scale,
"eigenvalue {er}+{ei}i matched only to {}+{}i (distance {best_distance:e})\ncomputed re = {re:?}\ncomputed im = {im:?}",
re[j],
im[j]
);
}
}
/// `A v = lambda v` for every column, unit norm, largest component real and positive, and
/// conjugate pairs adjacent with the positive imaginary part first.
#[track_caller]
fn check_eig(a: &[f64], n: usize, tol: f64) -> (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>) {
let (re, im, vr, vi) = eig(a, n).unwrap();
let scale = max_abs(a).max(1.0);
let arr = mm(a, &vr, n, n, n);
let ari = mm(a, &vi, n, n, n);
for j in 0..n {
let mut norm_sq = 0.0;
let mut largest = (0.0, 0.0, 0.0);
for r in 0..n {
let (x, y) = (vr[r * n + j], vi[r * n + j]);
norm_sq += x * x + y * y;
if x * x + y * y > largest.0 {
largest = (x * x + y * y, x, y);
}
let want_re = re[j] * x - im[j] * y;
let want_im = re[j] * y + im[j] * x;
assert!(
(arr[r * n + j] - want_re).abs() <= tol * scale
&& (ari[r * n + j] - want_im).abs() <= tol * scale,
"A v != lambda v for eigenvalue {j} ({}+{}i), row {r}",
re[j],
im[j]
);
}
assert!(
(norm_sq - 1.0).abs() < 1e-12,
"column {j} has squared norm {norm_sq}"
);
assert!(
largest.1 > 0.0 && largest.2 == 0.0,
"column {j}: largest component {largest:?}"
);
if im[j] > 0.0 {
assert!(
j + 1 < n && re[j + 1] == re[j] && im[j + 1] == -im[j],
"pair not adjacent at {j}"
);
for r in 0..n {
assert_eq!(
vr[r * n + j + 1],
vr[r * n + j],
"partner vector is the conjugate"
);
assert_eq!(
vi[r * n + j + 1],
-vi[r * n + j],
"partner vector is the conjugate"
);
}
}
}
(re, im, vr, vi)
}
#[test]
fn eig_real_spectrum_of_nonsymmetric_3x3() {
let (re, im, _, _) = check_eig(&E3, 3, 1e-11);
assert_eigs_match(&re, &im, &E3_RE, &E3_IM, 1e-10);
assert_eigs_match(&re, &im, &[1.0, 2.0, 3.0], &[0.0; 3], 1e-10);
}
#[test]
fn eig_4x4_with_two_complex_pairs_matches_numpy() {
assert!(E4_IM.iter().filter(|x| **x != 0.0).count() == 4);
let (re, im, _, _) = check_eig(&E4, 4, 1e-12);
assert_eigs_match(&re, &im, &E4_RE, &E4_IM, 1e-11);
}
#[test]
fn eig_5x5_with_one_complex_pair_matches_numpy() {
assert!(E5_IM.iter().filter(|x| **x != 0.0).count() == 2);
let (re, im, _, _) = check_eig(&E5, 5, 1e-12);
assert_eigs_match(&re, &im, &E5_RE, &E5_IM, 1e-11);
}
#[test]
fn eig_6x6_matches_numpy() {
let (re, im, _, _) = check_eig(&E6, 6, 1e-12);
assert_eigs_match(&re, &im, &E6_RE, &E6_IM, 1e-11);
}
#[test]
fn eig_complex_pairs_are_adjacent_with_positive_imaginary_part_first() {
let (_, im, _, _) = eig(&E6, 6).unwrap();
let mut j = 0;
while j < 6 {
if im[j] != 0.0 {
assert!(im[j] > 0.0 && im[j + 1] == -im[j], "{im:?}");
j += 2;
} else {
j += 1;
}
}
}
#[test]
fn eig_rotation_matrix_has_eigenvalues_plus_minus_i() {
let (re, im, vr, vi) = check_eig(&E_ROT, 2, 1e-15);
assert_eigs_match(&re, &im, &E_ROT_RE, &E_ROT_IM, 1e-15);
// For lambda = +i the eigenvector is (1, -i)/sqrt(2), as NumPy reports.
assert!(im[0] > 0.0);
let h = core::f64::consts::FRAC_1_SQRT_2;
assert_close(&[vr[0], vr[2], vi[0], vi[2]], &[h, 0.0, 0.0, -h], 1e-15);
}
#[test]
fn eig_cyclic_permutation_has_the_fourth_roots_of_unity() {
let (re, im, _, _) = check_eig(&E_PERM4, 4, 1e-14);
assert_eigs_match(
&re,
&im,
&[1.0, -1.0, 0.0, 0.0],
&[0.0, 0.0, 1.0, -1.0],
1e-14,
);
assert_eigs_match(&re, &im, &E_PERM4_RE, &E_PERM4_IM, 1e-14);
}
#[test]
fn eig_upper_triangular_returns_the_diagonal() {
let (re, im, _, _) = check_eig(&E_UPPER4, 4, 1e-13);
assert_eigs_match(&re, &im, &E_UPPER4_RE, &E_UPPER4_IM, 1e-14);
assert_eigs_match(&re, &im, &[4.0, -2.0, 7.0, 0.5], &[0.0; 4], 1e-14);
}
#[test]
fn eig_companion_matrix_recovers_roots_of_a_polynomial() {
// (x-1)(x-2)...(x-6): eigenvalues of the companion matrix are the roots (mildly ill-conditioned).
let (re, im, _, _) = eig(&E_COMP6, 6).unwrap();
assert_eigs_match(&re, &im, &[1.0, 2.0, 3.0, 4.0, 5.0, 6.0], &[0.0; 6], 1e-10);
}
#[test]
fn eig_defective_jordan_block_still_satisfies_the_eigen_equation() {
let a = [2.0, 1.0, 0.0, 2.0];
let (re, im, vr, _) = check_eig(&a, 2, 1e-14);
assert_eq!(re, vec![2.0, 2.0]);
assert_eq!(im, vec![0.0, 0.0]);
// The only eigenvector direction is e1; both columns must be parallel to it.
assert_close(&[vr[0], vr[2], vr[1], vr[3]], &[1.0, 0.0, 1.0, 0.0], 1e-14);
}
#[test]
fn eig_identity_zero_and_1x1() {
let (re, im, vr, vi) = eig(&eye(4), 4).unwrap();
assert_eq!((re, im), (vec![1.0; 4], vec![0.0; 4]));
assert_eq!(vr, eye(4));
assert!(vi.iter().all(|x| *x == 0.0));
let (re, _, vr, _) = eig(&[0.0; 9], 3).unwrap();
assert_eq!(re, vec![0.0; 3]);
assert_eq!(vr, eye(3));
let (re, im, vr, vi) = eig(&[-3.0], 1).unwrap();
assert_eq!(
(re, im, vr, vi),
(vec![-3.0], vec![0.0], vec![1.0], vec![0.0])
);
let (re, _, _, _) = eig::<f64>(&[], 0).unwrap();
assert!(re.is_empty());
}
#[test]
fn eig_of_symmetric_matrix_agrees_with_eigh() {
let (re, im, vr, vi) = check_eig(&E_SYM5, 5, 1e-12);
assert!(im.iter().all(|x| *x == 0.0) && vi.iter().all(|x| *x == 0.0));
let mut sorted = re.clone();
sorted.sort_by(f64::total_cmp);
assert_close(&sorted, &E_SYM5_W, 1e-12);
// Eigenvectors of a symmetric matrix with distinct eigenvalues are orthonormal.
assert_orthonormal_columns(&vr, 5, 5, 1e-10);
}
#[test]
fn eig_30x30_pseudo_random_matrix_matches_numpy() {
let a = lcg_matrix(30);
let (re, im, _, _) = check_eig(&a, 30, 1e-11);
assert_eigs_match(&re, &im, &E30_RE, &E30_IM, 1e-11);
// Trace = sum of eigenvalues; determinant = product (complex arithmetic).
assert!((re.iter().sum::<f64>() - E30_TRACE).abs() < 1e-11);
let (mut pr, mut pi) = (1.0_f64, 0.0_f64);
for (x, y) in re.iter().zip(&im) {
(pr, pi) = (pr * x - pi * y, pr * y + pi * x);
}
assert_rel(pr, E30_DET, 1e-9);
assert!(pi.abs() < 1e-9 * E30_DET.abs());
}
#[test]
fn eig_is_scale_invariant() {
let (re, im, _, _) = eig(&scaled(&E4, 1e-20), 4).unwrap();
assert_eigs_match(
&scaled(&re, 1e20),
&scaled(&im, 1e20),
&E4_RE,
&E4_IM,
1e-11,
);
assert_eigs_match(
&scaled(&re, 1e20),
&scaled(&im, 1e20),
&scaled(&E_C4_SCALED_RE, 1e20),
&scaled(&E_C4_SCALED_IM, 1e20),
1e-11,
);
let (re, im, _, _) = eig(&scaled(&E4, 1e20), 4).unwrap();
assert_eigs_match(
&scaled(&re, 1e-20),
&scaled(&im, 1e-20),
&E4_RE,
&E4_IM,
1e-11,
);
}
#[test]
fn eig_f32_complex_spectrum() {
let a = to_f32(&E4);
let (re, im, _, _) = eig(&a, 4).unwrap();
let (re, im): (Vec<f64>, Vec<f64>) = (
re.iter().map(|&x| f64::from(x)).collect(),
im.iter().map(|&x| f64::from(x)).collect(),
);
assert_eigs_match(&re, &im, &E4_RE, &E4_IM, 1e-4);
}
#[test]
fn eig_rejects_nan_inf_and_bad_length() {
assert!(matches!(
eig(&[1.0, f64::NAN, 0.0, 1.0], 2),
Err(KernelError::NoConvergence)
));
assert!(matches!(
eig(&[1.0, f64::INFINITY, 0.0, 1.0], 2),
Err(KernelError::NoConvergence)
));
assert!(matches!(eig(&[1.0; 5], 2), Err(KernelError::Shape)));
}
#[test]
fn eigvals_matches_numpy_and_eig() {
let (re, im) = eigvals(&E6, 6).unwrap();
assert_eigs_match(&re, &im, &E6_RE, &E6_IM, 1e-11);
let (re2, im2, _, _) = eig(&E6, 6).unwrap();
assert_eq!((re, im), (re2, im2));
}
#[test]
fn eigvalsh_matches_numpy_and_is_ascending() {
let w = eigvalsh(&EH_A6, 6).unwrap();
assert_close(&w, &EH_W6, 1e-12);
assert!(w.windows(2).all(|p| p[0] <= p[1]));
assert_close(&eigvalsh(&EH_ONES4, 4).unwrap(), &EH_ONES4_W, 1e-13);
}
#[test]
fn slogdet_matches_numpy() {
for (a, n, sign, log) in [
(&SLD_A3[..], 3, SLD_A3_SIGN, SLD_A3_LOG),
(&SLD_A4[..], 4, SLD_A4_SIGN, SLD_A4_LOG),
(&SLD_NEG[..], 2, SLD_NEG_SIGN, SLD_NEG_LOG),
] {
let (s, l) = slogdet(a, n).unwrap();
assert_eq!(s, sign);
assert!((l - log).abs() < 1e-13 * log.abs().max(1.0), "{l} vs {log}");
}
assert_eq!(SLD_NEG_SIGN, -1.0);
}
#[test]
fn slogdet_of_singular_matrix_is_zero_and_minus_infinity() {
assert_eq!(
slogdet(&[1.0, 2.0, 2.0, 4.0], 2).unwrap(),
(0.0, f64::NEG_INFINITY)
);
assert_eq!(slogdet(&[0.0; 9], 3).unwrap(), (0.0, f64::NEG_INFINITY));
}
#[test]
fn slogdet_survives_determinants_that_overflow() {
// diag(1e200, 1e200, 1e-5): det = 1e395 overflows f64, the log does not.
let a = [1e200, 0.0, 0.0, 0.0, 1e200, 0.0, 0.0, 0.0, 1e-5];
assert_eq!(det_of(&a, 3), f64::INFINITY);
let (s, l) = slogdet(&a, 3).unwrap();
assert_eq!(s, SLD_BIG_SIGN);
assert!((l - SLD_BIG_LOG).abs() < 1e-10, "{l} vs {SLD_BIG_LOG}");
}
#[test]
fn slogdet_edge_cases_and_consistency_with_det() {
assert_eq!(slogdet(&[-2.0], 1).unwrap(), (-1.0, 2.0_f64.ln()));
assert_eq!(slogdet::<f64>(&[], 0).unwrap(), (1.0, 0.0));
let (s, l) = slogdet(&DET_A4, 4).unwrap();
assert_rel(s * l.exp(), DET_A4_V, 1e-12);
// Permutation sign flows through.
assert_eq!(slogdet(&DET_P4, 4).unwrap(), (-1.0, 0.0));
assert!(matches!(slogdet(&[1.0; 3], 2), Err(KernelError::Shape)));
}
#[test]
fn cond_two_norm_matches_numpy_square_tall_and_wide() {
assert_rel(cond(&DET_A4, 4, 4).unwrap(), COND_A4_2, 1e-11);
assert_rel(cond(&DET_A3, 3, 3).unwrap(), COND_A3_2, 1e-11);
assert_rel(cond(&SV_TALL, 4, 3).unwrap(), COND_TALL, 1e-11);
assert_rel(cond(&SV_WIDE, 3, 5).unwrap(), COND_WIDE, 1e-11);
}
#[test]
fn cond_other_norms_match_numpy() {
for (a, n, one, inf, fro, neg2) in [
(
&DET_A4[..],
4,
COND_A4_1,
COND_A4_INF,
COND_A4_FRO,
COND_A4_NEG2,
),
(
&DET_A3[..],
3,
COND_A3_1,
COND_A3_INF,
COND_A3_FRO,
COND_A3_NEG2,
),
] {
assert_rel(cond_ord(a, n, n, NormOrd::One).unwrap(), one, 1e-11);
assert_rel(cond_ord(a, n, n, NormOrd::Inf).unwrap(), inf, 1e-11);
assert_rel(cond_ord(a, n, n, NormOrd::Fro).unwrap(), fro, 1e-11);
assert_rel(cond_ord(a, n, n, NormOrd::NegTwo).unwrap(), neg2, 1e-11);
let two = cond_ord(a, n, n, NormOrd::Two).unwrap();
assert_rel(two * neg2, 1.0, 1e-12);
}
}
#[test]
fn cond_of_singular_matrix_is_huge_or_infinite() {
// NumPy: 2.5e16 for the 2-norm (rounding noise in s_min) and inf for ord=1 (LinAlgError caught).
assert!(cond(&[1.0, 2.0, 2.0, 4.0], 2, 2).unwrap() > 1e15);
assert_eq!(COND_SING_1, f64::INFINITY);
assert_eq!(
cond_ord(&[1.0, 2.0, 2.0, 4.0], 2, 2, NormOrd::One).unwrap(),
f64::INFINITY
);
assert_eq!(cond(&[0.0; 4], 2, 2).unwrap(), f64::INFINITY);
}
#[test]
fn cond_of_orthogonal_and_identity_matrices_is_one() {
assert_rel(cond(&eye(5), 5, 5).unwrap(), 1.0, 1e-14);
let (q, _, _) = qr(&DET_A4, 4, 4).unwrap();
assert_rel(cond(&q, 4, 4).unwrap(), 1.0, 1e-12);
}
#[test]
fn cond_hilbert_8_matches_numpy_within_conditioning() {
assert_rel(cond(&hilbert(8), 8, 8).unwrap(), COND_H8, 1e-6);
// Another way: ||H||_1 ||H^-1||_1 with the exact inverse, to ~7 digits.
let exact_one = {
let h = matrix_norm(&hilbert(8), 8, 8, NormOrd::One).unwrap();
h * matrix_norm(&INV_H8_EXACT, 8, 8, NormOrd::One).unwrap()
};
assert_rel(
cond_ord(&hilbert(8), 8, 8, NormOrd::One).unwrap(),
exact_one,
1e-4,
);
}
#[test]
fn cond_rejects_unsupported_requests() {
assert!(matches!(
cond_ord(&[1.0; 6], 2, 3, NormOrd::One),
Err(KernelError::Shape)
));
assert!(matches!(
cond_ord(&eye(2), 2, 2, NormOrd::Zero),
Err(KernelError::Shape)
));
assert!(matches!(cond::<f64>(&[], 0, 0), Err(KernelError::Shape)));
}
fn chain<'a>(items: &'a [(&'a [f64], usize, usize)]) -> Vec<(&'a [f64], usize, usize)> {
items.to_vec()
}
#[test]
fn multi_dot_three_matrices_matches_numpy() {
let items = [(&MD3_0[..], 4, 12), (&MD3_1[..], 12, 2), (&MD3_2[..], 2, 6)];
let (out, rows, cols) = multi_dot(&chain(&items)).unwrap();
assert_eq!((rows, cols), (4, 6));
assert_close(&out, &MD3_OUT, 1e-13);
}
#[test]
fn multi_dot_five_matrices_matches_numpy() {
let items = [
(&MD5_0[..], 3, 8),
(&MD5_1[..], 8, 2),
(&MD5_2[..], 2, 9),
(&MD5_3[..], 9, 4),
(&MD5_4[..], 4, 5),
];
let (out, rows, cols) = multi_dot(&chain(&items)).unwrap();
assert_eq!((rows, cols), (3, 5));
assert_close(&out, &MD5_OUT, 1e-12);
}
#[test]
fn multi_dot_integer_chain_is_exact() {
let items: [(&[i64], usize, usize); 4] = [
(&MDI_0, 2, 5),
(&MDI_1, 5, 3),
(&MDI_2, 3, 4),
(&MDI_3, 4, 2),
];
let (out, rows, cols) = multi_dot(&items).unwrap();
assert_eq!((rows, cols), (2, 2));
assert_eq!(out, MDI_OUT.to_vec());
}
#[test]
fn multi_dot_two_matrices_is_a_plain_product() {
let items = [(&G34_A[..], 3, 4), (&G42_B[..], 4, 2)];
let (out, _, _) = multi_dot(&chain(&items)).unwrap();
assert_close(&out, &G32_C, 1e-14);
}
#[test]
fn matrix_chain_order_finds_the_documented_optimal_costs() {
// Classic textbook chain (CLRS 15.2): 6 matrices, optimum 15125 as ((A1 (A2 A3)) ((A4 A5) A6)).
let (cost, splits) = matrix_chain_order(&[30, 35, 15, 5, 10, 20, 25]);
assert_eq!(cost, 15125);
assert_eq!(splits[5], 2); // A1..A3 | A4..A6
assert_eq!(splits[8], 1); // A2 | A3
assert_eq!(splits[23], 4); // (A4 A5) | A6
// Four matrices 40x20, 20x30, 30x10, 10x30: optimum 26000.
assert_eq!(matrix_chain_order(&[40, 20, 30, 10, 30]).0, 26000);
// Our 4x12, 12x2, 2x6 example: (AB)C costs 144, A(BC) costs 432.
let (cost, splits) = matrix_chain_order(&[4, 12, 2, 6]);
assert_eq!(cost, 144);
assert_eq!(splits[2], 1);
// A single matrix costs nothing.
assert_eq!(matrix_chain_order(&[3, 4]).0, 0);
// Running left to right on 10x30 @ 30x5 @ 5x60 costs 4500; the other order 27000.
assert_eq!(matrix_chain_order(&[10, 30, 5, 60]).0, 4500);
}
#[test]
fn multi_dot_rejects_bad_chains() {
let a = (&G34_A[..], 3, 4);
assert!(matches!(multi_dot(&[a]), Err(KernelError::Shape)));
assert!(matches!(multi_dot::<f64>(&[]), Err(KernelError::Shape)));
assert!(matches!(multi_dot(&[a, a]), Err(KernelError::Shape))); // 4 != 3
assert!(matches!(
multi_dot(&[a, (&G42_B[..], 3, 3)]),
Err(KernelError::Shape)
)); // wrong length
}
#[test]
fn permute_axes_matches_numpy_transpose() {
let out = permute_axes(&PA_IN, &[2, 3, 4, 2], &[2, 0, 3, 1]).unwrap();
assert_eq!(out, PA_OUT.to_vec());
let out = permute_axes(&TD1_A, &[2, 3, 4], &[2, 0, 1]).unwrap();
assert_eq!(out, PA_OUT2.to_vec());
// Identity permutation and a 2-d transpose.
assert_eq!(
permute_axes(&[1, 2, 3, 4, 5, 6], &[2, 3], &[0, 1]).unwrap(),
vec![1, 2, 3, 4, 5, 6]
);
assert_eq!(
permute_axes(&[1, 2, 3, 4, 5, 6], &[2, 3], &[1, 0]).unwrap(),
vec![1, 4, 2, 5, 3, 6]
);
// 0-d and empty buffers.
assert_eq!(permute_axes(&[7], &[], &[]).unwrap(), vec![7]);
assert!(
permute_axes::<f64>(&[], &[2, 0], &[1, 0])
.unwrap()
.is_empty()
);
}
#[test]
fn permute_axes_rejects_bad_permutations() {
assert!(matches!(
permute_axes(&[1, 2], &[2], &[1]),
Err(KernelError::Shape)
));
assert!(matches!(
permute_axes(&[1, 2, 3, 4], &[2, 2], &[0, 0]),
Err(KernelError::Shape)
));
assert!(matches!(
permute_axes(&[1, 2, 3, 4], &[2, 2], &[0]),
Err(KernelError::Shape)
));
assert!(matches!(
permute_axes(&[1, 2, 3], &[2, 2], &[1, 0]),
Err(KernelError::Shape)
));
}
#[test]
fn tensordot_two_axes_of_3d_arrays_matches_numpy() {
let (out, shape) = tensordot(&TD1_A, &[2, 3, 4], &TD1_B, &[4, 3, 2], &[1, 2], &[1, 0]).unwrap();
assert_eq!(shape, vec![2, 2]);
assert_close(&out, &TD1_OUT, 1e-13);
}
#[test]
fn tensordot_single_axis_is_a_matrix_product() {
let (out, shape) = tensordot(&TD2_A, &[3, 4], &TD2_B, &[4, 2], &[1], &[0]).unwrap();
assert_eq!(shape, vec![3, 2]);
assert_close(&out, &TD2_OUT, 1e-14);
}
#[test]
fn tensordot_without_contraction_is_the_outer_product() {
let (out, shape) = tensordot(&TD3_A, &[2, 3], &TD3_B, &[4], &[], &[]).unwrap();
assert_eq!(shape, vec![2, 3, 4]);
assert_close(&out, &TD3_OUT, 1e-15);
}
#[test]
fn tensordot_full_contraction_gives_a_scalar() {
let (out, shape) = tensordot(&TD3_A, &[2, 3], &TD4_B, &[2, 3], &[0, 1], &[0, 1]).unwrap();
assert!(shape.is_empty());
assert_eq!(out.len(), 1);
assert_rel(out[0], TD4_OUT, 1e-13);
}
#[test]
fn tensordot_contracting_first_axes_transposes_operands() {
let (out, shape) = tensordot(&TD5_A, &[3, 4], &TD5_B, &[3, 5], &[0], &[0]).unwrap();
assert_eq!(shape, vec![4, 5]);
assert_close(&out, &TD5_OUT, 1e-14);
}
#[test]
fn tensordot_non_adjacent_axes_of_two_3d_arrays() {
let (out, shape) = tensordot(&TD6_A, &[2, 3, 4], &TD6_B, &[4, 2, 5], &[0, 2], &[1, 0]).unwrap();
assert_eq!(shape, vec![3, 5]);
assert_close(&out, &TD6_OUT, 1e-13);
}
#[test]
fn tensordot_trailing_by_leading_axis_keeps_all_free_axes() {
// axes=1 in NumPy: last axis of a with first axis of b.
let (out, shape) = tensordot(&TD8_A, &[2, 3, 4], &TD8_B, &[4, 2, 2], &[2], &[0]).unwrap();
assert_eq!(shape, vec![2, 3, 2, 2]);
assert_close(&out, &TD8_OUT, 1e-13);
}
#[test]
fn tensordot_integer_arrays_are_exact() {
let (out, shape) = tensordot(&TD7_A, &[2, 3, 2], &TD7_B, &[2, 2, 3], &[0, 2], &[1, 0]).unwrap();
assert_eq!(shape, vec![3, 3]);
assert_eq!(out, TD7_OUT.to_vec());
}
#[test]
fn tensordot_with_zero_length_contraction_gives_zeros() {
let empty: [f64; 0] = [];
let (out, shape) = tensordot(&empty, &[2, 0], &empty, &[0, 3], &[1], &[0]).unwrap();
assert_eq!(shape, vec![2, 3]);
assert_eq!(out, vec![0.0; 6]);
}
#[test]
fn tensordot_rejects_inconsistent_axes() {
let a = [1.0; 6];
let b = [1.0; 6];
// Dimension mismatch (3 vs 2), repeated axis, out-of-range axis, unequal axis list lengths, bad buffer length.
assert!(matches!(
tensordot(&a, &[2, 3], &b, &[2, 3], &[1], &[0]),
Err(KernelError::Shape)
));
assert!(matches!(
tensordot(&a, &[2, 3], &b, &[3, 2], &[1, 1], &[0, 0]),
Err(KernelError::Shape)
));
assert!(matches!(
tensordot(&a, &[2, 3], &b, &[3, 2], &[2], &[0]),
Err(KernelError::Shape)
));
assert!(matches!(
tensordot(&a, &[2, 3], &b, &[3, 2], &[1], &[]),
Err(KernelError::Shape)
));
assert!(matches!(
tensordot(&a, &[2, 4], &b, &[4, 2], &[1], &[0]),
Err(KernelError::Shape)
));
}
#[test]
fn vdot_flattens_both_operands() {
// (2,3) . (3,2) flattened: shapes differ but sizes agree, exactly like np.vdot.
assert_rel(vdot(&VD_A, &VD_B).unwrap(), VD_OUT, 1e-14);
assert_eq!(vdot(&[1, 2, 3], &[4, 5, 6]).unwrap(), 32);
assert_eq!(vdot::<f64>(&[], &[]).unwrap(), 0.0);
assert_eq!(vdot(&[1.0_f32, 2.0], &[0.5, 0.25]).unwrap(), 1.0);
}
#[test]
fn vdot_rejects_different_lengths() {
assert!(matches!(vdot(&[1.0, 2.0], &[1.0]), Err(KernelError::Shape)));
}
#[test]
fn matrix_power_positive_exponents_match_numpy() {
assert_close(&matrix_power(&MP_A, 3, 1).unwrap(), &MP_A, 1e-15);
assert_close(&matrix_power(&MP_A, 3, 2).unwrap(), &MP_P2, 1e-13);
assert_close(&matrix_power(&MP_A, 3, 5).unwrap(), &MP_P5, 1e-12);
assert_close(&matrix_power(&MP_A, 3, 10).unwrap(), &MP_P10, 1e-11);
}
#[test]
fn matrix_power_negative_exponents_use_the_inverse() {
assert_close(&matrix_power(&MP_A, 3, -1).unwrap(), &MP_PM1, 1e-12);
assert_close(&matrix_power(&MP_A, 3, -3).unwrap(), &MP_PM3, 1e-11);
assert_close(&matrix_power(&MP_A, 3, -7).unwrap(), &MP_PM7, 1e-10);
// A^-3 A^3 = I.
let product = mm(
&matrix_power(&MP_A, 3, -3).unwrap(),
&matrix_power(&MP_A, 3, 3).unwrap(),
3,
3,
3,
);
assert_close(&product, &eye(3), 1e-11);
}
#[test]
fn matrix_power_zero_is_identity_even_for_singular_matrices() {
assert_eq!(matrix_power(&MP_A, 3, 0).unwrap(), eye(3));
assert_eq!(matrix_power(&[1.0, 2.0, 2.0, 4.0], 2, 0).unwrap(), eye(2));
// Positive powers of a singular matrix are fine: [[1,2],[2,4]]^2 = 5 * itself.
assert_close(
&matrix_power(&[1.0, 2.0, 2.0, 4.0], 2, 2).unwrap(),
&[5.0, 10.0, 10.0, 20.0],
1e-15,
);
}
#[test]
fn matrix_power_negative_of_singular_matrix_is_an_error() {
assert!(matches!(
matrix_power(&[1.0, 2.0, 2.0, 4.0], 2, -1),
Err(KernelError::Singular)
));
assert!(matches!(
matrix_power(&[1.0; 5], 2, 2),
Err(KernelError::Shape)
));
}
#[test]
fn matrix_power_of_fibonacci_matrix_is_exact() {
let fib = [1.0, 1.0, 1.0, 0.0];
assert_eq!(matrix_power(&fib, 2, 30).unwrap(), MP_FIB30.to_vec());
assert_eq!(MP_FIB30[0], 1_346_269.0);
assert_close(&matrix_power(&fib, 2, -5).unwrap(), &MP_FIBM5, 1e-13);
}
#[test]
fn eigh_and_svd_report_non_finite_input_instead_of_returning_garbage() {
let bad = [1.0, f64::NAN, f64::NAN, 1.0];
assert!(matches!(eigh(&bad, 2), Err(KernelError::NoConvergence)));
assert!(matches!(svd(&bad, 2, 2), Err(KernelError::NoConvergence)));
assert!(matches!(
eigvalsh(&[f64::INFINITY, 0.0, 0.0, 1.0], 2),
Err(KernelError::NoConvergence)
));
assert!(matches!(
matrix_rank(&bad, 2, 2),
Err(KernelError::NoConvergence)
));
}