algebra_linear 1.2.0

//! Native BLAS: the level-1, -2 and -3 operations the LAPACK routines need.
//!
//! # Column-major indexing
//!
//! Every matrix here is column-major with `nrow` rows, so element `(i, j)` is at
//! `data[i + j * nrow]`. Getting that detail wrong is the classic source of a
//! plausible-looking but wrong kernel, so each routine states its indexing in a
//! comment and the tests exercise non-square shapes where a transposed index
//! cannot coincidentally be right.

/// `y := alpha * x + y` (daxpy), where `y` and `x` are the same length.
///
/// # Panics
///
/// Panics if the slices differ in length.
pub fn daxpy(alpha: f64, x: &[f64], y: &mut [f64]) {
    assert_eq!(x.len(), y.len(), "daxpy: vectors must have the same length");
    for i in 0..y.len() {
        y[i] += alpha * x[i];
    }
}

/// `x := alpha * x` (dscal).
pub fn dscal(alpha: f64, x: &mut [f64]) {
    for v in x.iter_mut() {
        *v *= alpha;
    }
}

/// `y := x` (dcopy).
///
/// # Panics
///
/// Panics if the slices differ in length.
pub fn dcopy(x: &[f64], y: &mut [f64]) {
    assert_eq!(x.len(), y.len(), "dcopy: vectors must have the same length");
    y.copy_from_slice(x);
}

/// `x . y` (ddot).
///
/// BLAS's `ddot` returns 0 for mismatched lengths by reading the shorter one;
/// this does the same rather than panicking, since that is the documented
/// behavior callers may rely on.
pub fn ddot(x: &[f64], y: &[f64]) -> f64 {
    x.iter().zip(y.iter()).map(|(a, b)| a * b).sum()
}

/// Euclidean norm `sqrt(x . x)` (dnrm2).
///
/// Scaled to avoid overflow and underflow for very large or very small vectors,
/// as the reference implementation is.
pub fn dnrm2(x: &[f64]) -> f64 {
    // Blue's scaling: find the largest magnitude first, then accumulate in units
    // of it, so a vector of 1e200s does not overflow in the squares.
    let mut scale = 0.0f64;
    for &v in x {
        let a = v.abs();
        if a > scale {
            scale = a;
        }
    }
    if scale == 0.0 || !scale.is_finite() {
        // All zeros, or a non-finite entry (which propagates as it should).
        return x.iter().map(|v| v * v).sum::<f64>().sqrt();
    }
    let inv = 1.0 / scale;
    let ss: f64 = x.iter().map(|v| (v * inv) * (v * inv)).sum();
    scale * ss.sqrt()
}

/// Index of the first entry with the largest absolute value (idamax).
///
/// 0-based, unlike BLAS's `idamax`, which is 1-based. The *first* maximum wins,
/// which matters when equal magnitudes appear and is what BLAS specifies. An
/// empty slice yields 0.
pub fn idamax(x: &[f64]) -> usize {
    let mut best = 0usize;
    let mut best_abs = f64::NEG_INFINITY;
    for (i, &v) in x.iter().enumerate() {
        let a = v.abs();
        if a > best_abs {
            best_abs = a;
            best = i;
        }
    }
    best
}

/// Sum of absolute values `|x|_1` (dasum).
pub fn dasum(x: &[f64]) -> f64 {
    x.iter().map(|v| v.abs()).sum()
}

/// `y := alpha * A * x + beta * y` (dgemv, trans = 'N').
///
/// `a` is `m x n` column-major with `m` rows.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn dgemv(alpha: f64, a: &[f64], x: &[f64], beta: f64, y: &mut [f64], m: usize, n: usize) {
    assert_eq!(a.len(), m * n, "dgemv: a must be m*n");
    assert_eq!(x.len(), n, "dgemv: x must have n entries");
    assert_eq!(y.len(), m, "dgemv: y must have m entries");

    // Scale y first so it can be read and written in place.
    if beta == 0.0 {
        for v in y.iter_mut() {
            *v = 0.0;
        }
    } else if beta != 1.0 {
        for v in y.iter_mut() {
            *v *= beta;
        }
    }

    for j in 0..n {
        let xj = x[j];
        if xj == 0.0 {
            continue;
        }
        let c = alpha * xj;
        // Column j occupies a[j*m .. (j+1)*m].
        let col = &a[j * m..(j + 1) * m];
        for i in 0..m {
            y[i] += c * col[i];
        }
    }
}

/// `y := alpha * A^T * x + beta * y` (dgemv, trans = 'T').
///
/// `a` is `m x n` column-major with `m` rows, so `A^T` is `n x m`: `x` has `m`
/// entries and `y` has `n`.
///
/// Implemented directly rather than through [`dgemv`], because transposing into a
/// temporary would cost a full copy of the matrix for what is a change of loop
/// order.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn dgemv_trans(alpha: f64, a: &[f64], x: &[f64], beta: f64, y: &mut [f64], m: usize, n: usize) {
    assert_eq!(a.len(), m * n, "dgemv_trans: a must be m*n");
    assert_eq!(x.len(), m, "dgemv_trans: x must have m entries");
    assert_eq!(y.len(), n, "dgemv_trans: y must have n entries");

    for j in 0..n {
        // Row j of A^T is column j of A, so this is a dot product per column.
        let col = &a[j * m..(j + 1) * m];
        let mut sum = 0.0;
        for i in 0..m {
            sum += col[i] * x[i];
        }
        y[j] = if beta == 0.0 { alpha * sum } else { alpha * sum + beta * y[j] };
    }
}

/// `A := alpha * x * y^T + A` (dger), with `A` `m x n` column-major.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn dger(alpha: f64, x: &[f64], y: &[f64], a: &mut [f64], m: usize, n: usize) {
    assert_eq!(a.len(), m * n, "dger: a must be m*n");
    assert_eq!(x.len(), m, "dger: x must have m entries");
    assert_eq!(y.len(), n, "dger: y must have n entries");

    for j in 0..n {
        let c = alpha * y[j];
        if c == 0.0 {
            continue;
        }
        let col = &mut a[j * m..(j + 1) * m];
        for i in 0..m {
            col[i] += c * x[i];
        }
    }
}

/// Whether an operand is used as stored or transposed.
///
/// BLAS spells this with a character argument. Taking it as a type here keeps the
/// decision out of the packing loops' inner bodies (each variant's branch is
/// taken once per block) and makes a mistyped flag a compile error rather than a
/// silently different product.
#[derive(Clone, Copy, PartialEq, Eq)]
enum Trans {
    No,
    Yes,
}

/// `C := alpha * A * B + beta * C` (dgemm, both trans = 'N').
///
/// `a` is `m x k`, `b` is `k x n`, `c` is `m x n`, all column-major.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn dgemm(alpha: f64, a: &[f64], b: &[f64], beta: f64, c: &mut [f64], m: usize, n: usize, k: usize) {
    assert_eq!(a.len(), m * k, "dgemm: a must be m*k");
    assert_eq!(b.len(), k * n, "dgemm: b must be k*n");
    assert_eq!(c.len(), m * n, "dgemm: c must be m*n");
    gemm_core(alpha, a, b, beta, c, m, n, k, Trans::No, Trans::No);
}

/// `C := alpha * A^T * B + beta * C` (dgemm, transa = 'T', transb = 'N').
///
/// `a` is `k x m`, `b` is `k x n`, `c` is `m x n`, all column-major. This is the
/// operation `A^T B`, and it is a kernel of its own rather than a transpose into
/// a temporary because transposing an `m x k` operand costs
/// memory traffic proportional to `m*k` for a product that does `m*n*k` flops:
/// harmless when `n` is large and dominant when it is not. See `gemm_core` for
/// where the packing differs.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn dgemm_tn(alpha: f64, a: &[f64], b: &[f64], beta: f64, c: &mut [f64], m: usize, n: usize, k: usize) {
    assert_eq!(a.len(), k * m, "dgemm_tn: a must be k*m");
    assert_eq!(b.len(), k * n, "dgemm_tn: b must be k*n");
    assert_eq!(c.len(), m * n, "dgemm_tn: c must be m*n");
    gemm_core(alpha, a, b, beta, c, m, n, k, Trans::Yes, Trans::No);
}

/// `C := alpha * A * B^T + beta * C` (dgemm, transa = 'N', transb = 'T').
///
/// `a` is `m x k`, `b` is `n x k`, `c` is `m x n`, all column-major. Note the
/// asymmetry with [`dgemm_tn`]: there the contracted dimension is the *row* count
/// of both operands, here it is the *column* count. Handing `dgemm_tn` shapes to
/// this one is silently a different product, not an error.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn dgemm_nt(alpha: f64, a: &[f64], b: &[f64], beta: f64, c: &mut [f64], m: usize, n: usize, k: usize) {
    assert_eq!(a.len(), m * k, "dgemm_nt: a must be m*k");
    assert_eq!(b.len(), n * k, "dgemm_nt: b must be n*k");
    assert_eq!(c.len(), m * n, "dgemm_nt: c must be m*n");
    gemm_core(alpha, a, b, beta, c, m, n, k, Trans::No, Trans::Yes);
}

/// `C := alpha * A^T * B^T + beta * C` (dgemm, transa = 'T', transb = 'T').
///
/// `a` is `k x m`, `b` is `n x k`, `c` is `m x n`, all column-major, so the
/// product is `(m x k) * (k x n)` and both operands are read transposed.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn dgemm_tt(alpha: f64, a: &[f64], b: &[f64], beta: f64, c: &mut [f64], m: usize, n: usize, k: usize) {
    assert_eq!(a.len(), k * m, "dgemm_tt: a must be k*m");
    assert_eq!(b.len(), n * k, "dgemm_tt: b must be n*k");
    assert_eq!(c.len(), m * n, "dgemm_tt: c must be m*n");
    gemm_core(alpha, a, b, beta, c, m, n, k, Trans::Yes, Trans::Yes);
}

/// The blocked kernel behind the four real `dgemm` variants.
///
/// `trans_a` and `trans_b` change only *where a packed element is read from*: the
/// panels, the block sizes and the micro-kernel are the same, which is what makes
/// one kernel worth carrying instead of a temporary transpose plus [`dgemm`].
///
/// The published shape is always `C` `m x n` computed as `A_op (m x k) * B_op
/// (k x n)`. What transposition changes is the shape *as stored*: `A` is `k x m`
/// when transposed and `m x k` otherwise, and likewise `B`. So the packed element
/// `ap[i + p*MC]` is
///
/// ```text
///   not transposed:  a[ic + i + (pc + p) * m]    -- a column of A, contiguous in i
///   transposed:      a[pc + p + (ic + i) * k]    -- a row of A, contiguous in p
/// ```
///
/// and that is the whole difference; the second form reads the *transpose* of the
/// block, which is exactly what is wanted. Both orientations are packed into the
/// same layout, so the inner loop keeps its contiguous vector load.
///
/// # Panics
///
/// Panics if a slice does not match the shape implied by `m`, `n`, `k`.
#[allow(clippy::too_many_arguments)]
fn gemm_core(alpha: f64, a: &[f64], b: &[f64], beta: f64, c: &mut [f64], m: usize, n: usize, k: usize, trans_a: Trans, trans_b: Trans) {
    // The stored shape is `m x k` or `k x m` depending on `trans_a`, but the
    // *length* is `m*k` either way -- the same reason the flag cannot be inferred
    // from the shapes at the call site.
    assert!(a.len() >= m * k, "gemm: a is too short");
    assert!(b.len() >= k * n, "gemm: b is too short");
    assert!(c.len() >= m * n, "gemm: c is too short");

    if m == 0 || n == 0 {
        return;
    }

    // Scale C once, up front. The micro-kernel then only accumulates, so `beta`
    // does not have to reach inside the hot loop.
    match beta {
        0.0 => c.fill(0.0),
        1.0 => {}
        _ => {
            for x in c.iter_mut() {
                *x *= beta;
            }
        }
    }
    if k == 0 || alpha == 0.0 {
        return;
    }

    // Row counts of the operands as they are *stored*, which is also their
    // column-major leading dimension.
    let ada = if trans_a == Trans::No { m } else { k };
    let adb = if trans_b == Trans::No { k } else { n };

    // Block sizes. The packed panels are the only large allocation; everything
    // inside them is reused across the whole inner loop, which is the point --
    // the naive kernel re-reads `a` and `b` from memory for every element of
    // `c`, so it is bound by cache misses rather than by arithmetic.
    const MC: usize = 128;
    const KC: usize = 256;
    const NC: usize = 256;

    let mut ap = vec![0.0f64; MC * KC];
    let mut bp = vec![0.0f64; KC * NC];

    #[cfg(all(target_arch = "x86_64", not(target_arch = "wasm32")))]
    let use_avx = has_avx2_fma();
    #[cfg(not(all(target_arch = "x86_64", not(target_arch = "wasm32"))))]
    let use_avx = false;

    for jc in (0..n).step_by(NC) {
        let nc = NC.min(n - jc);
        for pc in (0..k).step_by(KC) {
            let kc = KC.min(k - pc);

            // Pack B's panel: bp[p + j*kc] = B_op[p][j] of the block at (pc, jc).
            match trans_b {
                // Not transposed, B's column jc+j is contiguous in p, so this is
                // a straight copy -- the inner-loop direction is already the
                // memory direction.
                Trans::No => {
                    for j in 0..nc {
                        let b_col = &b[pc + (jc + j) * adb..pc + kc + (jc + j) * adb];
                        let dst = &mut bp[j * kc..j * kc + kc];
                        dst.copy_from_slice(b_col);
                    }
                }
                // Transposed, B is `n x k`, so row `p` of `B^T` is row jc+j's
                // *column* p read down B -- contiguous in j, not in p. The loop
                // therefore runs p outer and reads B's row in one run, writing
                // bp with stride kc. Reads are the expensive direction, so this
                // keeps them sequential and pays the stride on the store.
                Trans::Yes => {
                    for p in 0..kc {
                        for j in 0..nc {
                            bp[j * kc + p] = b[jc + j + (pc + p) * adb];
                        }
                    }
                }
            }

            for ic in (0..m).step_by(MC) {
                let mc = MC.min(m - ic);
                // Pack A's panel: ap[i + p*MC] = A_op[i][p] of the block at
                // (ic, pc), contiguous in i so the micro-kernel's vector load is
                // a straight 32-byte read.
                //
                // `alpha` is folded in here rather than in the kernels: the
                // micro-kernel keeps C in registers and accumulates, so applying
                // a scale per element there would defeat the point. Scaling the
                // packed panel costs one multiply per element of A per block and
                // leaves both kernels as pure accumulators.
                match trans_a {
                    Trans::No => {
                        for p in 0..kc {
                            let a_col = &a[ic + (pc + p) * ada..ic + mc + (pc + p) * ada];
                            let dst = &mut ap[p * MC..p * MC + mc];
                            if alpha == 1.0 {
                                dst.copy_from_slice(a_col);
                            } else {
                                for (d, s) in dst.iter_mut().zip(a_col) {
                                    *d = alpha * s;
                                }
                            }
                        }
                    }
                    Trans::Yes => {
                        // A is `k x m`: `A^T[i][p]` is A's row ic+i, entry pc+p,
                        // so a fixed i reads contiguously in p. Same trade as
                        // B's: sequential reads, strided stores into the panel.
                        for i in 0..mc {
                            for p in 0..kc {
                                let v = a[pc + p + (ic + i) * ada];
                                ap[p * MC + i] = if alpha == 1.0 { v } else { alpha * v };
                            }
                        }
                    }
                }
                gemm_panel(&ap, &bp, &mut c[ic + jc * m..], m, mc, nc, kc, MC, kc, use_avx);
            }
        }
    }
}

/// Whether this CPU has AVX2 and FMA, checked once.
///
/// The micro-kernel is compiled with `#[target_feature(enable = "avx2,fma")]`,
/// so calling it on a CPU without them is undefined behavior. A build with
/// `-C target-feature=+avx2` would make this redundant, but the crate cannot
/// assume that -- it is built for whatever target the caller chooses, and
/// x86-64 baseline does not include AVX2.
#[cfg(all(target_arch = "x86_64", not(target_arch = "wasm32")))]
fn has_avx2_fma() -> bool {
    use std::sync::OnceLock;
    static CACHED: OnceLock<bool> = OnceLock::new();
    *CACHED.get_or_init(|| std::arch::is_x86_feature_detected!("avx2") && std::arch::is_x86_feature_detected!("fma"))
}

/// Accumulates one packed panel product into `c`, whose element `(i, j)` is at
/// `c[i + j * ldc]` with the panel starting at row `ic` and column `jc`.
///
/// Split out so the safe wrapper above stays readable and the `unsafe` block is
/// confined to one place.
#[allow(clippy::too_many_arguments)]
fn gemm_panel(ap: &[f64], bp: &[f64], c: &mut [f64], ldc: usize, mc: usize, nc: usize, kc: usize, lda: usize, ldb: usize, use_avx: bool) {
    // The panel's leading dimensions are the packed strides, so the micro-kernel
    // reads `ap[i + p*lda]` and `bp[p + j*ldb]`.
    let mut j = 0;
    while j < nc {
        let mut i = 0;
        if j + 4 <= nc && use_avx {
            #[cfg(all(target_arch = "x86_64", not(target_arch = "wasm32")))]
            {
                while i + 4 <= mc {
                    // Safety: interior indices are in bounds by construction
                    // (i + 4 <= mc <= lda rows, j + 4 <= nc). The loads are
                    // unaligned: `ap` starts at i, which is a multiple of 4 but
                    // not necessarily 32-byte aligned.
                    unsafe {
                        gemm_micro_4x4(
                            ap.as_ptr().add(i),
                            lda,
                            bp.as_ptr().add(j * ldb),
                            ldb,
                            c.as_mut_ptr().add(i + j * ldc),
                            ldc,
                            kc,
                        );
                    }
                    i += 4;
                }
            }
            // Remainder rows of this 4-wide column strip.
            while i < mc {
                for jj in j..j + 4 {
                    let mut sum = 0.0;
                    for p in 0..kc {
                        sum += ap[i + p * lda] * bp[p + jj * ldb];
                    }
                    c[i + jj * ldc] += sum;
                }
                i += 1;
            }
            j += 4;
        } else {
            // Fewer than four columns left: no 4x4 block fits.
            for jj in j..nc {
                for ii in 0..mc {
                    let mut sum = 0.0;
                    for p in 0..kc {
                        sum += ap[ii + p * lda] * bp[p + jj * ldb];
                    }
                    c[ii + jj * ldc] += sum;
                }
            }
            j = nc;
        }
    }
}

/// The 4x4 AVX2+FMA micro-kernel: four rows by four columns of the packed
/// panels, accumulated in four vector registers.
///
/// Register pressure is the reason for 4x4. Four `c` accumulators plus one `a`
/// vector is five of the sixteen YMM registers, so the loop does not spill; the
/// `b` scalars are folded in as memory operands. Each iteration of `p` performs
/// one 32-byte load of `a` and four broadcasts, against 32 flops -- so this is
/// arithmetic-bound, unlike the naive kernel.
///
/// # Safety
///
/// Requires AVX2 and FMA. `ap` must have `4 + (kc-1)*lda` elements in reach,
/// `bp` must have `4 + (kc-1)*ldb`, and `c` must have `4 + 3*ldc` from the
/// pointer, with `kc >= 1`.
#[cfg(all(target_arch = "x86_64", not(target_arch = "wasm32")))]
#[target_feature(enable = "avx2,fma")]
unsafe fn gemm_micro_4x4(ap: *const f64, lda: usize, bp: *const f64, ldb: usize, c: *mut f64, ldc: usize, kc: usize) {
    use core::arch::x86_64::{_mm256_add_pd, _mm256_broadcast_sd, _mm256_fmadd_pd, _mm256_loadu_pd, _mm256_setzero_pd, _mm256_storeu_pd};

    // The whole body is unsafe: raw-pointer reads on the packed panels and
    // the AVX2/FMA intrinsics. Written as one block so the pre-2024 `unsafe
    // fn` allowance does not hide which operations are actually unsafe.
    //
    // SAFETY: the function's `# Safety` section states the contract -- `ap` and `bp`
    // hold four rows (columns) by `kc` elements at the packed strides `lda`/`ldb`,
    // and `c` holds four rows by four columns at `ldc`. `gemm_panel` is the only
    // caller and enters here only when `i + 4 <= mc <= lda`, `j + 4 <= nc` and
    // `kc >= 1`, which is that contract.
    unsafe {
        let mut c0 = _mm256_setzero_pd();
        let mut c1 = _mm256_setzero_pd();
        let mut c2 = _mm256_setzero_pd();
        let mut c3 = _mm256_setzero_pd();

        // The four columns of B this block advances. `ldb` is the packed column
        // stride, so column `j` occupies `bp[j*ldb .. j*ldb + kc]` and its p-th
        // element is `bp[j*ldb + p]`. Getting this the wrong way round -- reading
        // `bp[p*ldb + j]`, i.e. B's row p -- silently computes a different product,
        // which is the same wrong-axis mistake this port has produced before.
        let b0 = bp;
        let b1 = bp.add(ldb);
        let b2 = bp.add(2 * ldb);
        let b3 = bp.add(3 * ldb);
        for p in 0..kc {
            let av = _mm256_loadu_pd(ap.add(p * lda));
            c0 = _mm256_fmadd_pd(av, _mm256_broadcast_sd(&*b0.add(p)), c0);
            c1 = _mm256_fmadd_pd(av, _mm256_broadcast_sd(&*b1.add(p)), c1);
            c2 = _mm256_fmadd_pd(av, _mm256_broadcast_sd(&*b2.add(p)), c2);
            c3 = _mm256_fmadd_pd(av, _mm256_broadcast_sd(&*b3.add(p)), c3);
        }

        // Fold the four accumulators into C. Read-modify-write, so that this kernel
        // composes over splits of `k` without needing beta applied per panel.
        let mut t = [0.0f64; 4];
        _mm256_storeu_pd(t.as_mut_ptr(), _mm256_add_pd(_mm256_loadu_pd(c), c0));
        core::ptr::copy_nonoverlapping(t.as_ptr(), c, 4);
        _mm256_storeu_pd(t.as_mut_ptr(), _mm256_add_pd(_mm256_loadu_pd(c.add(ldc)), c1));
        core::ptr::copy_nonoverlapping(t.as_ptr(), c.add(ldc), 4);
        _mm256_storeu_pd(t.as_mut_ptr(), _mm256_add_pd(_mm256_loadu_pd(c.add(2 * ldc)), c2));
        core::ptr::copy_nonoverlapping(t.as_ptr(), c.add(2 * ldc), 4);
        _mm256_storeu_pd(t.as_mut_ptr(), _mm256_add_pd(_mm256_loadu_pd(c.add(3 * ldc)), c3));
        core::ptr::copy_nonoverlapping(t.as_ptr(), c.add(3 * ldc), 4);
    }
}

/// `C := alpha * A * A^T + beta * C` (dsyrk, uplo = 'L' or 'U').
///
/// `a` is `n x k` column-major and `c` is `n x n`. Only the requested triangle of
/// `c` is written; the other is left alone, as BLAS specifies, so a caller that
/// only reads one triangle is not surprised by the rest containing stale data.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn dsyrk(alpha: f64, a: &[f64], beta: f64, c: &mut [f64], n: usize, k: usize, upper: bool) {
    assert_eq!(a.len(), n * k, "dsyrk: a must be n*k");
    assert_eq!(c.len(), n * n, "dsyrk: c must be n*n");

    for j in 0..n {
        // Only the requested triangle of column j. Both branches must be the
        // same type, so the inclusive end is expressed as `j + 1`.
        let i_range = if upper { 0..j + 1 } else { j..n };
        for i in i_range {
            let mut sum = 0.0;
            for l in 0..k {
                sum += a[i + l * n] * a[j + l * n];
            }
            let idx = i + j * n;
            c[idx] = if beta == 0.0 { alpha * sum } else { alpha * sum + beta * c[idx] };
        }
    }
}

/// `nrm := |A|` for the one-norm, infinity-norm or Frobenius norm (dlange).
///
/// `a` is `m x n` column-major. `norm` is the BLAS character: `'1'`, `'I'` (or
/// `'i'`) for infinity, `'F'` (or `'E'`) for Frobenius, and `'M'` for the maximum
/// absolute entry.
///
/// # Panics
///
/// Panics on an unrecognized norm character, which in Fortran would be an
/// argument error and here should be a compile-time-visible mistake.
pub fn dlange(norm: char, a: &[f64], m: usize, n: usize) -> f64 {
    assert_eq!(a.len(), m * n, "dlange: a must be m*n");
    match norm {
        'M' | 'm' => a.iter().fold(0.0f64, |acc, v| acc.max(v.abs())),
        // 'O' is LAPACK's other spelling of the one-norm, and is what the C shim
        // sends for `Norm::One`.
        '1' | 'O' | 'o' => {
            // Maximum column sum.
            (0..n)
                .map(|j| a[j * m..(j + 1) * m].iter().map(|v| v.abs()).sum::<f64>())
                .fold(0.0f64, f64::max)
        }
        'I' | 'i' => {
            // Maximum row sum.
            (0..m).map(|i| (0..n).map(|j| a[i + j * m].abs()).sum::<f64>()).fold(0.0f64, f64::max)
        }
        'F' | 'f' | 'E' | 'e' => dnrm2(a),
        other => panic!("dlange: unrecognized norm '{other}'; expected one of 1, O, I, F, E, M"),
    }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod tests {
    use super::*;

    /// A non-square matrix, so a transposed index cannot accidentally be right.
    /// Column-major storage of
    ///   [ 1  2  3 ]
    ///   [ 4  5  6 ]
    fn a23() -> Vec<f64> {
        vec![1.0, 4.0, 2.0, 5.0, 3.0, 6.0]
    }

    #[test]
    fn level1_operations() {
        let x = vec![1.0, 2.0, 3.0];
        let mut y = vec![10.0, 20.0, 30.0];
        daxpy(2.0, &x, &mut y);
        assert_eq!(y, vec![12.0, 24.0, 36.0]);

        dscal(0.5, &mut y);
        assert_eq!(y, vec![6.0, 12.0, 18.0]);

        let mut z = vec![0.0; 3];
        dcopy(&y, &mut z);
        assert_eq!(z, y);

        assert_eq!(ddot(&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0]), 32.0);
        // Mismatched lengths read the shorter one, as BLAS does.
        assert_eq!(ddot(&[1.0, 2.0, 3.0], &[4.0, 5.0]), 14.0);

        assert_eq!(dnrm2(&[3.0, 4.0]), 5.0);
        assert_eq!(dnrm2(&[]), 0.0);
        assert_eq!(dasum(&[-1.0, 2.0, -3.0]), 6.0);
    }

    #[test]
    fn dnrm2_does_not_overflow_or_underflow() {
        // A naive sum of squares overflows to infinity here.
        let big = vec![1e200, 1e200];
        let n = dnrm2(&big);
        assert!(n.is_finite(), "dnrm2 overflowed: {n}");
        assert!((n / (1e200 * 2.0f64.sqrt()) - 1.0).abs() < 1e-12);

        // ...and underflows to zero here.
        let tiny = vec![1e-200, 1e-200];
        let n = dnrm2(&tiny);
        assert!(n > 0.0, "dnrm2 underflowed to zero");
    }

    #[test]
    fn dgemv_matches_hand_computation() {
        // y := A * [1, 1, 1] over the 2x3 matrix, so y is the row sums: [6, 15].
        let a = a23();
        let x = vec![1.0, 1.0, 1.0];
        let mut y = vec![0.0; 2];
        dgemv(1.0, &a, &x, 0.0, &mut y, 2, 3);
        assert_eq!(y, vec![6.0, 15.0]);

        // beta must scale the existing y.
        let mut y2 = vec![1.0, 1.0];
        dgemv(2.0, &a, &x, 3.0, &mut y2, 2, 3);
        assert_eq!(y2, vec![2.0 * 6.0 + 3.0, 2.0 * 15.0 + 3.0]);
    }

    /// The naive triple loop, used as the reference for the blocked kernel.
    fn naive_dgemm(alpha: f64, a: &[f64], b: &[f64], beta: f64, c: &mut [f64], m: usize, n: usize, k: usize) {
        for j in 0..n {
            for i in 0..m {
                let mut sum = 0.0;
                for l in 0..k {
                    sum += a[i + l * m] * b[l + j * k];
                }
                let idx = i + j * m;
                c[idx] = if beta == 0.0 { alpha * sum } else { alpha * sum + beta * c[idx] };
            }
        }
    }

    /// The blocked kernel must agree with the naive one on every shape.
    ///
    /// The sizes are chosen to land on each branch the kernel has, because the
    /// small hand-written cases above never reach any of them: the 4x4
    /// micro-kernel needs `mc, nc >= 4`; the scalar row and column remainders
    /// need `mc % 4 != 0` and `nc % 4 != 0`; and the packing blocks are MC=128,
    /// KC=256, NC=256, so `m > 128`, `k > 256` or `n > 256` are needed before
    /// the block loop iterates more than once.
    #[test]
    fn dgemm_blocked_matches_naive_on_every_branch() {
        let shapes = [
            // tiny / degenerate
            (1, 1, 1),
            (1, 5, 1),
            (5, 1, 1),
            (1, 1, 5),
            (2, 3, 4),
            // micro-kernel boundaries and remainders
            (3, 3, 3),
            (4, 4, 4),
            (5, 5, 5),
            (7, 9, 11),
            (8, 8, 8),
            (13, 6, 17),
            // block boundaries: MC=128, KC=256, NC=256
            (128, 8, 8),
            (129, 8, 8),
            (8, 128, 8),
            (8, 129, 8),
            (8, 8, 256),
            (8, 8, 257),
            (200, 300, 260),
            (260, 130, 300),
        ];
        for (m, n, k) in shapes {
            let a: Vec<f64> = (0..m * k).map(|i| ((i * 7 % 19) as f64) - 9.0).collect();
            let b: Vec<f64> = (0..k * n).map(|i| ((i * 11 % 23) as f64) - 11.0).collect();
            let c0: Vec<f64> = (0..m * n).map(|i| ((i % 5) as f64) * 0.25).collect();

            for (alpha, beta) in [(1.0, 1.0), (1.0, 0.0), (2.5, -1.5), (0.0, 3.0)] {
                let mut want = c0.clone();
                naive_dgemm(alpha, &a, &b, beta, &mut want, m, n, k);
                let mut got = c0.clone();
                dgemm(alpha, &a, &b, beta, &mut got, m, n, k);
                for idx in 0..m * n {
                    let scale = want[idx].abs().max(1.0);
                    assert!(
                        (got[idx] - want[idx]).abs() / scale < 1e-11,
                        "m={m} n={n} k={k} alpha={alpha} beta={beta}: element {idx} is {} not {}",
                        got[idx],
                        want[idx]
                    );
                }
            }
        }
    }

    #[test]
    fn dgemm_matches_hand_computation() {
        // A (2x3) * B (3x2) = C (2x2), all from a23-like data.
        //   A = [[1,2,3],[4,5,6]]
        //   B = [[7,8],[9,10],[11,12]]  stored column-major as 3x2
        let a = a23();
        let b = vec![7.0, 9.0, 11.0, 8.0, 10.0, 12.0]; // column-major 3x2
        let mut c = vec![0.0; 4];
        dgemm(1.0, &a, &b, 0.0, &mut c, 2, 2, 3);
        // C[0][0] = 1*7+2*9+3*11 = 58 ; C[1][0] = 4*7+5*9+6*11 = 139
        // C[0][1] = 1*8+2*10+3*12 = 64 ; C[1][1] = 4*8+5*10+6*12 = 154
        assert_eq!(c, vec![58.0, 139.0, 64.0, 154.0]);
    }

    #[test]
    fn dgemm_beta_scales_c() {
        let a = a23();
        let b = vec![7.0, 9.0, 11.0, 8.0, 10.0, 12.0];
        let mut c = vec![1.0, 1.0, 1.0, 1.0];
        dgemm(1.0, &a, &b, 2.0, &mut c, 2, 2, 3);
        assert_eq!(c, vec![58.0 + 2.0, 139.0 + 2.0, 64.0 + 2.0, 154.0 + 2.0]);
    }

    #[test]
    fn dsyrk_computes_both_triangles_and_leaves_the_other_alone() {
        // A = [[1,2],[3,4]] as 2x2 column-major: A A^T = [[5,11],[11,25]].
        let a = vec![1.0, 3.0, 2.0, 4.0];
        let n = 2;
        let k = 2;

        let mut lower = vec![0.0; 4];
        dsyrk(1.0, &a, 0.0, &mut lower, n, k, false);
        // Lower triangle: (0,0)=5, (1,0)=11, (1,1)=25; (0,1) untouched.
        assert_eq!(lower[0], 5.0);
        assert_eq!(lower[1], 11.0);
        assert_eq!(lower[3], 25.0);
        assert_eq!(lower[2], 0.0, "the upper triangle must be left alone");

        let mut upper = vec![0.0; 4];
        dsyrk(1.0, &a, 0.0, &mut upper, n, k, true);
        assert_eq!(upper[0], 5.0);
        assert_eq!(upper[2], 11.0);
        assert_eq!(upper[3], 25.0);
        assert_eq!(upper[1], 0.0, "the lower triangle must be left alone");
    }

    #[test]
    fn dlange_norms() {
        let a = a23(); // [[1,2,3],[4,5,6]]
        assert_eq!(dlange('M', &a, 2, 3), 6.0);
        // Column sums: 5, 7, 9 -> max 9.
        assert_eq!(dlange('1', &a, 2, 3), 9.0);
        // Row sums: 6, 15 -> max 15.
        assert_eq!(dlange('I', &a, 2, 3), 15.0);
        // Frobenius: sqrt(1+4+9+16+25+36) = sqrt(91).
        assert!((dlange('F', &a, 2, 3) - 91.0f64.sqrt()).abs() < 1e-12);
        // 'E' is an alias of 'F'.
        assert_eq!(dlange('E', &a, 2, 3), dlange('F', &a, 2, 3));
    }

    #[test]
    #[should_panic(expected = "unrecognized norm")]
    fn dlange_rejects_an_unknown_norm() {
        dlange('X', &[1.0], 1, 1);
    }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Complex BLAS
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// `Complex64` is `#[repr(C)] { re: f64, im: f64 }`, so a slice of it is a
// sequence of interleaved real/imaginary pairs and can be driven by the real
// kernels applied to the two components separately. That is exact rather than an
// approximation: BLAS's complex operations are defined componentwise in the
// same way.
//
// The alternative -- viewing the buffer as `&[f64]` of twice the length and
// treating the complex layout as a stride-2 real matrix -- would work too, but it
// makes the *matrix* stride depend on the complex stride and is much easier to get
// subtly wrong. Adapting the routines below, which already index a matrix, is
// clearer and exercises code that is already tested.

use num_complex::Complex64;

/// `x := alpha * x` for a complex vector (zscal).
///
/// A complex `alpha` mixes the components, so this cannot decompose into two real
/// scalings; it is a genuine complex multiply.
pub fn zscal(alpha: Complex64, x: &mut [Complex64]) {
    for v in x.iter_mut() {
        *v = alpha * *v;
    }
}

/// `y := alpha * x + y` for complex vectors (zaxpy).
///
/// # Panics
///
/// Panics if the slices differ in length.
pub fn zaxpy(alpha: Complex64, x: &[Complex64], y: &mut [Complex64]) {
    assert_eq!(x.len(), y.len(), "zaxpy: vectors must have the same length");
    for i in 0..y.len() {
        y[i] += alpha * x[i];
    }
}

/// `y := x` for complex vectors (zcopy).
///
/// # Panics
///
/// Panics if the slices differ in length.
pub fn zcopy(x: &[Complex64], y: &mut [Complex64]) {
    assert_eq!(x.len(), y.len(), "zcopy: vectors must have the same length");
    y.copy_from_slice(x);
}

/// Euclidean norm of a complex vector (dznrm2), which is real.
pub fn dznrm2(x: &[Complex64]) -> f64 {
    // Sum of squared moduli, scaled to avoid overflow as in dnrm2.
    let mut scale = 0.0f64;
    for v in x {
        let a = v.norm();
        if a > scale {
            scale = a;
        }
    }
    if scale == 0.0 || !scale.is_finite() {
        return x.iter().map(|v| v.norm_sqr()).sum::<f64>().sqrt();
    }
    let inv = 1.0 / scale;
    let ss: f64 = x.iter().map(|v| (v.norm() * inv) * (v.norm() * inv)).sum();
    scale * ss.sqrt()
}

/// `A := alpha * x * y^T + A` for complex vectors (zger, i.e. CBLAS's `zgeru`).
///
/// `A` is `m x n` column-major. **Not conjugated**: BLAS's `zger` is the
/// unconjugated rank-1 update (`zgeru`), while `zgerc` conjugates `y`, and the C
/// shim routes this to `cblas_zgeru` accordingly.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zger(alpha: Complex64, x: &[Complex64], y: &[Complex64], a: &mut [Complex64], m: usize, n: usize) {
    assert_eq!(a.len(), m * n, "zger: a must be m*n");
    assert_eq!(x.len(), m, "zger: x must have m entries");
    assert_eq!(y.len(), n, "zger: y must have n entries");

    for j in 0..n {
        let c = alpha * y[j];
        if c == Complex64::new(0.0, 0.0) {
            continue;
        }
        let col = &mut a[j * m..(j + 1) * m];
        for i in 0..m {
            col[i] += c * x[i];
        }
    }
}

/// `y := alpha * A * x + beta * y` for complex data (zgemv, trans = 'N').
///
/// `a` is `m x n` column-major.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zgemv(alpha: Complex64, a: &[Complex64], x: &[Complex64], beta: Complex64, y: &mut [Complex64], m: usize, n: usize) {
    assert_eq!(a.len(), m * n, "zgemv: a must be m*n");
    assert_eq!(x.len(), n, "zgemv: x must have n entries");
    assert_eq!(y.len(), m, "zgemv: y must have m entries");

    for v in y.iter_mut() {
        *v = beta * *v;
    }
    for j in 0..n {
        let c = alpha * x[j];
        if c == Complex64::new(0.0, 0.0) {
            continue;
        }
        let col = &a[j * m..(j + 1) * m];
        for i in 0..m {
            y[i] += c * col[i];
        }
    }
}

/// `C := alpha * A * B + beta * C` for complex data (zgemm, both 'N').
///
/// `a` is `m x k`, `b` is `k x n`, `c` is `m x n`, all column-major.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zgemm(alpha: Complex64, a: &[Complex64], b: &[Complex64], beta: Complex64, c: &mut [Complex64], m: usize, n: usize, k: usize) {
    assert_eq!(a.len(), m * k, "zgemm: a must be m*k");
    assert_eq!(b.len(), k * n, "zgemm: b must be k*n");
    assert_eq!(c.len(), m * n, "zgemm: c must be m*n");
    zgemm_core(alpha, a, b, beta, c, m, n, k, Trans::No, Trans::No);
}

/// `C := alpha * A^T * B + beta * C` for complex data (zgemm, transa = 'T', transb = 'N').
///
/// `a` is `k x m`, `b` is `k x n`, `c` is `m x n`, all column-major.
///
/// **Not conjugated.** BLAS spells the conjugated operation `'C'` and it is a
/// different product; the complex `A^T B` product passes `CblasTrans`, so this is
/// its twin, not [`zherk`]'s conjugation.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zgemm_tn(alpha: Complex64, a: &[Complex64], b: &[Complex64], beta: Complex64, c: &mut [Complex64], m: usize, n: usize, k: usize) {
    assert_eq!(a.len(), k * m, "zgemm_tn: a must be k*m");
    assert_eq!(b.len(), k * n, "zgemm_tn: b must be k*n");
    assert_eq!(c.len(), m * n, "zgemm_tn: c must be m*n");
    zgemm_core(alpha, a, b, beta, c, m, n, k, Trans::Yes, Trans::No);
}

/// `C := alpha * A * B^T + beta * C` for complex data (zgemm, transa = 'N', transb = 'T').
///
/// `a` is `m x k`, `b` is `n x k`, `c` is `m x n`, all column-major. Not
/// conjugated, as in [`zgemm_tn`].
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zgemm_nt(alpha: Complex64, a: &[Complex64], b: &[Complex64], beta: Complex64, c: &mut [Complex64], m: usize, n: usize, k: usize) {
    assert_eq!(a.len(), m * k, "zgemm_nt: a must be m*k");
    assert_eq!(b.len(), n * k, "zgemm_nt: b must be n*k");
    assert_eq!(c.len(), m * n, "zgemm_nt: c must be m*n");
    zgemm_core(alpha, a, b, beta, c, m, n, k, Trans::No, Trans::Yes);
}

/// `C := alpha * A^T * B^T + beta * C` for complex data (zgemm, transa = 'T', transb = 'T').
///
/// `a` is `k x m`, `b` is `n x k`, `c` is `m x n`, all column-major. Not
/// conjugated, as in [`zgemm_tn`].
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zgemm_tt(alpha: Complex64, a: &[Complex64], b: &[Complex64], beta: Complex64, c: &mut [Complex64], m: usize, n: usize, k: usize) {
    assert_eq!(a.len(), k * m, "zgemm_tt: a must be k*m");
    assert_eq!(b.len(), n * k, "zgemm_tt: b must be n*k");
    assert_eq!(c.len(), m * n, "zgemm_tt: c must be m*n");
    zgemm_core(alpha, a, b, beta, c, m, n, k, Trans::Yes, Trans::Yes);
}

/// The kernel behind the four complex `zgemm` variants.
///
/// Same index algebra as `gemm_core`, without the blocking: this is the
/// straightforward triple loop, and its transposed arms read `a[l + i*k]` (A
/// stored `k x m`) instead of `a[i + l*m]`. The complex entry point is not on a
/// hot path -- the complex `A^T B` product is its only caller -- so it is left
/// unblocked deliberately rather than by omission.
#[allow(clippy::too_many_arguments)]
fn zgemm_core(
    alpha: Complex64,
    a: &[Complex64],
    b: &[Complex64],
    beta: Complex64,
    c: &mut [Complex64],
    m: usize,
    n: usize,
    k: usize,
    trans_a: Trans,
    trans_b: Trans,
) {
    // Row counts as stored, i.e. the leading dimensions.
    let ada = if trans_a == Trans::No { m } else { k };
    let adb = if trans_b == Trans::No { k } else { n };
    for j in 0..n {
        for i in 0..m {
            let mut sum = Complex64::new(0.0, 0.0);
            for l in 0..k {
                let av = match trans_a {
                    Trans::No => a[i + l * ada],
                    Trans::Yes => a[l + i * ada],
                };
                let bv = match trans_b {
                    Trans::No => b[l + j * adb],
                    Trans::Yes => b[j + l * adb],
                };
                sum += av * bv;
            }
            let idx = i + j * m;
            c[idx] = alpha * sum + beta * c[idx];
        }
    }
}

/// `C := alpha * A * A^H + beta * C` for complex data (zherk).
///
/// `a` is `n x k` column-major, `c` is `n x n`. `c` is Hermitian, so only the
/// requested triangle is written and the diagonal is real by construction; the
/// imaginary part of a diagonal entry is set to zero rather than left to rounding.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zherk(alpha: f64, a: &[Complex64], beta: f64, c: &mut [Complex64], n: usize, k: usize, upper: bool) {
    assert_eq!(a.len(), n * k, "zherk: a must be n*k");
    assert_eq!(c.len(), n * n, "zherk: c must be n*n");

    for j in 0..n {
        let i_range = if upper { 0..j + 1 } else { j..n };
        for i in i_range {
            let mut sum = Complex64::new(0.0, 0.0);
            for l in 0..k {
                // a[i][l] * conj(a[j][l])
                sum += a[i + l * n] * a[j + l * n].conj();
            }
            let idx = i + j * n;
            let mut v = Complex64::new(alpha * sum.re, alpha * sum.im) + Complex64::new(beta * c[idx].re, beta * c[idx].im);
            if i == j {
                // Hermitian diagonal must be exactly real.
                v.im = 0.0;
            }
            c[idx] = v;
        }
    }
}

/// `C := alpha * A * A^T + beta * C` for complex data (zsyrk).
///
/// Like [`zherk`] but without conjugation, as BLAS defines `zsyrk` for complex
/// symmetric matrices.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zsyrk(alpha: Complex64, a: &[Complex64], beta: Complex64, c: &mut [Complex64], n: usize, k: usize, upper: bool) {
    assert_eq!(a.len(), n * k, "zsyrk: a must be n*k");
    assert_eq!(c.len(), n * n, "zsyrk: c must be n*n");

    for j in 0..n {
        let i_range = if upper { 0..j + 1 } else { j..n };
        for i in i_range {
            let mut sum = Complex64::new(0.0, 0.0);
            for l in 0..k {
                sum += a[i + l * n] * a[j + l * n];
            }
            let idx = i + j * n;
            c[idx] = alpha * sum + beta * c[idx];
        }
    }
}

/// `|A|` for a complex matrix (zlange), returning a real norm.
///
/// `norm` is `'1'`/`'O'`, `'I'`, `'F'`/`'E'` or `'M'`, as in [`dlange`]. The one- and
/// infinity-norms use the modulus of each entry, as LAPACK specifies.
///
/// # Panics
///
/// Panics on an unrecognized norm character.
pub fn zlange(norm: char, a: &[Complex64], m: usize, n: usize) -> f64 {
    assert_eq!(a.len(), m * n, "zlange: a must be m*n");
    match norm {
        'M' | 'm' => a.iter().fold(0.0f64, |acc, v| acc.max(v.norm())),
        '1' | 'O' | 'o' => {
            // Maximum column sum of moduli. 'O' is LAPACK's other one-norm
            // spelling and is what the C shim sends.
            (0..n)
                .map(|j| a[j * m..(j + 1) * m].iter().map(|v| v.norm()).sum::<f64>())
                .fold(0.0f64, f64::max)
        }
        'I' | 'i' => {
            // Maximum row sum of moduli.
            (0..m).map(|i| (0..n).map(|j| a[i + j * m].norm()).sum::<f64>()).fold(0.0f64, f64::max)
        }
        'F' | 'f' | 'E' | 'e' => {
            // Frobenius: sqrt of the sum of squared moduli.
            let sum: f64 = a.iter().map(|v| v.norm_sqr()).sum();
            sum.sqrt()
        }
        other => panic!("zlange: unrecognized norm '{other}'; expected one of 1, O, I, F, E, M"),
    }
}

/// The symmetric rank-k update BLAS calls `dsyrk`, in the orientation the caller
/// actually has.
///
/// `dsyrk` is specified as `C := alpha A A^T + beta C` with `A` of `(n, k)`, or as
/// `C := alpha A^T A + beta C` with `A` of `(k, n)`. Both are the same
/// computation with the two operands exchanged, so rather than duplicate the
/// kernel this adapts the operands and defers to [`dsyrk`].
///
/// Only the `upper` or lower triangle of `c` is written.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
#[allow(clippy::too_many_arguments)]
pub fn dsyrk_rank(alpha: f64, a: &[f64], beta: f64, c: &mut [f64], n: usize, k: usize, upper: bool, transposed: bool) {
    if !transposed {
        // C := alpha A A^T with A of (n, k): the contraction is over A's columns,
        // which is exactly what `dsyrk` does.
        dsyrk(alpha, a, beta, c, n, k, upper);
        return;
    }

    // C := alpha A^T A with A of (k, n). The *output* is still (n, n), but the
    // contraction runs over A's k rows -- so `k` here is the row count, and
    // delegating to `dsyrk` with the same arguments would contract over the wrong
    // length. This was wrong once: it produced 39 where 15 was expected.
    assert_eq!(a.len(), k * n, "dsyrk_rank: a must be k*n");
    assert_eq!(c.len(), n * n, "dsyrk_rank: c must be n*n");

    for j in 0..n {
        let i_range = if upper { 0..j + 1 } else { j..n };
        for i in i_range {
            let mut sum = 0.0;
            for l in 0..k {
                // a[l][i] * a[l][j], with a of (k, n) column-major.
                sum += a[l + i * k] * a[l + j * k];
            }
            let idx = i + j * n;
            c[idx] = if beta == 0.0 { alpha * sum } else { alpha * sum + beta * c[idx] };
        }
    }
}

/// `y := alpha * A^T * x + beta * y` for complex data (zgemv, trans = 'T').
///
/// `a` is `m x n` column-major, so `A^T` is `n x m`: `x` has `m` entries and `y`
/// has `n`. Not conjugated; BLAS spells that `CblasConjTrans` and it is a
/// different operation.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zgemv_trans(alpha: Complex64, a: &[Complex64], x: &[Complex64], beta: Complex64, y: &mut [Complex64], m: usize, n: usize) {
    assert_eq!(a.len(), m * n, "zgemv_trans: a must be m*n");
    assert_eq!(x.len(), m, "zgemv_trans: x must have m entries");
    assert_eq!(y.len(), n, "zgemv_trans: y must have n entries");

    for j in 0..n {
        let col = &a[j * m..(j + 1) * m];
        let mut sum = Complex64::new(0.0, 0.0);
        for i in 0..m {
            sum += col[i] * x[i];
        }
        y[j] = alpha * sum + beta * y[j];
    }
}

/// `C := alpha * A^T * A + beta * C` for complex data (zsyrk, trans = 'T').
///
/// `a` is `k x n` column-major and `c` is `n x n`, so the contraction runs over
/// `a`'s `k` rows. Not the same call as [`zsyrk`], which contracts over `a`'s
/// columns; passing the same shapes to both is silently wrong.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zsyrk_rank(alpha: Complex64, a: &[Complex64], beta: Complex64, c: &mut [Complex64], n: usize, k: usize, upper: bool) {
    assert_eq!(a.len(), k * n, "zsyrk_rank: a must be k*n");
    assert_eq!(c.len(), n * n, "zsyrk_rank: c must be n*n");

    for j in 0..n {
        let i_range = if upper { 0..j + 1 } else { j..n };
        for i in i_range {
            let mut sum = Complex64::new(0.0, 0.0);
            for l in 0..k {
                sum += a[l + i * k] * a[l + j * k];
            }
            let idx = i + j * n;
            c[idx] = alpha * sum + beta * c[idx];
        }
    }
}

/// `C := alpha * A^H * A + beta * C` for complex data (zherk, trans = 'C').
///
/// `a` is `k x n` column-major, `c` is `n x n` and Hermitian. The diagonal is
/// forced real, as it must be for a Hermitian matrix.
///
/// **Which operand is conjugated matters.** The entry is
/// `sum_l conj(a[l][i]) * a[l][j]`: the conjugation is on the **first** index,
/// which is the one that is contracted with `A^H`'s row. Conjugating the second
/// instead gives the complex conjugate of every off-diagonal entry -- a
/// different Hermitian matrix, and the one this function used to compute. It is
/// a convincing near-miss: the shape and the real parts are right and only the
/// imaginary parts differ, so it takes a genuinely complex test matrix to see it.
/// Verified against `zgemm('C','N')` and `zherk('L','C')` on the test data, where
/// `C(2,1)` must be `15-1j` and the second-index form gives `15+5j`.
///
/// # Panics
///
/// Panics if the shapes are inconsistent.
pub fn zherk_rank(alpha: f64, a: &[Complex64], beta: f64, c: &mut [Complex64], n: usize, k: usize, upper: bool) {
    assert_eq!(a.len(), k * n, "zherk_rank: a must be k*n");
    assert_eq!(c.len(), n * n, "zherk_rank: c must be n*n");

    for j in 0..n {
        let i_range = if upper { 0..j + 1 } else { j..n };
        for i in i_range {
            let mut sum = Complex64::new(0.0, 0.0);
            for l in 0..k {
                // conj(a[l][i]) * a[l][j]: the conjugation is on `i`, not `j`.
                sum += a[l + i * k].conj() * a[l + j * k];
            }
            let idx = i + j * n;
            let mut v = Complex64::new(alpha * sum.re, alpha * sum.im) + Complex64::new(beta * c[idx].re, beta * c[idx].im);
            if i == j {
                // A Hermitian diagonal is real by definition.
                v.im = 0.0;
            }
            c[idx] = v;
        }
    }
}

/// `dtrsv`: solves `A x = b` for a triangular `A`, in place.
///
/// `a` is `n x n`, packed column-major. `lower` selects the triangle that is
/// stored; `unit` declares a unit diagonal, whose stored entries are ignored.
///
/// This is BLAS `DTRSV` without the stride arguments. The transposed solve is
/// [`dtrsv_transpose`], and the matrix right-hand side is [`dtrsm`].
///
/// # Panics
///
/// Panics if `a.len() != n * n` or `x.len() != n`.
pub fn dtrsv(a: &[f64], x: &mut [f64], n: usize, lower: bool, unit: bool) {
    assert_eq!(a.len(), n * n, "dtrsv: a must be n*n");
    assert_eq!(x.len(), n, "dtrsv: x must be n");
    tri_solve_f64(a, x, n, 1, lower, unit, false);
}

/// `dtrsv` with `A` transposed: solves `Aᵀ x = b` in place.
///
/// # Panics
///
/// Panics if `a.len() != n * n` or `x.len() != n`.
pub fn dtrsv_transpose(a: &[f64], x: &mut [f64], n: usize, lower: bool, unit: bool) {
    assert_eq!(a.len(), n * n, "dtrsv_transpose: a must be n*n");
    assert_eq!(x.len(), n, "dtrsv_transpose: x must be n");
    tri_solve_f64(a, x, n, 1, lower, unit, true);
}

/// `dtrsm`: solves a triangular system with a matrix right-hand side.
///
/// With `left`, `B := alpha A⁻¹ B` and `a` is `m x m`; otherwise
/// `B := alpha B A⁻¹` and `a` is `n x n`. `b` is `m x n`, packed column-major.
/// `lower` and `unit` have the same meaning as in [`dtrsv`].
///
/// This is BLAS `DTRSM` for the no-transpose, one-triangle case.
///
/// # Panics
///
/// Panics if `b.len() != m * n`, or if `a` does not have the order the chosen
/// side requires.
pub fn dtrsm(alpha: f64, a: &[f64], b: &mut [f64], m: usize, n: usize, left: bool, lower: bool, unit: bool) {
    assert_eq!(b.len(), m * n, "dtrsm: b must be m*n");
    if left {
        assert_eq!(a.len(), m * m, "dtrsm: a must be m*m for the left side");
        for j in 0..n {
            let x = &mut b[j * m..(j + 1) * m];
            for v in x.iter_mut() {
                *v *= alpha;
            }
            tri_solve_f64(a, x, m, 1, lower, unit, false);
        }
    } else {
        assert_eq!(a.len(), n * n, "dtrsm: a must be n*n for the right side");
        for i in 0..m {
            let x = &mut b[i..];
            for k in 0..n {
                x[k * m] *= alpha;
            }
            tri_solve_f64(a, x, n, m, lower, unit, true);
        }
    }
}

/// The triangular solve behind [`dtrsv`] and [`dtrsm`].
///
/// `x` is read and written with `stride`; `trans` swaps the triangle that the
/// forward/backward sweep walks.
fn tri_solve_f64(a: &[f64], x: &mut [f64], n: usize, stride: usize, lower: bool, unit: bool, trans: bool) {
    let effective_lower = if trans { !lower } else { lower };
    if effective_lower {
        for i in 0..n {
            let mut sum = x[i * stride];
            for k in 0..i {
                let a_ik = if trans { a[k + i * n] } else { a[i + k * n] };
                sum -= a_ik * x[k * stride];
            }
            if !unit {
                sum /= a[i + i * n];
            }
            x[i * stride] = sum;
        }
    } else {
        for i in (0..n).rev() {
            let mut sum = x[i * stride];
            for k in i + 1..n {
                let a_ik = if trans { a[k + i * n] } else { a[i + k * n] };
                sum -= a_ik * x[k * stride];
            }
            if !unit {
                sum /= a[i + i * n];
            }
            x[i * stride] = sum;
        }
    }
}

/// `ztrsv`: solves `A x = b` for a complex triangular `A`, in place.
///
/// Complex twin of [`dtrsv`].
///
/// # Panics
///
/// Panics if `a.len() != n * n` or `x.len() != n`.
pub fn ztrsv(a: &[Complex64], x: &mut [Complex64], n: usize, lower: bool, unit: bool) {
    assert_eq!(a.len(), n * n, "ztrsv: a must be n*n");
    assert_eq!(x.len(), n, "ztrsv: x must be n");
    tri_solve_z(a, x, n, 1, lower, unit, false);
}

/// `ztrsv` with `A` transposed: solves `Aᵀ x = b` in place.
///
/// Complex twin of [`dtrsv_transpose`]. Note that `Aᵀ` is not `Aᴴ`; this is the
/// plain transpose, as in BLAS `'T'`.
///
/// # Panics
///
/// Panics if `a.len() != n * n` or `x.len() != n`.
pub fn ztrsv_transpose(a: &[Complex64], x: &mut [Complex64], n: usize, lower: bool, unit: bool) {
    assert_eq!(a.len(), n * n, "ztrsv_transpose: a must be n*n");
    assert_eq!(x.len(), n, "ztrsv_transpose: x must be n");
    tri_solve_z(a, x, n, 1, lower, unit, true);
}

/// `ztrsm`: complex triangular solve with a matrix right-hand side.
///
/// Complex twin of [`dtrsm`]; no conjugation is applied.
///
/// # Panics
///
/// Panics if `b.len() != m * n`, or if `a` does not have the order the chosen
/// side requires.
pub fn ztrsm(alpha: Complex64, a: &[Complex64], b: &mut [Complex64], m: usize, n: usize, left: bool, lower: bool, unit: bool) {
    assert_eq!(b.len(), m * n, "ztrsm: b must be m*n");
    if left {
        assert_eq!(a.len(), m * m, "ztrsm: a must be m*m for the left side");
        for j in 0..n {
            let x = &mut b[j * m..(j + 1) * m];
            for v in x.iter_mut() {
                *v *= alpha;
            }
            tri_solve_z(a, x, m, 1, lower, unit, false);
        }
    } else {
        assert_eq!(a.len(), n * n, "ztrsm: a must be n*n for the right side");
        for i in 0..m {
            let x = &mut b[i..];
            for k in 0..n {
                x[k * m] *= alpha;
            }
            tri_solve_z(a, x, n, m, lower, unit, true);
        }
    }
}

/// The complex triangular solve behind [`ztrsv`] and [`ztrsm`].
fn tri_solve_z(a: &[Complex64], x: &mut [Complex64], n: usize, stride: usize, lower: bool, unit: bool, trans: bool) {
    let effective_lower = if trans { !lower } else { lower };
    if effective_lower {
        for i in 0..n {
            let mut sum = x[i * stride];
            for k in 0..i {
                let a_ik = if trans { a[k + i * n] } else { a[i + k * n] };
                sum -= a_ik * x[k * stride];
            }
            if !unit {
                sum /= a[i + i * n];
            }
            x[i * stride] = sum;
        }
    } else {
        for i in (0..n).rev() {
            let mut sum = x[i * stride];
            for k in i + 1..n {
                let a_ik = if trans { a[k + i * n] } else { a[i + k * n] };
                sum -= a_ik * x[k * stride];
            }
            if !unit {
                sum /= a[i + i * n];
            }
            x[i * stride] = sum;
        }
    }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod complex_tests {
    use super::*;

    /// A `2 x 4` genuinely complex matrix: with real or Hermitian data the
    /// conjugation order in `zherk_rank` makes no difference, so this has to be
    /// complex to test anything at all.
    fn complex_2x4() -> Vec<Complex64> {
        vec![
            Complex64::new(1.0, -1.0),
            Complex64::new(3.0, 1.0), // column 0
            Complex64::new(2.0, 0.0),
            Complex64::new(1.0, 0.0), // column 1
            Complex64::new(1.0, 0.0),
            Complex64::new(3.0, 0.0), // column 2
            Complex64::new(1.0, 0.0),
            Complex64::new(1.0, 2.0), // column 3
        ]
    }

    fn assert_close(got: Complex64, want: Complex64, what: &str) {
        assert!((got - want).norm() < 1e-13, "{what}: got {got} but the reference says {want}");
    }

    /// `C := A^H A` lower triangle, compared entry by entry with reference BLAS.
    ///
    /// The expected values are the lower triangle of `zgemm('C','N')` on the same
    /// data, cross-checked against `zherk('L','C')` with `alpha=1, beta=0`.
    /// Column-major lower triangle, column by column: (0,0), (1,0), (2,0), (3,0),
    /// (1,1), (2,1), (3,1), (2,2), (3,2), (3,3).
    #[test]
    fn zherk_rank_conjugates_the_first_index() {
        let a = complex_2x4();
        let (k, n) = (2usize, 4usize);
        let mut c = vec![Complex64::new(0.0, 0.0); n * n];
        zherk_rank(1.0, &a, 0.0, &mut c, n, k, false);

        #[rustfmt::skip]
        let lower = [
            (0, 0, 12.0,  0.0),
            (1, 0,  5.0, -1.0),
            (2, 0, 10.0,  2.0),
            (3, 0,  6.0, -6.0),
            (1, 1,  5.0,  0.0),
            (2, 1,  5.0,  0.0),
            (3, 1,  3.0, -2.0),
            (2, 2, 10.0,  0.0),
            (3, 2,  4.0, -6.0),
            (3, 3,  6.0,  0.0),
        ];
        for (i, j, re, im) in lower {
            assert_close(c[i + j * n], Complex64::new(re, im), &format!("({i},{j})"));
        }

        // The upper triangle must be untouched.
        for j in 1..n {
            for i in 0..j {
                assert_eq!(c[i + j * n], Complex64::new(0.0, 0.0), "({i},{j}) must be left alone");
            }
        }

        // Conjugating the second index instead would give the complex conjugate
        // of every off-diagonal entry; this is the assertion that fails if the
        // conjugation order is ever "fixed" back.
        // (row 1, column 0) is `1 + 0 * n`; written out plainly because the
        // `0 * n` term is what `erasing_op` (rightly) objects to.
        assert_eq!(c[1].im, -1.0, "the off-diagonal imaginary part must not be negated");
    }

    /// The same product in the upper triangle, which writes the mirror image.
    #[test]
    fn zherk_rank_upper_triangle() {
        let a = complex_2x4();
        let (k, n) = (2usize, 4usize);
        let mut c = vec![Complex64::new(0.0, 0.0); n * n];
        zherk_rank(1.0, &a, 0.0, &mut c, n, k, true);
        assert_close(c[0 + 1 * n], Complex64::new(5.0, 1.0), "(0,1)");
        assert_close(c[1 + 2 * n], Complex64::new(5.0, 0.0), "(1,2)");
        // The lower triangle must be untouched -- (row 1, column 0) is `c[1]`.
        assert_eq!(c[1], Complex64::new(0.0, 0.0));
    }

    /// `beta` must scale the existing `c`, and the diagonal must come out exactly
    /// real even when `beta * c` has an imaginary part there.
    #[test]
    fn zherk_rank_applies_beta_and_clears_the_diagonal() {
        let a = complex_2x4();
        let (k, n) = (2usize, 4usize);
        let mut c = vec![Complex64::new(0.0, 0.0); n * n];
        for v in c.iter_mut() {
            *v = Complex64::new(1.0, 1.0);
        }
        zherk_rank(2.0, &a, 3.0, &mut c, n, k, false);
        // C(0,0) = 2*12 + 3*(1+1j) -> real 27, imaginary forced to zero.
        assert_close(c[0], Complex64::new(27.0, 0.0), "(0,0)");
        // C(1,0) = 2*(5-1j) + 3*(1+1j) = 13 + 1j.
        assert_close(c[1], Complex64::new(13.0, 1.0), "(1,0)");
    }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/// The transposed `dgemm`/`zgemm` kernels against values captured from CBLAS.
///
/// # Why these are not compared against the other kernels in this file
///
/// The tests above compare a kernel with a reference this crate wrote, which
/// catches a kernel that drifted but not a mistake shared by both. These compare
/// with an *external* implementation.
///
/// # Provenance
///
/// A throwaway C driver (its source and full output are in the commit message
/// that added these kernels) built these same deterministic, non-repeating input
/// matrices and called `cblas_dgemm` in each of the three transposed
/// combinations, and `cblas_zgemm` likewise, with `alpha = 2.5, beta = -1.5` and
/// with the complex pair `alpha = 2 - i, beta = 0.5 + 0.25i`. Every matrix is
/// column-major. The driver printed its *inputs* as well as its outputs, so what
/// is embedded below is a transcription of output rather than of source.
///
/// ```text
/// gcc -O2 -o ref_openblas ref.c -I/usr/include/openblas -fopenmp -lopenblas
/// source /opt/intel/oneapi/setvars.sh
/// gcc -O2 -o ref_mkl ref.c -I$MKLROOT/include -L$MKLROOT/lib/intel64 \
///     -lmkl_intel_lp64 -lmkl_core -lmkl_sequential -lpthread -lm
/// ```
///
/// OpenBLAS and MKL 2026.0 agree on this data: of the 50 printed arrays, 37 are
/// bit-for-bit identical, and the worst disagreement in the remaining 13 is
/// 1.5e-14 relative (in `dgemm_tt`, 473.55555555555287 against 473.5555555555456).
/// *That measured spread is the reason the comparison below is 1e-13 relative
/// and not exact*: it is a real, bounded library-to-library difference, while an
/// index or leading-dimension mistake is wrong by O(1) and cannot hide under it.
/// The small-integer case is exact in both libraries and is asserted exactly.
#[cfg(test)]
mod reference_tests {
    use super::*;
    use crate::cpx;

    /// 1e-13 relative, scaled by the magnitude of the reference value. The
    /// reference libraries themselves disagree by up to 1.5e-14 on this data (see
    /// the module comment), so this is tight enough to fail on a wrong matrix and
    /// loose enough not to fail on a different summation order.
    fn assert_close(got: f64, want: f64, what: &str) {
        let scale = want.abs().max(1.0);
        assert!((got - want).abs() / scale < 1e-13, "{what}: got {got} but the reference says {want}");
    }

    /// The complex twin, comparing the modulus of the difference.
    fn assert_close_cpx(got: Complex64, want: Complex64, what: &str) {
        let scale = want.norm().max(1.0);
        assert!((got - want).norm() / scale < 1e-13, "{what}: got {got} but the reference says {want}");
    }

    /// `C := alpha A^T B + beta C`, A is `5 x 3`, B is `5 x 4`, C is `3 x 4`.
    ///
    /// The shapes are all distinct in magnitude, so a leading dimension taken from
    /// the wrong axis reads a different matrix rather than the same one shifted.
    #[test]
    fn dgemm_tn_matches_cblas() {
        let a: &[f64] = &[
            -266.66666666666669,
            298.33333333333331,
            -21.0,
            286.66666666666669,
            -39.666666666666664,
            104.66666666666667,
            -117.33333333333333,
            -180.33333333333334,
            -221.66666666666666,
            100.0,
            243.33333333333334,
            103.33333333333333,
            -14.0,
            4.333333333333333,
            188.66666666666666,
        ];

        let b: &[f64] = &[
            39.666666666666664,
            -154.33333333333334,
            -316.66666666666669,
            -141.33333333333334,
            -16.666666666666668,
            252.66666666666666,
            272.0,
            -70.666666666666671,
            -96.666666666666671,
            -67.666666666666671,
            -315.66666666666669,
            -222.0,
            -9.0,
            -296.0,
            -27.0,
            -152.33333333333334,
            -323.33333333333331,
            -23.333333333333332,
            219.0,
            -93.666666666666671,
        ];

        let c0: &[f64] = &[
            188.0,
            -183.33333333333334,
            -303.66666666666669,
            -45.333333333333336,
            69.0,
            190.0,
            -173.0,
            20.666666666666668,
            162.0,
            -106.33333333333333,
            198.0,
            -211.0,
        ];

        let want: &[f64] = &[
            -89825.0,
            109028.0,
            -5619.1111111111113,
            -9774.111111111124,
            21935.777777777781,
            77392.888888888891,
            -65645.555555555533,
            57544.555555555547,
            -106002.88888888891,
            11146.555555555573,
            -31710.333333333336,
            -86875.0,
        ];

        let want_ab: &[f64] = &[
            -224844.5,
            272845.0,
            -13592.277777777777,
            -24367.27777777781,
            54735.944444444453,
            193197.22222222222,
            -163854.38888888882,
            143830.38888888888,
            -265250.22222222225,
            28025.888888888934,
            -79572.833333333343,
            -216871.0,
        ];

        let (m, n, k) = (3usize, 4usize, 5usize);

        let mut c = c0.to_vec();
        dgemm_tn(1.0, a, b, 0.0, &mut c, m, n, k);
        for i in 0..m * n {
            assert_close(c[i], want[i], &format!("dgemm_tn [{i}]"));
        }

        let mut c = c0.to_vec();
        dgemm_tn(2.5, a, b, -1.5, &mut c, m, n, k);
        for i in 0..m * n {
            assert_close(c[i], want_ab[i], &format!("dgemm_tn alpha/beta [{i}]"));
        }
    }

    /// `C := alpha A B^T + beta C`, A is `3 x 5`, B is `4 x 5`, C is `3 x 4`.
    #[test]
    fn dgemm_nt_matches_cblas() {
        let a: &[f64] = &[
            242.66666666666666,
            55.0,
            287.66666666666669,
            144.0,
            -254.33333333333334,
            284.0,
            105.0,
            78.666666666666671,
            -87.666666666666671,
            -48.333333333333336,
            -295.33333333333331,
            -227.33333333333334,
            80.0,
            -188.33333333333334,
            -198.0,
        ];

        let b: &[f64] = &[
            113.66666666666667,
            26.0,
            23.333333333333332,
            240.33333333333334,
            221.0,
            221.33333333333334,
            49.333333333333336,
            -219.66666666666666,
            171.33333333333334,
            -87.0,
            -59.0,
            -216.33333333333334,
            179.0,
            -14.0,
            -133.33333333333334,
            -96.666666666666671,
            289.66666666666669,
            -286.0,
            114.33333333333333,
            -228.0,
        ];

        let c0: &[f64] = &[
            -15.333333333333334,
            -3.0,
            -241.0,
            -330.66666666666669,
            306.66666666666669,
            -118.66666666666667,
            271.33333333333331,
            -128.33333333333334,
            -237.0,
            -126.0,
            -212.0,
            72.0,
        ];

        let want: &[f64] = &[
            91918.777777777766,
            -143896.33333333334,
            -17604.777777777785,
            6843.0000000000036,
            -3708.4444444444521,
            137775.66666666669,
            22162.333333333336,
            1939.8888888888862,
            33568.333333333336,
            -9593.8888888888869,
            123557.55555555555,
            92835.333333333343,
        ];

        let want_ab: &[f64] = &[
            229819.94444444441,
            -359736.33333333337,
            -43650.44444444446,
            17603.500000000007,
            -9731.1111111111295,
            344617.16666666674,
            54998.833333333343,
            5042.2222222222154,
            84276.333333333343,
            -23795.722222222219,
            309211.88888888888,
            231980.33333333337,
        ];

        let (m, n, k) = (3usize, 4usize, 5usize);

        let mut c = c0.to_vec();
        dgemm_nt(1.0, a, b, 0.0, &mut c, m, n, k);
        for i in 0..m * n {
            assert_close(c[i], want[i], &format!("dgemm_nt [{i}]"));
        }

        let mut c = c0.to_vec();
        dgemm_nt(2.5, a, b, -1.5, &mut c, m, n, k);
        for i in 0..m * n {
            assert_close(c[i], want_ab[i], &format!("dgemm_nt alpha/beta [{i}]"));
        }
    }

    /// `C := alpha A^T B^T + beta C`, A is `5 x 3`, B is `4 x 5`, C is `3 x 4`.
    #[test]
    fn dgemm_tt_matches_cblas() {
        let a: &[f64] = &[
            -192.33333333333334,
            -188.33333333333334,
            319.0,
            1.3333333333333333,
            198.0,
            73.333333333333329,
            327.0,
            -329.66666666666669,
            -231.0,
            80.666666666666671,
            -167.0,
            109.0,
            -103.33333333333333,
            286.33333333333331,
            82.333333333333329,
        ];

        let b: &[f64] = &[
            -321.33333333333331,
            -217.33333333333334,
            54.666666666666664,
            97.666666666666671,
            6.333333333333333,
            -266.66666666666669,
            271.66666666666669,
            39.333333333333336,
            28.0,
            -235.33333333333334,
            69.666666666666671,
            120.0,
            -4.0,
            70.666666666666671,
            147.0,
            -207.33333333333334,
            -71.0,
            -222.33333333333334,
            200.66666666666666,
            233.0,
        ];

        let c0: &[f64] = &[
            59.0,
            177.66666666666666,
            99.0,
            51.0,
            -122.66666666666667,
            -150.33333333333334,
            -228.66666666666666,
            -277.33333333333331,
            -246.0,
            -145.33333333333334,
            45.0,
            77.333333333333329,
        ];

        let want: &[f64] = &[
            55479.0,
            -35527.444444444445,
            44468.666666666664,
            -26976.444444444445,
            -59815.111111111109,
            33474.555555555555,
            473.55555555555287,
            52107.222222222219,
            71896.0,
            57945.222222222219,
            47153.555555555555,
            -64605.777777777774,
        ];

        let want_ab: &[f64] = &[
            138609.0,
            -89085.111111111109,
            111023.16666666666,
            -67517.611111111109,
            -149353.77777777778,
            83911.888888888891,
            1526.8888888888821,
            130684.05555555555,
            180109.0,
            145081.05555555556,
            117816.38888888889,
            -161630.44444444444,
        ];

        let (m, n, k) = (3usize, 4usize, 5usize);

        let mut c = c0.to_vec();
        dgemm_tt(1.0, a, b, 0.0, &mut c, m, n, k);
        for i in 0..m * n {
            assert_close(c[i], want[i], &format!("dgemm_tt [{i}]"));
        }

        let mut c = c0.to_vec();
        dgemm_tt(2.5, a, b, -1.5, &mut c, m, n, k);
        for i in 0..m * n {
            assert_close(c[i], want_ab[i], &format!("dgemm_tt alpha/beta [{i}]"));
        }
    }

    /// The one case that must match **exactly**: small integer data, `alpha = 1`,
    /// `beta = 0`, so every partial sum is an integer below 2^53 and no rounding
    /// order can change the result. Both reference libraries print these values
    /// bit for bit, and so must the kernel.
    ///
    /// It is worth having both this and the tolerant comparison: an exact case
    /// catches a kernel that is "close but subtly reordered" and, more to the
    /// point, it is the only kind of assertion that would notice a kernel dropping
    /// a term from an integer-magnitude sum.
    #[test]
    fn dgemm_tn_matches_cblas_exactly_on_integers() {
        let a: &[f64] = &[-3.0, 0.0, 3.0, -1.0, 2.0, -2.0, 1.0, -3.0, 0.0, 3.0, -1.0, 2.0, -2.0, 1.0, -3.0];

        let b: &[f64] = &[
            -4.0, 1.0, -3.0, 2.0, -2.0, 3.0, -1.0, 4.0, 0.0, -4.0, 1.0, -3.0, 2.0, -2.0, 3.0, -1.0, 4.0, 0.0, -4.0, 1.0,
        ];

        let want: &[f64] = &[-3.0, 12.0, 20.0, -5.0, -31.0, -1.0, 11.0, -2.0, -22.0, 9.0, 9.0, 2.0];

        let (m, n, k) = (3usize, 4usize, 5usize);
        let mut c = vec![0.0; m * n];
        dgemm_tn(1.0, a, b, 0.0, &mut c, m, n, k);
        assert_eq!(c.as_slice(), want, "dgemm_tn must be exact on integer data");
    }

    /// `C := alpha A^T B + beta C` for complex data, A is `4 x 2`, B is `4 x 3`.
    ///
    /// Genuinely complex everywhere, since with real data the transpose and the
    /// conjugate transpose are the same matrix and nothing distinguishes them.
    #[test]
    fn zgemm_tn_matches_cblas() {
        let a: &[Complex64] = &[
            cpx!(-266.66666666666669, -59.166666666666664),
            cpx!(298.33333333333331, 134.66666666666666),
            cpx!(-21.0, -4.0),
            cpx!(286.66666666666669, -142.0),
            cpx!(-39.666666666666664, -115.66666666666667),
            cpx!(104.66666666666667, 21.0),
            cpx!(-117.33333333333333, -86.333333333333329),
            cpx!(-180.33333333333334, 94.166666666666671),
        ];

        let b: &[Complex64] = &[
            cpx!(-192.33333333333334, -123.66666666666667),
            cpx!(-188.33333333333334, 120.16666666666667),
            cpx!(319.0, -136.16666666666666),
            cpx!(1.3333333333333333, -94.0),
            cpx!(198.0, -72.833333333333329),
            cpx!(73.333333333333329, -149.0),
            cpx!(327.0, 24.5),
            cpx!(-329.66666666666669, 139.83333333333334),
            cpx!(-231.0, 148.0),
            cpx!(80.666666666666671, -127.5),
            cpx!(-167.0, 163.16666666666666),
            cpx!(109.0, -132.16666666666666),
        ];

        let c0: &[Complex64] = &[
            cpx!(-160.66666666666666, 60.666666666666664),
            cpx!(-108.66666666666667, 13.75),
            cpx!(27.333333333333332, 71.916666666666671),
            cpx!(48.833333333333336, 36.0),
            cpx!(3.1666666666666665, -63.583333333333336),
            cpx!(-133.33333333333334, 71.0),
        ];

        let want: &[Complex64] = &[
            cpx!(-68015.111111111095, 107181.2222222222),
            cpx!(-97738.770833333314, 152039.93055555556),
            cpx!(-134964.09027777781, 213039.76388888893),
            cpx!(-112307.58333333333, -247393.70833333334),
            cpx!(147381.42361111112, -346457.55555555562),
            cpx!(166229.97222222222, 13265.805555555547),
        ];

        let (m, n, k) = (2usize, 3usize, 4usize);
        let alpha = cpx!(2.0, -1.0);
        let beta = cpx!(0.5, 0.25);
        let mut c = c0.to_vec();
        zgemm_tn(alpha, a, b, beta, &mut c, m, n, k);
        for i in 0..m * n {
            assert_close_cpx(c[i], want[i], &format!("zgemm_tn [{i}]"));
        }
    }

    /// `C := alpha A B^T + beta C` for complex data, A is `2 x 4`, B is `3 x 4`.
    #[test]
    fn zgemm_nt_matches_cblas() {
        let a: &[Complex64] = &[
            cpx!(39.666666666666664, 94.0),
            cpx!(-154.33333333333334, -91.666666666666671),
            cpx!(-316.66666666666669, -151.83333333333334),
            cpx!(-141.33333333333334, -22.666666666666668),
            cpx!(-16.666666666666668, 34.5),
            cpx!(252.66666666666666, 95.0),
            cpx!(272.0, -86.5),
            cpx!(-70.666666666666671, 10.333333333333334),
        ];

        let b: &[Complex64] = &[
            cpx!(113.66666666666667, -7.666666666666667),
            cpx!(26.0, -1.5),
            cpx!(23.333333333333332, -120.5),
            cpx!(240.33333333333334, -165.33333333333334),
            cpx!(221.0, 153.33333333333334),
            cpx!(221.33333333333334, -59.333333333333336),
            cpx!(49.333333333333336, 135.66666666666666),
            cpx!(-219.66666666666666, -64.166666666666671),
            cpx!(171.33333333333334, -118.5),
            cpx!(-87.0, -63.0),
            cpx!(-59.0, -106.0),
            cpx!(-216.33333333333334, 36.0),
        ];

        let c0: &[Complex64] = &[
            cpx!(-160.66666666666666, 60.666666666666664),
            cpx!(-108.66666666666667, 13.75),
            cpx!(27.333333333333332, 71.916666666666671),
            cpx!(48.833333333333336, 36.0),
            cpx!(3.1666666666666665, -63.583333333333336),
            cpx!(-133.33333333333334, 71.0),
        ];

        let want: &[Complex64] = &[
            cpx!(-245210.5, 162737.38888888891),
            cpx!(-48025.215277777796, 151967.04166666666),
            cpx!(-239712.14583333334, -155012.73611111112),
            cpx!(-211119.41666666669, -41983.347222222219),
            cpx!(-223677.24305555562, 159284.97222222225),
            cpx!(45678.305555555562, -19420.722222222219),
        ];

        let (m, n, k) = (2usize, 3usize, 4usize);
        let alpha = cpx!(2.0, -1.0);
        let beta = cpx!(0.5, 0.25);
        let mut c = c0.to_vec();
        zgemm_nt(alpha, a, b, beta, &mut c, m, n, k);
        for i in 0..m * n {
            assert_close_cpx(c[i], want[i], &format!("zgemm_nt [{i}]"));
        }
    }

    /// `C := alpha A^T B^T + beta C` for complex data, A is `4 x 2`, B is `3 x 4`.
    #[test]
    fn zgemm_tt_matches_cblas() {
        let a: &[Complex64] = &[
            cpx!(-173.0, -96.166666666666671),
            cpx!(-246.0, -94.166666666666671),
            cpx!(-209.66666666666666, 159.5),
            cpx!(193.66666666666666, 0.66666666666666663),
            cpx!(92.0, 99.0),
            cpx!(-329.33333333333331, 36.666666666666664),
            cpx!(216.33333333333334, 163.5),
            cpx!(130.66666666666666, -164.83333333333334),
        ];

        let b: &[Complex64] = &[
            cpx!(59.0, -160.66666666666666),
            cpx!(177.66666666666666, -108.66666666666667),
            cpx!(99.0, 27.333333333333332),
            cpx!(51.0, 48.833333333333336),
            cpx!(-122.66666666666667, 3.1666666666666665),
            cpx!(-150.33333333333334, -133.33333333333334),
            cpx!(-228.66666666666666, 135.83333333333334),
            cpx!(-277.33333333333331, 19.666666666666668),
            cpx!(-246.0, 14.0),
            cpx!(-145.33333333333334, -117.66666666666667),
            cpx!(45.0, 34.833333333333336),
            cpx!(77.333333333333329, 60.0),
        ];

        let c0: &[Complex64] = &[
            cpx!(-160.66666666666666, 60.666666666666664),
            cpx!(-108.66666666666667, 13.75),
            cpx!(27.333333333333332, 71.916666666666671),
            cpx!(48.833333333333336, 36.0),
            cpx!(3.1666666666666665, -63.583333333333336),
            cpx!(-133.33333333333334, 71.0),
        ];

        let want: &[Complex64] = &[
            cpx!(-153416.0, -129677.22222222225),
            cpx!(-237261.88194444444, 62147.180555555577),
            cpx!(76879.965277777781, -111139.90277777778),
            cpx!(-10298.583333333325, -99573.763888888876),
            cpx!(150643.70138888888, -69828.0),
            cpx!(59107.805555555562, -8056.1666666666533),
        ];

        let (m, n, k) = (2usize, 3usize, 4usize);
        let alpha = cpx!(2.0, -1.0);
        let beta = cpx!(0.5, 0.25);
        let mut c = c0.to_vec();
        zgemm_tt(alpha, a, b, beta, &mut c, m, n, k);
        for i in 0..m * n {
            assert_close_cpx(c[i], want[i], &format!("zgemm_tt [{i}]"));
        }
    }

    /// Column-major transpose of a `rows x cols` matrix, into a `cols x rows` one.
    fn transpose(src: &[f64], rows: usize, cols: usize) -> Vec<f64> {
        let mut dst = vec![0.0; rows * cols];
        for j in 0..cols {
            for i in 0..rows {
                dst[j + i * cols] = src[i + j * rows];
            }
        }
        dst
    }

    /// The complex twin of [`transpose`].
    fn transpose_cpx(src: &[Complex64], rows: usize, cols: usize) -> Vec<Complex64> {
        let mut dst = vec![Complex64::new(0.0, 0.0); rows * cols];
        for j in 0..cols {
            for i in 0..rows {
                dst[j + i * cols] = src[i + j * rows];
            }
        }
        dst
    }

    /// `C := alpha A B + beta C` from the definition, with `A` `m x k` and `B`
    /// `k x n`: the plain triple loop, over the already-transposed operands.
    fn naive_mat_mul(alpha: f64, a: &[f64], b: &[f64], beta: f64, c: &mut [f64], m: usize, n: usize, k: usize) {
        for j in 0..n {
            for i in 0..m {
                let mut sum = 0.0;
                for l in 0..k {
                    sum += a[i + l * m] * b[l + j * k];
                }
                let idx = i + j * m;
                c[idx] = if beta == 0.0 { alpha * sum } else { alpha * sum + beta * c[idx] };
            }
        }
    }

    /// The complex twin of [`naive_mat_mul`].
    fn naive_zmat_mul(alpha: Complex64, a: &[Complex64], b: &[Complex64], beta: Complex64, c: &mut [Complex64], m: usize, n: usize, k: usize) {
        for j in 0..n {
            for i in 0..m {
                let mut sum = Complex64::new(0.0, 0.0);
                for l in 0..k {
                    sum += a[i + l * m] * b[l + j * k];
                }
                let idx = i + j * m;
                c[idx] = alpha * sum + beta * c[idx];
            }
        }
    }

    type RealGemm = fn(f64, &[f64], &[f64], f64, &mut [f64], usize, usize, usize);
    type ComplexGemm = fn(Complex64, &[Complex64], &[Complex64], Complex64, &mut [Complex64], usize, usize, usize);

    /// Every real transposed variant over shapes that reach each branch of the
    /// blocking: the 4x4 micro-kernel needs `mc, nc >= 4`, the scalar remainders
    /// need `mc % 4 != 0` and `nc % 4 != 0`, and MC=128, KC=256, NC=256 are the
    /// block boundaries, so `m > 128`, `k > 256` or `n > 256` makes the block loop
    /// repeat.
    ///
    /// This is the shape coverage the fixtures above cannot give: they are a single
    /// 3x4x5 case, which never enters the blocking at all. The reference here
    /// applies the orientation by *actually transposing* the stored operand and
    /// then calling the plain product, so it shares no index algebra with the
    /// kernel under test -- which matters, because a mis-derived reference that did
    /// share it would pass while the kernel was wrong.
    #[test]
    fn transposed_dgemm_matches_the_definition_on_every_branch() {
        let shapes = [
            // tiny / degenerate
            (1, 1, 1),
            (1, 5, 1),
            (5, 1, 1),
            (1, 1, 5),
            (2, 3, 4),
            // micro-kernel boundaries and remainders
            (3, 3, 3),
            (4, 4, 4),
            (5, 5, 5),
            (7, 9, 11),
            (8, 8, 8),
            (13, 6, 17),
            // block boundaries: MC=128, KC=256, NC=256
            (128, 8, 8),
            (129, 8, 8),
            (8, 128, 8),
            (8, 129, 8),
            (8, 8, 256),
            (8, 8, 257),
            (200, 300, 260),
            (260, 130, 300),
        ];
        let variants: [(bool, bool, RealGemm); 3] = [(true, false, dgemm_tn), (false, true, dgemm_nt), (true, true, dgemm_tt)];
        for (m, n, k) in shapes {
            // The operands as the *product* sees them; each variant gets them in
            // the orientation it expects.
            let a_op: Vec<f64> = (0..m * k).map(|i| ((i * 7 % 19) as f64) - 9.0).collect();
            let b_op: Vec<f64> = (0..k * n).map(|i| ((i * 11 % 23) as f64) - 11.0).collect();
            let a_t = transpose(&a_op, m, k);
            let b_t = transpose(&b_op, k, n);
            for (ta, tb, kernel) in variants {
                let a_stored = if ta { &a_t } else { &a_op };
                let b_stored = if tb { &b_t } else { &b_op };
                let c0: Vec<f64> = (0..m * n).map(|i| ((i % 5) as f64) * 0.25).collect();
                for (alpha, beta) in [(1.0, 1.0), (1.0, 0.0), (2.5, -1.5), (0.0, 3.0)] {
                    let mut want = c0.clone();
                    naive_mat_mul(alpha, &a_op, &b_op, beta, &mut want, m, n, k);
                    let mut got = c0.clone();
                    kernel(alpha, a_stored, b_stored, beta, &mut got, m, n, k);
                    for idx in 0..m * n {
                        let scale = want[idx].abs().max(1.0);
                        assert!(
                            (got[idx] - want[idx]).abs() / scale < 1e-11,
                            "ta={ta} tb={tb} m={m} n={n} k={k} alpha={alpha} beta={beta}: element {idx} is {} not {}",
                            got[idx],
                            want[idx]
                        );
                    }
                }
            }
        }
    }

    /// The same sweep for the complex twins, at the sizes where the naive kernel's
    /// index arithmetic branches (it has no blocking to reach).
    #[test]
    fn transposed_zgemm_matches_the_definition_on_every_shape() {
        let shapes = [(1, 1, 1), (1, 4, 1), (3, 2, 5), (2, 3, 4), (5, 5, 5), (7, 11, 13)];
        let variants: [(bool, bool, ComplexGemm); 3] = [(true, false, zgemm_tn), (false, true, zgemm_nt), (true, true, zgemm_tt)];
        let alpha = cpx!(0.5, -1.5);
        let beta = cpx!(2.0, 0.25);
        for (m, n, k) in shapes {
            let a_op: Vec<Complex64> = (0..m * k).map(|i| cpx!((i % 3) as f64 - 1.0, (i % 5) as f64 - 2.0)).collect();
            let b_op: Vec<Complex64> = (0..k * n).map(|i| cpx!((i % 4) as f64 - 1.5, (i % 7) as f64 - 3.0)).collect();
            let a_t = transpose_cpx(&a_op, m, k);
            let b_t = transpose_cpx(&b_op, k, n);
            for (ta, tb, kernel) in variants {
                let a_stored = if ta { &a_t } else { &a_op };
                let b_stored = if tb { &b_t } else { &b_op };
                let c0: Vec<Complex64> = (0..m * n).map(|i| cpx!((i % 2) as f64, 0.5)).collect();
                let mut want = c0.clone();
                naive_zmat_mul(alpha, &a_op, &b_op, beta, &mut want, m, n, k);
                let mut got = c0.clone();
                kernel(alpha, a_stored, b_stored, beta, &mut got, m, n, k);
                for idx in 0..m * n {
                    assert_close_cpx(got[idx], want[idx], &format!("ta={ta} tb={tb} m={m} n={n} k={k} [{idx}]"));
                }
            }
        }
    }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/// The routines the hand-written cases above do not reach, each checked against
/// its mathematical definition rather than against a captured fixture.
///
/// The definitions are written as plain indexed loops with the layout spelled out
/// (`data[i + j * m]`), so they share no code with the kernels they judge. Every
/// shape is non-square and every datum is distinct, which is what makes an index,
/// a contracted length or a leading dimension taken from the wrong axis visible
/// instead of absorbed.
#[cfg(test)]
mod definition_tests {
    use super::*;
    use crate::cpx;

    /// Column-major `m x n` complex matrix whose entries are all distinct and
    /// genuinely complex, so a conjugation error changes the result.
    fn complex_matrix(m: usize, n: usize) -> Vec<Complex64> {
        (0..m * n)
            .map(|idx| {
                let (i, j) = (idx % m, idx / m);
                cpx!(i as f64 + 0.5 * j as f64, 1.0 - 2.0 * i as f64 + 0.25 * j as f64)
            })
            .collect()
    }

    /// Column-major `m x n` real matrix with distinct entries.
    fn real_matrix(m: usize, n: usize) -> Vec<f64> {
        (0..m * n).map(|idx| ((idx * 7 % 19) as f64) - 9.0 + idx as f64 * 0.125).collect()
    }

    fn assert_close(got: f64, want: f64, what: &str) {
        let scale = want.abs().max(1.0);
        assert!((got - want).abs() / scale < 1e-13, "{what}: got {got} but the definition says {want}");
    }

    fn assert_close_cpx(got: Complex64, want: Complex64, what: &str) {
        let scale = want.norm().max(1.0);
        assert!((got - want).norm() / scale < 1e-13, "{what}: got {got} but the definition says {want}");
    }

    /// `zaxpy`, `zscal`, `zcopy`: the complex level-1 updates, on values simple
    /// enough to check by hand as well as by the loop.
    #[test]
    fn complex_level1_operations() {
        let x = [cpx!(1.0, 1.0), cpx!(2.0, -1.0)];
        let mut y = [cpx!(3.0, 0.0), cpx!(0.0, 1.0)];
        zaxpy(cpx!(0.0, 1.0), &x, &mut y);
        // i(1+i) = -1+i and i(2-i) = 1+2i, added to the old y.
        assert_eq!(y, [cpx!(2.0, 1.0), cpx!(1.0, 3.0)]);

        zscal(cpx!(2.0, 0.0), &mut y);
        assert_eq!(y, [cpx!(4.0, 2.0), cpx!(2.0, 6.0)]);

        let mut z = [cpx!(0.0, 0.0); 2];
        zcopy(&y, &mut z);
        assert_eq!(z, y);
    }

    /// `y := alpha * A^T * x + beta * y`, with `a` `2 x 3`: `x` has 2 entries and
    /// `y` has 3, so a leading dimension taken from the wrong axis cannot agree.
    #[test]
    fn dgemv_trans_matches_the_definition() {
        let (m, n) = (2usize, 3usize);
        let a = real_matrix(m, n);
        let x = real_matrix(m, 1);

        for (alpha, beta) in [(1.0, 0.0), (-2.5, 1.75)] {
            let y0 = real_matrix(1, n);
            let mut want = y0.clone();
            for j in 0..n {
                let mut sum = 0.0;
                for i in 0..m {
                    sum += a[i + j * m] * x[i];
                }
                want[j] = alpha * sum + beta * want[j];
            }
            let mut got = y0.clone();
            dgemv_trans(alpha, &a, &x, beta, &mut got, m, n);
            for j in 0..n {
                assert_close(got[j], want[j], &format!("dgemv_trans alpha={alpha} beta={beta} [{j}]"));
            }
        }
    }

    /// `A := alpha * x * y^T + A`, `A` `2 x 3`, including the zero-scalar column
    /// that the kernel skips.
    #[test]
    fn dger_matches_the_definition() {
        let (m, n) = (2usize, 3usize);
        let x = real_matrix(m, 1);
        let y = real_matrix(1, n);
        for (label, y) in [("generic", y.clone()), ("one column zero", vec![y[0], 0.0, y[2]])] {
            let a0 = real_matrix(m, n);
            let mut want = a0.clone();
            for j in 0..n {
                for i in 0..m {
                    want[i + j * m] += 1.5 * x[i] * y[j];
                }
            }
            let mut got = a0.clone();
            dger(1.5, &x, &y, &mut got, m, n);
            for k in 0..m * n {
                assert_close(got[k], want[k], &format!("dger {label} [{k}]"));
            }
        }
    }

    /// Both orientations of the complex matrix-vector product, `a` `2 x 3`.
    #[test]
    fn zgemv_and_zgemv_trans_match_the_definition() {
        let (m, n) = (2usize, 3usize);
        let a = complex_matrix(m, n);
        let alpha = cpx!(0.5, -1.5);
        let beta = cpx!(2.0, 0.25);

        let x = complex_matrix(n, 1);
        let y0 = complex_matrix(m, 1);
        let mut want = y0.clone();
        for i in 0..m {
            let mut sum = Complex64::new(0.0, 0.0);
            for j in 0..n {
                sum += a[i + j * m] * x[j];
            }
            want[i] = alpha * sum + beta * want[i];
        }
        let mut got = y0.clone();
        zgemv(alpha, &a, &x, beta, &mut got, m, n);
        for i in 0..m {
            assert_close_cpx(got[i], want[i], &format!("zgemv [{i}]"));
        }

        // Transposed: x has m entries, y has n.
        let xt = complex_matrix(m, 1);
        let yt0 = complex_matrix(n, 1);
        let mut want = yt0.clone();
        for j in 0..n {
            let mut sum = Complex64::new(0.0, 0.0);
            for i in 0..m {
                sum += a[i + j * m] * xt[i];
            }
            want[j] = alpha * sum + beta * want[j];
        }
        let mut got = yt0.clone();
        zgemv_trans(alpha, &a, &xt, beta, &mut got, m, n);
        for j in 0..n {
            assert_close_cpx(got[j], want[j], &format!("zgemv_trans [{j}]"));
        }
    }

    /// `A := alpha * x * y^T + A` for complex data; not conjugated, so a shim that
    /// routed this to `zgerc` would show in every imaginary part. The second case
    /// zeroes one entry of `y`, which is the arm that skips a whole column.
    #[test]
    fn zger_matches_the_definition() {
        let (m, n) = (2usize, 3usize);
        let x = complex_matrix(m, 1);
        let alpha = cpx!(1.0, -0.5);
        let y = complex_matrix(1, n);
        for (label, y) in [("generic", y.clone()), ("one column zero", vec![y[0], cpx!(0.0, 0.0), y[2]])] {
            let a0 = complex_matrix(m, n);
            let mut want = a0.clone();
            for j in 0..n {
                for i in 0..m {
                    want[i + j * m] += alpha * x[i] * y[j];
                }
            }
            let mut got = a0.clone();
            zger(alpha, &x, &y, &mut got, m, n);
            for k in 0..m * n {
                assert_close_cpx(got[k], want[k], &format!("zger {label} [{k}]"));
            }
        }
    }

    /// `C := alpha * A * B + beta * C` for complex data (both 'N'), `A` `2 x 3` and
    /// `B` `3 x 4`, so all three dimensions differ.
    #[test]
    fn zgemm_matches_the_definition() {
        let (m, n, k) = (2usize, 3usize, 4usize);
        let a = complex_matrix(m, k);
        let b = complex_matrix(k, n);
        let alpha = cpx!(0.5, -1.5);
        let beta = cpx!(-2.0, 0.25);
        let c0 = complex_matrix(m, n);
        let mut want = c0.clone();
        for j in 0..n {
            for i in 0..m {
                let mut sum = Complex64::new(0.0, 0.0);
                for l in 0..k {
                    sum += a[i + l * m] * b[l + j * k];
                }
                let idx = i + j * m;
                want[idx] = alpha * sum + beta * want[idx];
            }
        }
        let mut got = c0.clone();
        zgemm(alpha, &a, &b, beta, &mut got, m, n, k);
        for idx in 0..m * n {
            assert_close_cpx(got[idx], want[idx], &format!("zgemm [{idx}]"));
        }
    }

    /// `zsyrk` (`A A^T`, `a` is `n x k`) and `zsyrk_rank` (`A^T A`, `a` is `k x n`),
    /// both triangles, with the untaken triangle left alone.
    #[test]
    fn zsyrk_routines_match_the_definition() {
        let (n, k) = (3usize, 4usize);
        let alpha = cpx!(2.5, -1.0);
        let beta = cpx!(-1.5, 0.25);
        let c0 = complex_matrix(n, n);

        for transposed in [false, true] {
            let a = complex_matrix(if transposed { k } else { n }, if transposed { n } else { k });
            for upper in [false, true] {
                let mut want = c0.clone();
                for j in 0..n {
                    for i in 0..n {
                        let in_triangle = if upper { i <= j } else { i >= j };
                        if !in_triangle {
                            continue;
                        }
                        let mut sum = Complex64::new(0.0, 0.0);
                        for l in 0..k {
                            sum += if transposed {
                                a[l + i * k] * a[l + j * k]
                            } else {
                                a[i + l * n] * a[j + l * n]
                            };
                        }
                        let idx = i + j * n;
                        want[idx] = alpha * sum + beta * want[idx];
                    }
                }

                let mut got = c0.clone();
                if transposed {
                    zsyrk_rank(alpha, &a, beta, &mut got, n, k, upper);
                } else {
                    zsyrk(alpha, &a, beta, &mut got, n, k, upper);
                }
                for j in 0..n {
                    for i in 0..n {
                        let idx = i + j * n;
                        assert_close_cpx(got[idx], want[idx], &format!("zsyrk transposed={transposed} upper={upper} ({i},{j})"));
                        let in_triangle = if upper { i <= j } else { i >= j };
                        if !in_triangle {
                            assert_eq!(got[idx], c0[idx], "the untaken triangle must be left alone");
                        }
                    }
                }
            }
        }
    }

    /// `dsyrk_rank` with `transposed = true`, which is the `A^T A` path; the
    /// `false` arm forwards to `dsyrk` and is covered by its own test.
    #[test]
    fn dsyrk_rank_transposed_matches_the_definition() {
        let (n, k) = (3usize, 4usize);
        let a = real_matrix(k, n);
        let (alpha, beta) = (2.5, -1.5);
        let c0 = real_matrix(n, n);
        for upper in [false, true] {
            let mut want = c0.clone();
            for j in 0..n {
                for i in 0..n {
                    let in_triangle = if upper { i <= j } else { i >= j };
                    if !in_triangle {
                        continue;
                    }
                    let mut sum = 0.0;
                    for l in 0..k {
                        sum += a[l + i * k] * a[l + j * k];
                    }
                    let idx = i + j * n;
                    want[idx] = alpha * sum + beta * want[idx];
                }
            }
            let mut got = c0.clone();
            dsyrk_rank(alpha, &a, beta, &mut got, n, k, upper, true);
            for idx in 0..n * n {
                assert_close(got[idx], want[idx], &format!("dsyrk_rank upper={upper} [{idx}]"));
            }
            // The forwarding arm must agree with `dsyrk` itself.
            let mut via_rank = c0.clone();
            let mut via_syrk = c0.clone();
            let a_nk = real_matrix(n, k);
            dsyrk_rank(alpha, &a_nk, beta, &mut via_rank, n, k, upper, false);
            dsyrk(alpha, &a_nk, beta, &mut via_syrk, n, k, upper);
            assert_eq!(via_rank, via_syrk, "dsyrk_rank(transposed = false) must forward to dsyrk");
        }
    }

    /// `zherk` (`A A^H`, `a` is `n x k`): the diagonal must come out exactly real
    /// and the untaken triangle must be untouched.
    #[test]
    fn zherk_matches_the_definition() {
        let (n, k) = (3usize, 4usize);
        let a = complex_matrix(n, k);
        let (alpha, beta) = (2.5, -1.5);
        let c0 = complex_matrix(n, n);
        for upper in [false, true] {
            let mut want = c0.clone();
            for j in 0..n {
                for i in 0..n {
                    let in_triangle = if upper { i <= j } else { i >= j };
                    if !in_triangle {
                        continue;
                    }
                    let mut sum = Complex64::new(0.0, 0.0);
                    for l in 0..k {
                        // a[i][l] * conj(a[j][l]): the second index is conjugated.
                        sum += a[i + l * n] * a[j + l * n].conj();
                    }
                    let idx = i + j * n;
                    let mut v = Complex64::new(alpha * sum.re, alpha * sum.im) + Complex64::new(beta * want[idx].re, beta * want[idx].im);
                    if i == j {
                        v.im = 0.0;
                    }
                    want[idx] = v;
                }
            }
            let mut got = c0.clone();
            zherk(alpha, &a, beta, &mut got, n, k, upper);
            for j in 0..n {
                for i in 0..n {
                    let idx = i + j * n;
                    let in_triangle = if upper { i <= j } else { i >= j };
                    if in_triangle {
                        assert_close_cpx(got[idx], want[idx], &format!("zherk upper={upper} ({i},{j})"));
                        if i == j {
                            assert_eq!(got[idx].im, 0.0, "the Hermitian diagonal must be exactly real");
                        }
                    } else {
                        assert_eq!(got[idx], c0[idx], "the untaken triangle must be left alone");
                    }
                }
            }
        }
    }

    /// `zlange` for every norm code, against column sums, row sums, the Frobenius
    /// norm and the maximum modulus; plus its rejection of anything else.
    #[test]
    fn zlange_matches_the_definition() {
        let (m, n) = (2usize, 3usize);
        let a = complex_matrix(m, n);
        let one = (0..n).map(|j| (0..m).map(|i| a[i + j * m].norm()).sum::<f64>()).fold(0.0f64, f64::max);
        let inf = (0..m).map(|i| (0..n).map(|j| a[i + j * m].norm()).sum::<f64>()).fold(0.0f64, f64::max);
        let frob = a.iter().map(|v| v.norm_sqr()).sum::<f64>().sqrt();
        let max = a.iter().map(|v| v.norm()).fold(0.0f64, f64::max);
        for norm in ['1', 'O', 'o', 'I', 'i', 'F', 'f', 'E', 'e', 'M', 'm'] {
            let want = match norm {
                '1' | 'O' | 'o' => one,
                'I' | 'i' => inf,
                'F' | 'f' | 'E' | 'e' => frob,
                _ => max,
            };
            assert_close(zlange(norm, &a, m, n), want, &format!("zlange '{norm}'"));
        }
    }

    #[test]
    #[should_panic(expected = "unrecognized norm")]
    fn zlange_rejects_an_unknown_norm() {
        zlange('X', &[cpx!(1.0, 0.0)], 1, 1);
    }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod triangular_tests {
    use super::*;

    fn assert_close(got: &[f64], want: &[f64], tol: f64) {
        assert_eq!(got.len(), want.len());
        for (i, (g, w)) in got.iter().zip(want.iter()).enumerate() {
            assert!((g - w).abs() <= tol * (1.0 + w.abs()), "entry {i}: {g} vs {w}");
        }
    }

    /// Dense `A x` for a packed column-major `A`.
    fn matvec(a: &[f64], x: &[f64], n: usize) -> Vec<f64> {
        (0..n).map(|i| (0..n).map(|k| a[i + k * n] * x[k]).sum()).collect()
    }

    /// Dense `alpha * A * B` for packed column-major operands.
    fn matmul(alpha: f64, a: &[f64], b: &[f64], m: usize, n: usize, k: usize) -> Vec<f64> {
        let mut c = vec![0.0; m * n];
        for j in 0..n {
            for i in 0..m {
                for l in 0..k {
                    c[i + j * m] += alpha * a[i + l * m] * b[l + j * k];
                }
            }
        }
        c
    }

    /// Lower triangular: `[[2, 0, 0], [1, 3, 0], [4, 5, 6]]`.
    fn lower3() -> Vec<f64> {
        vec![2.0, 1.0, 4.0, 0.0, 3.0, 5.0, 0.0, 0.0, 6.0]
    }

    /// Upper triangular: `[[2, 1, 4], [0, 3, 5], [0, 0, 6]]`.
    fn upper3() -> Vec<f64> {
        vec![2.0, 0.0, 0.0, 1.0, 3.0, 0.0, 4.0, 5.0, 6.0]
    }

    #[test]
    fn dtrsv_solves_lower_and_upper() {
        let x_exact = [1.0, -2.0, 0.5];
        let mut b = matvec(&lower3(), &x_exact, 3);
        dtrsv(&lower3(), &mut b, 3, true, false);
        assert_close(&b, &x_exact, 1e-14);

        let mut b = matvec(&upper3(), &x_exact, 3);
        dtrsv(&upper3(), &mut b, 3, false, false);
        assert_close(&b, &x_exact, 1e-14);
    }

    #[test]
    fn dtrsv_ignores_a_unit_diagonal() {
        // The same lower triangle with a diagonal of 7, declared unit.
        let a = [7.0, 1.0, 4.0, 0.0, 7.0, 5.0, 0.0, 0.0, 7.0];
        let x_exact = [1.0, -2.0, 0.5];
        // Unit diagonal: b = A x uses 1 where the stored entry is 7.
        let unit_a = [1.0, 1.0, 4.0, 0.0, 1.0, 5.0, 0.0, 0.0, 1.0];
        let mut b = matvec(&unit_a, &x_exact, 3);
        dtrsv(&a, &mut b, 3, true, true);
        assert_close(&b, &x_exact, 1e-14);
    }

    #[test]
    fn dtrsv_transpose_solves_the_transposed_system() {
        let x_exact = [1.0, -2.0, 0.5];
        // A is lower, so Aᵀ is upper.
        let at = {
            let a = lower3();
            let mut at = vec![0.0; 9];
            for j in 0..3 {
                for i in 0..3 {
                    at[i + j * 3] = a[j + i * 3];
                }
            }
            at
        };
        let mut b = matvec(&at, &x_exact, 3);
        dtrsv_transpose(&lower3(), &mut b, 3, true, false);
        assert_close(&b, &x_exact, 1e-14);
    }

    #[test]
    fn dtrsm_solves_both_sides() {
        let alpha = 2.0;
        let a = lower3();
        // Left: B = A * X, so the routine returns alpha * A^-1 * B = alpha * X.
        let x_exact = [1.0, 2.0, 3.0, -1.0, 0.5, 2.0];
        let mut b = matmul(1.0, &a, &x_exact, 3, 2, 3);
        dtrsm(alpha, &a, &mut b, 3, 2, true, true, false);
        let want: Vec<f64> = x_exact.iter().map(|v| alpha * v).collect();
        assert_close(&b, &want, 1e-14);

        // Right: B = X * A, so the routine returns alpha * B * A^-1 = alpha * X.
        let x_exact = [1.0, 4.0, 2.0, 5.0, 3.0, 6.0];
        let mut b = matmul(1.0, &x_exact, &a, 2, 3, 3);
        dtrsm(alpha, &a, &mut b, 2, 3, false, true, false);
        let want: Vec<f64> = x_exact.iter().map(|v| alpha * v).collect();
        assert_close(&b, &want, 1e-14);
    }

    #[test]
    fn ztrsv_and_ztrsm_solve_complex_systems() {
        let a = [
            Complex64::new(1.0, 1.0),
            Complex64::new(2.0, -1.0),
            Complex64::new(0.0, 0.0),
            Complex64::new(3.0, 0.5),
        ];
        let x_exact = [Complex64::new(1.0, -2.0), Complex64::new(0.5, 1.0)];
        let mut b = [Complex64::new(0.0, 0.0); 2];
        for i in 0..2 {
            b[i] = (0..2).map(|k| a[i + k * 2] * x_exact[k]).sum();
        }
        ztrsv(&a, &mut b, 2, true, false);
        for i in 0..2 {
            assert!((b[i] - x_exact[i]).norm() < 1e-14);
        }

        // Matrix right-hand side, left solve.
        let alpha = Complex64::new(1.5, -0.5);
        let x_exact = [
            Complex64::new(1.0, 0.0),
            Complex64::new(2.0, 1.0),
            Complex64::new(-1.0, 0.5),
            Complex64::new(0.0, 2.0),
        ];
        let mut b = [Complex64::new(0.0, 0.0); 4];
        for j in 0..2 {
            for i in 0..2 {
                b[i + j * 2] = (0..2).map(|k| a[i + k * 2] * x_exact[k + j * 2]).sum();
            }
        }
        ztrsm(alpha, &a, &mut b, 2, 2, true, true, false);
        for i in 0..4 {
            assert!((b[i] - alpha * x_exact[i]).norm() < 1e-14);
        }
    }

    #[test]
    fn empty_systems_are_noops() {
        dtrsv(&[], &mut [], 0, true, false);
        dtrsv(&[], &mut [], 0, false, false);
        dtrsv_transpose(&[], &mut [], 0, true, false);
        dtrsm(1.0, &[], &mut [], 0, 0, true, true, false);
        ztrsv(&[], &mut [], 0, true, false);
        ztrsm(Complex64::new(1.0, 0.0), &[], &mut [], 0, 0, true, true, false);
    }

    #[test]
    #[should_panic(expected = "dtrsv: x must be n")]
    fn dtrsv_rejects_a_wrong_rhs_length() {
        dtrsv(&[1.0], &mut [1.0, 2.0], 1, true, false);
    }
}