algebra_linear 1.2.0

//! Native real eigen-decompositions: `dgeev` and the pieces it is built from.
//!
//! # Why this is split into steps
//!
//! `dgeev` in LAPACK is a pipeline, and each stage is useful (and testable) on
//! its own:
//!
//! | Step | LAPACK | What it does |
//! |---|---|---|
//! | 1 | `dgebal` | permutes and scales to improve the conditioning |
//! | 2 | `dgehrd` | orthogonal reduction to upper Hessenberg form, `H = Q^T A Q` |
//! | 3 | `dhseqr` | the QR algorithm: real Schur form `T = Z^T H Z` |
//! | 4 | `dtrevc` | eigenvectors of `T` by back-substitution |
//! | 5 | `dgebak` | undoes the balancing |
//!
//! All five steps are here, plus the `dgeev` driver, in the reference's order:
//! balancing, Hessenberg reduction, Schur reduction, eigenvectors, and undoing
//! the balancing. Everything is column-major, so element `(i, j)` of an `n x n`
//! matrix lives at `data[i + j * n]`.
//!
//! # The Schur form is a *real* one
//!
//! A real matrix can have complex eigenvalues, so the real Schur form is only
//! *quasi*-triangular: 1x1 blocks on the diagonal for real eigenvalues and 2x2
//! blocks for complex conjugate pairs. Those 2x2 blocks are why the eigenvector
//! step needs `dlaln2`-style solving rather than a plain triangular substitution.

use super::LapackError;
use super::auxiliary::{Laln2, dlaln2};

/// A 2x2 block's eigenvalues, as [`dlanv2`] returns them.
///
/// `rt1r`/`rt1i` and `rt2r`/`rt2i` are the two eigenvalues. They are real when
/// `rt1i == rt2i == 0`, and otherwise a complex conjugate pair with `rt1i > 0`,
/// which is LAPACK's convention and is relied on by the 2x2-block handling.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Eigen2 {
    /// Real part of the first eigenvalue
    pub rt1r: f64,
    /// Imaginary part of the first eigenvalue (`> 0` for a conjugate pair)
    pub rt1i: f64,
    /// Real part of the second eigenvalue
    pub rt2r: f64,
    /// Imaginary part of the second eigenvalue (`-rt1i` for a pair)
    pub rt2i: f64,
}

/// `dlarfg`: an elementary reflector `H = I - tau * v * v^T` with `H x = beta e1`.
///
/// `alpha` is the entry to annihilate and `x` the rest of the column. The return
/// is `(beta, tau, v)` where `v[0] == 1` and the remaining entries are the tail;
/// `tau == 0` marks the identity (the tail was already zero) and `v` is then all
/// zeros.
///
/// The scaling follows LAPACK 3.12.1's `dlarfg`, because the naive form fails at
/// both ends of the exponent range:
///
/// * `beta = -sign(alpha) * hypot(alpha, ||x||)`. The *opposite* sign to `alpha`
///   avoids catastrophic cancellation when `x` is small; taking the same sign
///   loses most of the significant digits.
/// * If `beta` comes out subnormal then the norm and `beta` are themselves
///   inaccurate, so the vector is scaled up by `safmin` repeatedly (at most 20
///   times) and recomputed. The result is scaled back at the end, so `beta` and
///   `v` are consistent with the *original* input.
pub fn dlarfg(alpha: f64, x: &[f64]) -> (f64, f64, Vec<f64>) {
    let n = x.len() + 1;
    let mut v = vec![0.0; n];
    if n <= 1 {
        return (alpha, 0.0, v);
    }

    // ||x||_2, rescaled so that the sum of squares cannot overflow.
    let mut xscale = 0.0;
    for &xi in x {
        let a = xi.abs();
        if a > xscale {
            xscale = a;
        }
    }
    if xscale == 0.0 {
        // The tail is already zero, so no reflector is needed.
        return (alpha, 0.0, v);
    }
    let ss: f64 = x.iter().map(|&xi| (xi / xscale) * (xi / xscale)).sum();
    let mut xnorm = xscale * ss.sqrt();

    // `SAFMIN = DLAMCH('S') / DLAMCH('E')`, not `DLAMCH('S')`: the reference
    // divides by the relative precision, so the threshold is 1e-292 and not
    // 2.2e-308.
    let safmin = f64::MIN_POSITIVE / f64::EPSILON;
    let mut alpha = alpha;
    let mut beta = -alpha.signum() * f64::hypot(alpha, xnorm);

    if beta.abs() < safmin {
        // `XNORM` and `BETA` may be inaccurate: scale `X` and `ALPHA` *up* by
        // `RSAFMN = 1/SAFMN` and recompute, at most 20 times. Scaling them down
        // instead -- which this used to do -- drives the tail to zero, so `BETA`
        // ends up as `-0` and the reflector is lost.
        let rsafmn = 1.0 / safmin;
        let mut knt = 0usize;
        let mut xw = x.to_vec();
        loop {
            knt += 1;
            for xi in xw.iter_mut() {
                *xi *= rsafmn;
            }
            alpha *= rsafmn;
            // `XNORM = DNRM2(N-1, X, 1)`, rescaled so the squares cannot overflow.
            let mut sc = 0.0f64;
            for &xi in xw.iter() {
                let a = xi.abs();
                if a > sc {
                    sc = a;
                }
            }
            xnorm = if sc == 0.0 {
                0.0
            } else {
                let ss: f64 = xw.iter().map(|&xi| (xi / sc) * (xi / sc)).sum();
                sc * ss.sqrt()
            };
            beta = -alpha.signum() * f64::hypot(alpha, xnorm);
            if beta.abs() >= safmin || knt >= 20 {
                break;
            }
        }

        // TAU and V come out of the scaled frame, in which the frame factor
        // cancels: `(x*R^k) / ((alpha-beta)*R^k) = x/(alpha-beta)`. Only `BETA`
        // has to be taken back, by `SAFMN**KNT` -- multiplied in a loop, because
        // `SAFMN**KNT` itself overflows for KNT > 2.
        let tau = (beta - alpha) / beta;
        let s = 1.0 / (alpha - beta);
        for (i, xi) in xw.iter().enumerate() {
            v[i + 1] = xi * s;
        }
        for _ in 0..knt {
            beta *= safmin;
        }
        let _ = alpha;
        v[0] = 1.0;
        return (beta, tau, v);
    }

    let tau = (beta - alpha) / beta;
    let s = 1.0 / (alpha - beta);
    for (i, xi) in x.iter().enumerate() {
        v[i + 1] = xi * s;
    }
    v[0] = 1.0;
    (beta, tau, v)
}

/// `dlanv2`: the eigenvalues of a 2x2 block, plus the rotation that standardizes it.
///
/// Follows LAPACK 3.12.1's `dlanv2`. On input `a, b, c, d` are
///
/// ```text
/// [ a  b ]
/// [ c  d ]
/// ```
///
/// and `cs`/`sn` are the cosine and sine of the rotation applied. For real
/// eigenvalues the rotation makes the block upper triangular; for a complex pair
/// it makes the two diagonal entries equal, which is what lets the QR algorithm
/// deflate the block. `rt1i > 0` for a pair.
#[allow(clippy::too_many_arguments)]
pub fn dlanv2(a: &mut f64, b: &mut f64, c: &mut f64, d: &mut f64, cs: &mut f64, sn: &mut f64) -> Eigen2 {
    // LAPACK's MULTPL: below 4*eps the decision between real and complex
    // eigenvalues is not trustworthy, so the "equalize the diagonal" branch is
    // taken instead.
    const MULTPL: f64 = 4.0;
    let safmin = f64::MIN_POSITIVE;
    let eps = f64::EPSILON;
    let safmn2 = (safmin / eps).sqrt();
    let safmx2 = 1.0 / safmn2;

    let mut temp;

    if *c == 0.0 {
        // Already upper triangular.
        *cs = 1.0;
        *sn = 0.0;
    } else if *b == 0.0 {
        // Swap rows and columns: the block is lower triangular.
        *cs = 0.0;
        *sn = 1.0;
        temp = *d;
        *d = *a;
        *a = temp;
        *b = -*c;
        *c = 0.0;
    } else if (*a - *d) == 0.0 && b.signum() != c.signum() {
        // Already in standard complex form.
        *cs = 1.0;
        *sn = 0.0;
    } else {
        temp = *a - *d;
        let mut p = 0.5 * temp;
        let bcmax = b.abs().max(c.abs());
        let bcmis = b.abs().min(c.abs()) * b.signum() * c.signum();
        let scale = p.abs().max(bcmax);
        let mut z = (p / scale) * p + (bcmax / scale) * bcmis;

        if z >= MULTPL * eps {
            // Real eigenvalues.
            z = p + z.sqrt() * scale.sqrt() * p.signum();
            *a = *d + z;
            *d -= (bcmax / z) * bcmis;
            let tau = f64::hypot(*c, z);
            *cs = z / tau;
            *sn = *c / tau;
            *b -= *c;
            *c = 0.0;
        } else {
            // Complex eigenvalues, or real (almost) equal ones: equalize the
            // diagonal first.
            let mut count = 0;
            let mut sigma = *b + *c;
            loop {
                count += 1;
                let sc = temp.abs().max(sigma.abs());
                if sc >= safmx2 {
                    sigma *= safmn2;
                    temp *= safmn2;
                } else if sc <= safmn2 {
                    sigma *= safmx2;
                    temp *= safmx2;
                } else {
                    break;
                }
                if count > 20 {
                    break;
                }
            }
            p = 0.5 * temp;
            let tau = f64::hypot(sigma, temp);
            *cs = (0.5 * (1.0 + sigma.abs() / tau)).sqrt();
            *sn = -(p / (tau * *cs)) * sigma.signum();

            // [ AA BB ] = [ A B ] [ CS -SN ]
            // [ CC DD ]   [ C D ] [ SN  CS ]
            let aa = *a * *cs + *b * *sn;
            let bb = -*a * *sn + *b * *cs;
            let cc = *c * *cs + *d * *sn;
            let dd = -*c * *sn + *d * *cs;

            // [ A B ] = [ CS  SN ] [ AA BB ]
            // [ C D ]   [-SN  CS ] [ CC DD ]
            //
            // The products are bracketed the way LAPACK brackets them, to keep a
            // fused multiply-add from changing the result; see the FMA note in
            // dlanv2.f (upstream issue #1031).
            *a = aa * *cs + cc * *sn;
            *b = (bb * *cs) + (dd * *sn);
            *c = -(aa * *sn) + (cc * *cs);
            *d = -bb * *sn + dd * *cs;

            temp = 0.5 * (*a + *d);
            *a = temp;
            *d = temp;

            if *c != 0.0 {
                if *b != 0.0 {
                    if b.signum() == c.signum() {
                        // The pair turned out to be real after all: reduce to
                        // upper triangular form.
                        let sab = b.abs().sqrt();
                        let sac = c.abs().sqrt();
                        p = sab * sac * c.signum();
                        let tau = 1.0 / (b.abs() + c.abs()).sqrt();
                        *a = temp + p;
                        *d = temp - p;
                        *b -= *c;
                        *c = 0.0;
                        let cs1 = sab * tau;
                        let sn1 = sac * tau;
                        temp = *cs * cs1 - *sn * sn1;
                        *sn = *cs * sn1 + *sn * cs1;
                        *cs = temp;
                    }
                } else {
                    *b = -*c;
                    *c = 0.0;
                    temp = *cs;
                    *cs = -*sn;
                    *sn = temp;
                }
            }
        }
    }

    let rt1r = *a;
    let rt2r = *d;
    if *c == 0.0 {
        Eigen2 {
            rt1r,
            rt1i: 0.0,
            rt2r,
            rt2i: 0.0,
        }
    } else {
        let wi = b.abs().sqrt() * c.abs().sqrt();
        Eigen2 {
            rt1r,
            rt1i: wi,
            rt2r,
            rt2i: -wi,
        }
    }
}

/// `dgehrd` (the unblocked `dgehd2` path): reduce to upper Hessenberg form.
///
/// On exit `a` holds `H = Q^T A Q`, which is upper Hessenberg. `q`, if given,
/// receives the orthogonal `Q` accumulated from the reflectors, so that
/// `A = Q H Q^T`.
///
/// The reflectors are the same ones LAPACK packs into the zeroed part of `a`;
/// they are kept in a separate `Vec` here because this caller wants both `H` and
/// `Q`, and the packing buys nothing outside the blocked algorithm.
///
/// # The update is the literal `dgehd2` sequence
///
/// A step is `A := H A H` with `H = I - tau v v^T`, where `H` has order
/// `m = n - i - 1` and acts on rows/columns `i+1 .. n-1` (`v[0] == 1` sits at
/// row `i+1`). LAPACK forms this as two in-place applications against the same
/// array, in this order:
///
/// ```text
/// A(0:n, i+1:n)   := A(0:n, i+1:n) * H      (all rows, columns i+1..)
/// A(i+1:n, i+1:n) := H * A(i+1:n, i+1:n)
/// ```
///
/// The left application deliberately reads the values the right one just wrote,
/// because the target is `H * (A * H)`. Row `i` and column `i` are never touched.
/// Building `H` on rows `i .. n-1` (one row too early) or as a square reflector
/// on the wrong number of columns both break the similarity while still leaving
/// an orthogonal `Q`, so the mistake is easy to miss.
///
/// # Panics
///
/// Panics if `a.len() != n * n` or `q` is given with the wrong length.
pub fn dgehrd(a: &mut [f64], n: usize, q: Option<&mut [f64]>) {
    assert_eq!(a.len(), n * n, "dgehrd: a must be n*n");
    if let Some(qq) = q.as_ref() {
        assert_eq!(qq.len(), n * n, "dgehrd: q must be n*n");
    }
    if n < 3 {
        // Nothing to reduce; Q is the identity.
        if let Some(qq) = q {
            qq.fill(0.0);
            for i in 0..n {
                qq[i + i * n] = 1.0;
            }
        }
        return;
    }

    // A reflector is dropped (tau == 0) when a column is already in the required
    // form, so it is stored as an option.
    let mut vs: Vec<Option<(Vec<f64>, f64)>> = vec![None; n];

    for i in 0..(n - 2) {
        // Annihilate a[i+2..n, i] into a[i+1, i].
        let alpha = a[(i + 1) + i * n];
        let tail: Vec<f64> = ((i + 2)..n).map(|r| a[r + i * n]).collect();
        let (beta, tau, v) = dlarfg(alpha, &tail);
        a[(i + 1) + i * n] = beta;
        for r in (i + 2)..n {
            a[r + i * n] = 0.0;
        }
        if tau == 0.0 {
            continue;
        }
        vs[i] = Some((v.clone(), tau));

        // Apply H(i) from both sides, exactly as dgehd2 does. The reflector has
        // order m = n-i-1 and acts on rows/columns i+1..n-1; v[0] == 1 sits at
        // row i+1. Row i and column i are never touched.
        let m = n - i - 1;

        // Right: A(0:n, i+1:n) := A(0:n, i+1:n) * H, with
        // (A H)(r, j) = A(r, j) - tau * (sum_l A(r, i+1+l) v[l]) * v[j-(i+1)].
        // Every dot product is formed from the block before anything is written.
        let mut w = vec![0.0; n];
        for r in 0..n {
            let mut acc = 0.0;
            for l in 0..m {
                acc += a[r + (i + 1 + l) * n] * v[l];
            }
            w[r] = acc;
        }
        for r in 0..n {
            for l in 0..m {
                a[r + (i + 1 + l) * n] -= tau * w[r] * v[l];
            }
        }

        // Left: A(i+1:n, i+1:n) := H * A(i+1:n, i+1:n). This reads the values
        // the right application just wrote, giving H * (A * H) = H A H.
        for j in (i + 1)..n {
            let mut acc = 0.0;
            for l in 0..m {
                acc += v[l] * a[(i + 1 + l) + j * n];
            }
            for l in 0..m {
                a[(i + 1 + l) + j * n] -= tau * v[l] * acc;
            }
        }
    }

    if let Some(qq) = q {
        // Accumulate Q = H(0) H(1) ... H(n-3) into the identity.
        //
        // Applying the reflectors in reverse and multiplying on the left builds
        // the product: after processing i, `qq` holds H(i) H(i+1) ... H(n-3).
        // The reflector has order m = n-i-1, acts on rows i+1..n-1, and its
        // vector is `v` with v[0] == 1 at row i+1. Each column is read into a
        // scratch buffer first: forming the dot product from `qq` in place would
        // read rows the same update had already changed.
        qq.fill(0.0);
        for i in 0..n {
            qq[i + i * n] = 1.0;
        }
        for i in (0..(n - 2)).rev() {
            let (v, tau) = match &vs[i] {
                Some((v, tau)) => (v.as_slice(), *tau),
                None => continue,
            };
            let m = v.len();
            let base = i + 1;
            let mut col = vec![0.0; m];
            for j in 0..n {
                for r in 0..m {
                    col[r] = qq[(base + r) + j * n];
                }
                let mut acc = 0.0;
                for r in 0..m {
                    acc += v[r] * col[r];
                }
                let w = tau * acc;
                for r in 0..m {
                    qq[(base + r) + j * n] = col[r] - w * v[r];
                }
            }
        }
    }
}

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

/// `dlahqr` (LAPACK 3.12.1, `SRC/dlahqr.f`) with `WANTZ`.
///
/// This is the same literal transcription as before; the only addition is the
/// accumulation of the rotations into `z` when `wantz` is set, which is what
/// `dhseqr` with `COMPZ = 'V'` needs. `z` is `n x n` column-major and is only
/// referenced when `wantz` is true.
fn dlahqr(n: usize, ilo: usize, ihi: usize, h: &mut [f64], ldh: usize, wr: &mut [f64], wi: &mut [f64], wantz: bool, z: &mut [f64]) -> usize {
    // PARAMETER ( KEXSH = 10 ) and the exceptional-shift constants DAT1/DAT2.
    const KEXSH: usize = 10;
    const DAT1: f64 = 3.0 / 4.0;
    const DAT2: f64 = -0.4375;

    // Quick return if possible.
    if n == 0 {
        return 0;
    }
    if ilo == ihi {
        wr[ilo] = h[ilo + ilo * ldh];
        wi[ilo] = 0.0;
        return 0;
    }

    // ==== clear out the trash ====
    // DO 10 J = ILO, IHI - 3
    let mut j = ilo;
    while j + 3 <= ihi {
        h[(j + 2) + j * ldh] = 0.0;
        h[(j + 3) + j * ldh] = 0.0;
        j += 1;
    }
    if ilo + 2 <= ihi {
        h[ihi + (ihi - 2) * ldh] = 0.0;
    }

    let nh = ihi - ilo + 1;

    // Set machine-dependent constants for the stopping criterion.
    let safmin = f64::MIN_POSITIVE;
    let ulp = f64::EPSILON;
    let smlnum = safmin * (nh as f64 / ulp);

    // I1 and I2 are the indices of the first row and last column of H to which
    // transformations must be applied. With WANTT they are the whole matrix and
    // are fixed; otherwise they are set inside the main loop.
    let wantt = true;
    let mut i1 = 0usize;
    let mut i2 = n - 1;

    // ITMAX is the total number of QR iterations allowed.
    let itmax = 30 * nh.max(10);

    // KDEFL counts the number of iterations since a deflation.
    let mut kdefl: usize = 0;

    // The main loop begins here. I is the loop index and decreases from IHI to
    // ILO in steps of 1 or 2. Each iteration works with the active submatrix in
    // rows and columns L to I. Eigenvalues I+1 to IHI have already converged.
    // Either L = ILO or H(L,L-1) is negligible so that the matrix splits.
    let mut i = ihi;
    loop {
        // 20 CONTINUE / IF ( I.LT.ILO ) GO TO 160
        if i < ilo {
            return 0;
        }
        let mut l = ilo;

        // DO 140 ITS = 0, ITMAX. A completed loop falls through to INFO = I.
        let mut split_off = false;
        let mut iters = 0usize;
        while iters <= itmax {
            iters += 1;

            // Look for a single small subdiagonal element.
            // DO 30 K = I, L + 1, -1; running off the top leaves K = L.
            let mut k = i;
            while k > l {
                if h[k + (k - 1) * ldh].abs() <= smlnum {
                    break;
                }
                let mut tst = h[(k - 1) + (k - 1) * ldh].abs() + h[k + k * ldh].abs();
                if tst == 0.0 {
                    if k >= ilo + 2 {
                        tst += h[(k - 1) + (k - 2) * ldh].abs();
                    }
                    if k < ihi {
                        tst += h[(k + 1) + k * ldh].abs();
                    }
                }
                // ==== The following is a conservative small subdiagonal
                // .    deflation criterion due to Ahues & Tisseur (LAWN 122,
                // .    1997). It has better mathematical foundation and
                // .    improves accuracy in some cases. ====
                if h[k + (k - 1) * ldh].abs() <= ulp * tst {
                    let ab = h[k + (k - 1) * ldh].abs().max(h[(k - 1) + k * ldh].abs());
                    let ba = h[k + (k - 1) * ldh].abs().min(h[(k - 1) + k * ldh].abs());
                    let aa = h[k + k * ldh].abs().max((h[(k - 1) + (k - 1) * ldh] - h[k + k * ldh]).abs());
                    let bb = h[k + k * ldh].abs().min((h[(k - 1) + (k - 1) * ldh] - h[k + k * ldh]).abs());
                    let s = aa + ab;
                    if ba * (ab / s) <= smlnum.max(ulp * (bb * (aa / s))) {
                        break;
                    }
                }
                k -= 1;
            }
            // 40 CONTINUE
            l = k;
            if l > ilo {
                // H(L,L-1) is negligible.
                h[l + (l - 1) * ldh] = 0.0;
            }

            // Exit from loop if a submatrix of order 1 or 2 has split off.
            if l + 1 >= i {
                split_off = true;
                break;
            }
            kdefl += 1;

            // Now the active submatrix is in rows and columns L to I. If
            // eigenvalues only are being computed, only the active submatrix
            // need be transformed.
            if !wantt {
                i1 = l;
                i2 = i;
            }

            // Prepare the shift. The two exceptional shifts read a subdiagonal
            // two rows in from an end of the block, which is safe here because
            // L < I - 1 was checked above: the block always has at least 3 rows.
            let (mut h11, mut h12, mut h21, mut h22);
            if kdefl % (2 * KEXSH) == 0 {
                // Exceptional shift (from the bottom).
                let s = h[i + (i - 1) * ldh].abs() + h[(i - 1) + (i - 2) * ldh].abs();
                h11 = DAT1 * s + h[i + i * ldh];
                h12 = DAT2 * s;
                h21 = s;
                h22 = h11;
            } else if kdefl % KEXSH == 0 {
                // Exceptional shift (from the top).
                let s = h[(l + 1) + l * ldh].abs() + h[(l + 2) + (l + 1) * ldh].abs();
                h11 = DAT1 * s + h[l + l * ldh];
                h12 = DAT2 * s;
                h21 = s;
                h22 = h11;
            } else {
                // Prepare to use Francis' double shift
                // (i.e. 2nd degree generalized Rayleigh quotient).
                h11 = h[(i - 1) + (i - 1) * ldh];
                h21 = h[i + (i - 1) * ldh];
                h12 = h[(i - 1) + i * ldh];
                h22 = h[i + i * ldh];
            }
            let s = h11.abs() + h12.abs() + h21.abs() + h22.abs();
            let mut rt1r = 0.0;
            let mut rt1i = 0.0;
            let mut rt2r = 0.0;
            let mut rt2i = 0.0;
            if s == 0.0 {
                // All four entries are zero: the shifts are zero.
            } else {
                h11 /= s;
                h21 /= s;
                h12 /= s;
                h22 /= s;
                let tr = (h11 + h22) / 2.0;
                let det = (h11 - tr) * (h22 - tr) - h12 * h21;
                let rtdisc = det.abs().sqrt();
                if det >= 0.0 {
                    // ==== complex conjugate shifts ====
                    rt1r = tr * s;
                    rt2r = rt1r;
                    rt1i = rtdisc * s;
                    rt2i = -rt1i;
                } else {
                    // ==== real shifts (use only one of them) ====
                    rt1r = tr + rtdisc;
                    rt2r = tr - rtdisc;
                    if (rt1r - h22).abs() <= (rt2r - h22).abs() {
                        rt1r *= s;
                        rt2r = rt1r;
                    } else {
                        rt2r *= s;
                        rt1r = rt2r;
                    }
                    rt1i = 0.0;
                    rt2i = 0.0;
                }
            }

            // Look for two consecutive small subdiagonal elements.
            // DO 50 M = I - 2, L, -1; the M = L case always breaks out at 60.
            let mut m = i - 2;
            let mut v = [0.0f64; 3];
            loop {
                // Determine the effect of starting the double-shift QR
                // iteration at row M, and see if this would make H(M,M-1)
                // negligible. (The following uses scaling to avoid overflows and
                // most underflows.)
                let mut h21s = h[(m + 1) + m * ldh];
                let sv = (h[m + m * ldh] - rt2r).abs() + rt2i.abs() + h21s.abs();
                h21s = h[(m + 1) + m * ldh] / sv;
                v[0] = h21s * h[m + (m + 1) * ldh] + (h[m + m * ldh] - rt1r) * ((h[m + m * ldh] - rt2r) / sv) - rt1i * (rt2i / sv);
                v[1] = h21s * (h[m + m * ldh] + h[(m + 1) + (m + 1) * ldh] - rt1r - rt2r);
                v[2] = h21s * h[(m + 2) + (m + 1) * ldh];
                let sv = v[0].abs() + v[1].abs() + v[2].abs();
                v[0] /= sv;
                v[1] /= sv;
                v[2] /= sv;
                if m == l {
                    break;
                }
                if h[m + (m - 1) * ldh].abs() * (v[1].abs() + v[2].abs())
                    <= ulp * v[0].abs() * (h[(m - 1) + (m - 1) * ldh].abs() + h[m + m * ldh].abs() + h[(m + 1) + (m + 1) * ldh].abs())
                {
                    break;
                }
                m -= 1;
            }
            // 60 CONTINUE

            // Double-shift QR step.
            // DO 130 K = M, I - 1
            for k in m..=(i - 1) {
                // NR is the order of the reflector G, at most 3.
                let nr = 3.min(i - k + 1);
                if k > m {
                    // The subsequent iterations determine the reflector from
                    // the bulge in column K-1.
                    for t in 0..nr {
                        v[t] = h[(k + t) + (k - 1) * ldh];
                    }
                }
                let (beta, t1, vv) = dlarfg(v[0], &v[1..nr]);
                v[..nr].copy_from_slice(&vv[..nr]);
                if k > m {
                    h[k + (k - 1) * ldh] = beta;
                    h[(k + 1) + (k - 1) * ldh] = 0.0;
                    if k < i - 1 {
                        h[(k + 2) + (k - 1) * ldh] = 0.0;
                    }
                } else if m > l {
                    // ==== Use the following instead of
                    // .    H( K, K-1 ) = -H( K, K-1 ) to
                    // .    avoid a bug when v(2) and v(3)
                    // .    underflow. ====
                    h[k + (k - 1) * ldh] *= 1.0 - t1;
                }
                let v2 = v[1];
                let t2 = t1 * v2;
                if nr == 3 {
                    let v3 = v[2];
                    let t3 = t1 * v3;

                    // Apply G from the left to transform the rows of the matrix
                    // in columns K to I2.
                    for jj in k..=i2 {
                        let sum = h[k + jj * ldh] + v2 * h[(k + 1) + jj * ldh] + v3 * h[(k + 2) + jj * ldh];
                        h[k + jj * ldh] -= sum * t1;
                        h[(k + 1) + jj * ldh] -= sum * t2;
                        h[(k + 2) + jj * ldh] -= sum * t3;
                    }

                    // Apply G from the right to transform the columns of the
                    // matrix in rows I1 to min(K+3,I).
                    for jj in i1..=(k + 3).min(i) {
                        let sum = h[jj + k * ldh] + v2 * h[jj + (k + 1) * ldh] + v3 * h[jj + (k + 2) * ldh];
                        h[jj + k * ldh] -= sum * t1;
                        h[jj + (k + 1) * ldh] -= sum * t2;
                        h[jj + (k + 2) * ldh] -= sum * t3;
                    }

                    // Accumulate transformations in the matrix Z.
                    if wantz {
                        for jj in ilo..=ihi {
                            let sum = z[jj + k * n] + v2 * z[jj + (k + 1) * n] + v3 * z[jj + (k + 2) * n];
                            z[jj + k * n] -= sum * t1;
                            z[jj + (k + 1) * n] -= sum * t2;
                            z[jj + (k + 2) * n] -= sum * t3;
                        }
                    }
                } else if nr == 2 {
                    // Apply G from the left to transform the rows of the matrix
                    // in columns K to I2.
                    for jj in k..=i2 {
                        let sum = h[k + jj * ldh] + v2 * h[(k + 1) + jj * ldh];
                        h[k + jj * ldh] -= sum * t1;
                        h[(k + 1) + jj * ldh] -= sum * t2;
                    }

                    // Apply G from the right to transform the columns of the
                    // matrix in rows I1 to I.
                    for jj in i1..=i {
                        let sum = h[jj + k * ldh] + v2 * h[jj + (k + 1) * ldh];
                        h[jj + k * ldh] -= sum * t1;
                        h[jj + (k + 1) * ldh] -= sum * t2;
                    }

                    // Accumulate transformations in the matrix Z.
                    if wantz {
                        for jj in ilo..=ihi {
                            let sum = z[jj + k * n] + v2 * z[jj + (k + 1) * n];
                            z[jj + k * n] -= sum * t1;
                            z[jj + (k + 1) * n] -= sum * t2;
                        }
                    }
                }
            }
        }

        // Failure to converge in remaining number of iterations: INFO = I.
        if !split_off {
            return i;
        }

        // 150 CONTINUE
        if l == i {
            // H(I,I-1) is negligible: one eigenvalue has converged.
            wr[i] = h[i + i * ldh];
            wi[i] = 0.0;
        } else if l + 1 == i {
            // H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
            // Transform the 2-by-2 submatrix to standard Schur form, and compute
            // and store the eigenvalues.
            let (mut a, mut b, mut c, mut d) = (h[(i - 1) + (i - 1) * ldh], h[(i - 1) + i * ldh], h[i + (i - 1) * ldh], h[i + i * ldh]);
            let (mut cs, mut sn) = (0.0, 0.0);
            let e = dlanv2(&mut a, &mut b, &mut c, &mut d, &mut cs, &mut sn);
            h[(i - 1) + (i - 1) * ldh] = a;
            h[(i - 1) + i * ldh] = b;
            h[i + (i - 1) * ldh] = c;
            h[i + i * ldh] = d;
            wr[i - 1] = e.rt1r;
            wi[i - 1] = e.rt1i;
            wr[i] = e.rt2r;
            wi[i] = e.rt2i;

            if wantt {
                // Apply the transformation to the rest of H.
                //
                // DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH, CS, SN )
                if i2 > i {
                    for jj in (i + 1)..=i2 {
                        let x = h[(i - 1) + jj * ldh];
                        let y = h[i + jj * ldh];
                        let temp = cs * x + sn * y;
                        h[i + jj * ldh] = cs * y - sn * x;
                        h[(i - 1) + jj * ldh] = temp;
                    }
                }
                // DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN )
                let count = i - i1 - 1;
                for r in i1..(i1 + count) {
                    let x = h[r + (i - 1) * ldh];
                    let y = h[r + i * ldh];
                    let temp = cs * x + sn * y;
                    h[r + i * ldh] = cs * y - sn * x;
                    h[r + (i - 1) * ldh] = temp;
                }
            }
            if wantz {
                // Apply the transformation to Z.
                //
                // DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
                for r in ilo..=ihi {
                    let x = z[r + (i - 1) * n];
                    let y = z[r + i * n];
                    let temp = cs * x + sn * y;
                    z[r + i * n] = cs * y - sn * x;
                    z[r + (i - 1) * n] = temp;
                }
            }
        }

        // Reset the deflation counter and return to the start of the main loop
        // with the new value of I.
        kdefl = 0;
        if l == 0 {
            return 0;
        }
        i = l - 1;
    }
}

/// `dhseqr` with `JOB = 'S'`, `COMPZ = 'N'`: reduce an upper Hessenberg matrix to
/// real Schur form and return its eigenvalues.
///
/// On entry `t` is `n x n` and upper Hessenberg; on exit it is quasi-upper
/// triangular (zero below the first subdiagonal) and `(wr, wi)` holds the
/// eigenvalues in deflation order. A complex conjugate pair occupies two
/// adjacent positions with `wi > 0` for the first, matching the 2x2 block.
///
/// This calls `dlahqr` directly, which is what the reference `dhseqr` does for
/// `n <= 75`; the blocked `dlaqr0`/`dlaqr5` path is not implemented. `JOB = 'S'`
/// means `WANTT = .TRUE.` and `COMPZ = 'N'` means `WANTZ = .FALSE.`, and with
/// `ILO = 1, IHI = N` the whole matrix is reduced.
///
/// `dlahqr` reports non-convergence through `INFO`, which is returned as
/// [`LapackError::Convergence`]: a matrix that exhausts the iteration limit
/// leaves `t` partly reduced, and the partly converged eigenvalues are not
/// passed off as a Schur form.
///
/// The acceptance test is `parity::dhseqr_eigenvalues_agree_with_the_reference`,
/// which checks the quasi-triangular structure, the standard form of each 2x2
/// block, and the eigenvalue multiset against the reference `dgeev`.
pub fn dhseqr(t: &mut [f64], wr: &mut [f64], wi: &mut [f64], n: usize) -> Result<(), LapackError> {
    dhseqr_impl(t, wr, wi, None, n)
}

/// `dhseqr` with `JOB = 'S'`, `COMPZ = 'V'`: the real Schur form *and* the Schur
/// vectors.
///
/// This is the same reduction as [`dhseqr`], with `WANTZ = .TRUE.` instead of
/// `.FALSE.`, so the plane rotations that `dlahqr` applies are also accumulated
/// into `z`. On entry `z` must be an `n x n` orthogonal matrix (in `dgeev` it is
/// the `Q` from [`dgehrd`], so that `z` starts life as the Hessenberg
/// transformation); on exit it holds `z * Z_schur`, the Schur vectors of the
/// original matrix, so that `A = Z T Z^T` with `T` in `t`.
///
/// `t` and `(wr, wi)` are exactly as in the 4-argument form: the rotations
/// applied to `t` do not depend on whether they are recorded.
///
/// The test `dhseqr_with_schur_vectors_gives_the_schur_decomposition` checks
/// `A = Z T Z^T` and that `Z` is orthogonal.
pub(crate) fn dhseqr_with_schur_vectors(t: &mut [f64], wr: &mut [f64], wi: &mut [f64], z: &mut [f64], n: usize) -> Result<(), LapackError> {
    dhseqr_impl(t, wr, wi, Some(z), n)
}

fn dhseqr_impl(t: &mut [f64], wr: &mut [f64], wi: &mut [f64], z: Option<&mut [f64]>, n: usize) -> Result<(), LapackError> {
    assert_eq!(t.len(), n * n, "dhseqr: t must be n*n");
    assert_eq!(wr.len(), n, "dhseqr: wr must have n entries");
    assert_eq!(wi.len(), n, "dhseqr: wi must have n entries");
    if n == 0 {
        return Ok(());
    }
    // JOB = 'S' -> WANTT = .TRUE.; ILO = 1, IHI = N (0-based 0 and n - 1).
    let info = match z {
        Some(z) => {
            assert_eq!(z.len(), n * n, "dhseqr: z must be n*n");
            dlahqr(n, 0, n - 1, t, n, wr, wi, true, z)
        }
        None => dlahqr(n, 0, n - 1, t, n, wr, wi, false, &mut []),
    };
    // `INFO > 0` is dlahqr's non-convergence report. Returning it is the whole
    // point of the Result: a release build used to drop it on the floor.
    if info != 0 {
        return Err(LapackError::Convergence);
    }
    Ok(())
}

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

/// `dlartg` (LAPACK 3.12.1, `SRC/dlartg.f90`): a plane rotation with real cosine
/// and sine, returning `(c, s, r)` with
///
/// ```text
/// [  c  s ] [ f ]   [ r ]
/// [ -s  c ] [ g ] = [ 0 ]
/// ```
///
/// Used by `dgeev` to rotate the largest component of a complex eigenvector onto
/// the real axis. The three scaling regimes (tiny, normal, huge) are the
/// reference's, so that `f*f + g*g` can neither overflow nor underflow.
pub fn dlartg(f: f64, g: f64) -> (f64, f64, f64) {
    let safmin = f64::MIN_POSITIVE;
    let safmax = f64::MAX;
    let rtmin = safmin.sqrt();
    let rtmax = (safmax / 2.0).sqrt();
    let f1 = f.abs();
    let g1 = g.abs();
    if g == 0.0 {
        (1.0, 0.0, f)
    } else if f == 0.0 {
        (0.0, g.signum(), g1)
    } else if f1 > rtmin && f1 < rtmax && g1 > rtmin && g1 < rtmax {
        let d = (f * f + g * g).sqrt();
        let c = f1 / d;
        let r = d.copysign(f);
        let s = g / r;
        (c, s, r)
    } else {
        let u = safmax.min(safmin.max(f1.max(g1)));
        let fs = f / u;
        let gs = g / u;
        let d = (fs * fs + gs * gs).sqrt();
        let c = fs.abs() / d;
        let mut r = d.copysign(f);
        let s = gs / r;
        r *= u;
        (c, s, r)
    }
}

/// `dlaln2` on the path `dtrevc` needs: no transpose, `ca = 1`, `d1 = d2 = 1`.
///
/// Three of the general routine's arguments are constant here, and naming them
/// once keeps the fifteen call sites to the values that actually vary. This is a
/// forwarder, not a second implementation: the arithmetic lives in
/// [`super::auxiliary::dlaln2`].
#[allow(clippy::too_many_arguments)]
fn dlaln2_real(
    na: usize,
    nw: usize,
    smin: f64,
    a11: f64,
    a21: f64,
    a12: f64,
    a22: f64,
    b11: f64,
    b21: f64,
    b12: f64,
    b22: f64,
    wr: f64,
    wi: f64,
) -> Laln2 {
    dlaln2(false, na, nw, smin, 1.0, a11, a21, a12, a22, 1.0, 1.0, b11, b21, b12, b22, wr, wi)
}

// ---------------------------------------------------------------------------
// 1-based accessors into a column-major n x n array, so the transcription below
// reads like the Fortran.
// ---------------------------------------------------------------------------

fn tget(t: &[f64], n: usize, i: usize, j: usize) -> f64 {
    t[(i - 1) + (j - 1) * n]
}

fn vget(v: &[f64], n: usize, i: usize, j: usize) -> f64 {
    v[(i - 1) + (j - 1) * n]
}

fn vset(v: &mut [f64], n: usize, i: usize, j: usize, x: f64) {
    v[(i - 1) + (j - 1) * n] = x;
}

/// `dtrevc3` (LAPACK 3.12.1, `SRC/dtrevc3.f`) with `SIDE = 'R'`, `HOWMNY = 'B'`
/// and the unblocked `NB = 1` back-transformation.
///
/// On entry `t` is the quasi-upper-triangular Schur form and `vr` holds the
/// orthogonal Schur vectors `Z` (the `VR` of `dgeev`); on exit `vr` holds the
/// right eigenvectors of the original matrix in LAPACK's compact storage: for a
/// conjugate pair, columns `j` and `j + 1` are the real and imaginary parts of
/// `v(j)`, and the eigenvalue with the positive imaginary part comes first.
///
/// The BLAS-2 back-transformation is `dtrevc3`'s `NB = 1` branch, which is what
/// `dgeev.f` selects for its `lwork = 4n + 1` workspace.
/// `VCRIT = BIGNUM` is set again after each rescale, and the reference does the
/// same even though the next statement overwrites it (`VCRIT = BIGNUM / VMAX`):
/// the write is dead there too. Kept so a diff against `dtrevc3.f` still lines
/// up; `unused_assignments` is the only lint it trips.
#[allow(unused_assignments)]
fn dtrevc_right(t: &[f64], vr: &mut [f64], n: usize) {
    let unfl = f64::MIN_POSITIVE;
    let ulp = f64::EPSILON;
    let smlnum = unfl * (n as f64 / ulp);
    let bignum = (1.0 - ulp) / smlnum;

    // 1-norm of the strictly upper triangular part of each column of T, to
    // control overflow in the triangular solver.
    let mut colnorm = vec![0.0; n];
    for j in 2..=n {
        let mut s = 0.0;
        for i in 1..j {
            s += tget(t, n, i, j).abs();
        }
        colnorm[j - 1] = s;
    }

    let mut work1 = vec![0.0; n];
    let mut work2 = vec![0.0; n];
    let mut y1 = vec![0.0; n];
    let mut y2 = vec![0.0; n];

    // IP = 0: real eigenvalue; 1: first of a conjugate pair, already handled;
    // -1: second of a conjugate pair (the KI at which the pair is processed).
    let mut ip = 0i32;
    let mut ki = n; // 1-based, from N down to 1
    while ki >= 1 {
        if ip != 1 {
            // ---- decide whether KI is real or the lower half of a pair ----
            // The reference spells the first two cases as separate branches
            // (`KI == 1` and `T(KI,KI-1) == 0`), both setting `IP = 0`; merged
            // here because the two are indistinguishable.
            if ki == 1 || tget(t, n, ki, ki - 1) == 0.0 {
                ip = 0;
            } else {
                ip = -1;
            }

            // The KI-th eigenvalue (recomputed from T, as the reference does).
            let wr = tget(t, n, ki, ki);
            let wi = if ip != 0 {
                tget(t, n, ki, ki - 1).abs().sqrt() * tget(t, n, ki - 1, ki).abs().sqrt()
            } else {
                0.0
            };
            let smin = (ulp * (wr.abs() + wi.abs())).max(smlnum);

            if ip == 0 {
                // ================= real right eigenvector =================
                work1.fill(0.0);
                work1[ki - 1] = 1.0;
                for k in 1..ki {
                    work1[k - 1] = -tget(t, n, k, ki);
                }

                // Solve the upper quasi-triangular system
                //     (T(1:KI-1,1:KI-1) - WR) x = scale * work.
                let mut jnxt: isize = ki as isize - 1;
                let mut j: isize = ki as isize - 1;
                while j >= 1 {
                    if j <= jnxt {
                        let ju = j as usize;
                        let mut j1 = j;
                        let j2 = j;
                        jnxt = j - 1;
                        if j > 1 && tget(t, n, ju, ju - 1) != 0.0 {
                            j1 = j - 1;
                            jnxt = j - 2;
                        }
                        if j1 == j2 {
                            // 1x1 diagonal block.
                            let r = dlaln2_real(1, 1, smin, tget(t, n, ju, ju), 0.0, 0.0, 0.0, work1[ju - 1], 0.0, 0.0, 0.0, wr, 0.0);
                            let mut x11 = r.x11;
                            let mut scale = r.scale;
                            if r.xnorm > 1.0 && colnorm[ju - 1] > bignum / r.xnorm {
                                x11 /= r.xnorm;
                                scale /= r.xnorm;
                            }
                            if scale != 1.0 {
                                for l in 0..ki {
                                    work1[l] *= scale;
                                }
                            }
                            work1[ju - 1] = x11;
                            for l in 1..ju {
                                work1[l - 1] -= x11 * tget(t, n, l, ju);
                            }
                        } else {
                            // 2x2 diagonal block.
                            let r = dlaln2_real(
                                2,
                                1,
                                smin,
                                tget(t, n, ju - 1, ju - 1),
                                tget(t, n, ju, ju - 1),
                                tget(t, n, ju - 1, ju),
                                tget(t, n, ju, ju),
                                work1[ju - 2],
                                work1[ju - 1],
                                0.0,
                                0.0,
                                wr,
                                0.0,
                            );
                            let mut x11 = r.x11;
                            let mut x21 = r.x21;
                            let mut scale = r.scale;
                            if r.xnorm > 1.0 {
                                let beta = colnorm[ju - 2].max(colnorm[ju - 1]);
                                if beta > bignum / r.xnorm {
                                    x11 /= r.xnorm;
                                    x21 /= r.xnorm;
                                    scale /= r.xnorm;
                                }
                            }
                            if scale != 1.0 {
                                for l in 0..ki {
                                    work1[l] *= scale;
                                }
                            }
                            work1[ju - 2] = x11;
                            work1[ju - 1] = x21;
                            for l in 1..(ju - 1) {
                                work1[l - 1] -= x11 * tget(t, n, l, ju - 1);
                            }
                            for l in 1..(ju - 1) {
                                work1[l - 1] -= x21 * tget(t, n, l, ju);
                            }
                        }
                    }
                    j -= 1;
                }

                // Back-transform with the Schur vectors and normalize so that
                // the largest component has magnitude one:
                //     VR(:,KI) := Z * x.
                if ki > 1 {
                    for i in 1..=n {
                        let mut s = work1[ki - 1] * vget(vr, n, i, ki);
                        for l in 1..ki {
                            s += work1[l - 1] * vget(vr, n, i, l);
                        }
                        y1[i - 1] = s;
                    }
                    for i in 1..=n {
                        vset(vr, n, i, ki, y1[i - 1]);
                    }
                }
                let ii = crate::blas::idamax(&vr[(ki - 1) * n..ki * n]);
                let remax = 1.0 / vr[ii + (ki - 1) * n].abs();
                for i in 0..n {
                    vr[i + (ki - 1) * n] *= remax;
                }
            } else {
                // =================== complex right eigenvector ===================
                // Initial solve of the 2x2 block.
                if tget(t, n, ki - 1, ki).abs() >= tget(t, n, ki, ki - 1).abs() {
                    work1[ki - 2] = 1.0;
                    work2[ki - 1] = wi / tget(t, n, ki - 1, ki);
                } else {
                    work1[ki - 2] = -wi / tget(t, n, ki, ki - 1);
                    work2[ki - 1] = 1.0;
                }
                work1[ki - 1] = 0.0;
                work2[ki - 2] = 0.0;

                // Form the right-hand side.
                for k in 1..=(ki - 2) {
                    work1[k - 1] = -work1[ki - 2] * tget(t, n, k, ki - 1);
                    work2[k - 1] = -work2[ki - 1] * tget(t, n, k, ki);
                }

                // Solve (T(1:KI-2,1:KI-2) - (WR + i WI)) x = scale*(work1 + i work2).
                let mut jnxt: isize = ki as isize - 2;
                let mut j: isize = ki as isize - 2;
                while j >= 1 {
                    if j <= jnxt {
                        let ju = j as usize;
                        let mut j1 = j;
                        let j2 = j;
                        jnxt = j - 1;
                        if j > 1 && tget(t, n, ju, ju - 1) != 0.0 {
                            j1 = j - 1;
                            jnxt = j - 2;
                        }
                        if j1 == j2 {
                            // 1x1 diagonal block, two right-hand sides.
                            let r = dlaln2_real(
                                1,
                                2,
                                smin,
                                tget(t, n, ju, ju),
                                0.0,
                                0.0,
                                0.0,
                                work1[ju - 1],
                                0.0,
                                work2[ju - 1],
                                0.0,
                                wr,
                                wi,
                            );
                            let mut x11 = r.x11;
                            let mut x12 = r.x12;
                            let mut scale = r.scale;
                            if r.xnorm > 1.0 && colnorm[ju - 1] > bignum / r.xnorm {
                                x11 /= r.xnorm;
                                x12 /= r.xnorm;
                                scale /= r.xnorm;
                            }
                            if scale != 1.0 {
                                for l in 0..ki {
                                    work1[l] *= scale;
                                    work2[l] *= scale;
                                }
                            }
                            work1[ju - 1] = x11;
                            work2[ju - 1] = x12;
                            for l in 1..ju {
                                work1[l - 1] -= x11 * tget(t, n, l, ju);
                            }
                            for l in 1..ju {
                                work2[l - 1] -= x12 * tget(t, n, l, ju);
                            }
                        } else {
                            // 2x2 diagonal block, complex system.
                            let r = dlaln2_real(
                                2,
                                2,
                                smin,
                                tget(t, n, ju - 1, ju - 1),
                                tget(t, n, ju, ju - 1),
                                tget(t, n, ju - 1, ju),
                                tget(t, n, ju, ju),
                                work1[ju - 2],
                                work1[ju - 1],
                                work2[ju - 2],
                                work2[ju - 1],
                                wr,
                                wi,
                            );
                            let mut x11 = r.x11;
                            let mut x21 = r.x21;
                            let mut x12 = r.x12;
                            let mut x22 = r.x22;
                            let mut scale = r.scale;
                            if r.xnorm > 1.0 {
                                let beta = colnorm[ju - 2].max(colnorm[ju - 1]);
                                if beta > bignum / r.xnorm {
                                    let rec = 1.0 / r.xnorm;
                                    x11 *= rec;
                                    x12 *= rec;
                                    x21 *= rec;
                                    x22 *= rec;
                                    scale *= rec;
                                }
                            }
                            if scale != 1.0 {
                                for l in 0..ki {
                                    work1[l] *= scale;
                                    work2[l] *= scale;
                                }
                            }
                            work1[ju - 2] = x11;
                            work1[ju - 1] = x21;
                            work2[ju - 2] = x12;
                            work2[ju - 1] = x22;
                            for l in 1..(ju - 1) {
                                work1[l - 1] -= x11 * tget(t, n, l, ju - 1);
                            }
                            for l in 1..(ju - 1) {
                                work1[l - 1] -= x21 * tget(t, n, l, ju);
                            }
                            for l in 1..(ju - 1) {
                                work2[l - 1] -= x12 * tget(t, n, l, ju - 1);
                            }
                            for l in 1..(ju - 1) {
                                work2[l - 1] -= x22 * tget(t, n, l, ju);
                            }
                        }
                    }
                    j -= 1;
                }

                // Back-transform both columns and normalize so that the largest
                // (|re| + |im|) is one:
                //     VR(:,KI-1) := Z * work1,  VR(:,KI) := Z * work2.
                for i in 1..=n {
                    let mut s1 = work1[ki - 2] * vget(vr, n, i, ki - 1);
                    let mut s2 = work2[ki - 1] * vget(vr, n, i, ki);
                    for l in 1..(ki - 1) {
                        s1 += work1[l - 1] * vget(vr, n, i, l);
                        s2 += work2[l - 1] * vget(vr, n, i, l);
                    }
                    y1[i - 1] = s1;
                    y2[i - 1] = s2;
                }
                for i in 0..n {
                    vr[i + (ki - 2) * n] = y1[i];
                    vr[i + (ki - 1) * n] = y2[i];
                }
                let mut emax = 0.0f64;
                for i in 0..n {
                    emax = emax.max(vr[i + (ki - 2) * n].abs() + vr[i + (ki - 1) * n].abs());
                }
                let remax = 1.0 / emax;
                for i in 0..n {
                    vr[i + (ki - 2) * n] *= remax;
                    vr[i + (ki - 1) * n] *= remax;
                }
            }
        }

        // 130 CONTINUE
        if ip == 1 {
            ip = 0;
        } else if ip == -1 {
            ip = 1;
        }
        ki -= 1;
    }
}

/// Euclidean norm of column `j0` (0-based) of a column-major matrix.
///
/// A column of a column-major matrix is contiguous, so this is the crate's
/// [`dnrm2`](crate::blas::dnrm2) on a slice.
fn column_norm2(v: &[f64], n: usize, j0: usize) -> f64 {
    crate::blas::dnrm2(&v[j0 * n..(j0 + 1) * n])
}

/// `dgebal` with `JOB = 'B'` (LAPACK 3.12.1, `SRC/dgebal.f`): permute and scale a
/// general real matrix to improve the conditioning of its eigenvalues.
///
/// The permutation pushes rows (then columns) that isolate an eigenvalue to the
/// bottom (top), so that the matrix takes the block form
///
/// ```text
///       ( T1  X  Y )
/// P A P ( 0   B  Z )
///       ( 0   0  T2 )
/// ```
///
/// with `T1` and `T2` upper triangular and `B` the active block `ilo..=ihi`
/// (0-based). Scaling then makes the 1-norms of each row of `B` and its
/// corresponding column nearly equal by a diagonal similarity.
///
/// The returned `scale` vector follows LAPACK: entries inside `ilo..=ihi` hold
/// the scaling factors, entries outside hold the 1-based source index of the
/// permutation (used by `dgebak_right` to undo it).
fn dgebal(a: &mut [f64], n: usize) -> (usize, usize, Vec<f64>) {
    // PARAMETER ( SCLFAC = 2.0 ) and PARAMETER ( FACTOR = 0.95D+0 ).
    const SCLFAC: f64 = 2.0;
    const FACTOR: f64 = 0.95;

    let mut scale = vec![1.0f64; n];
    if n == 0 {
        return (0, 0, scale);
    }

    // K and L are LAPACK's ILO and IHI, 0-based here.
    let mut k = 0usize;
    let mut l = n - 1;

    // ==== Permutation to isolate eigenvalues if possible. ====
    //
    // Search for rows isolating an eigenvalue and push them down.
    let mut noconv = true;
    while noconv {
        noconv = false;
        // Fortran fixes the DO bound when the loop is entered; L is decremented
        // inside the body, so capture the bound first.
        let mut i = l as isize;
        while i >= 0 {
            let iu = i as usize;
            let mut canswap = true;
            for j in 0..=l {
                if iu != j && a[iu + j * n] != 0.0 {
                    canswap = false;
                    break;
                }
            }
            if canswap {
                scale[l] = (iu + 1) as f64;
                if iu != l {
                    // Swap rows (DSWAP(L, A(1,I), 1, A(1,L), 1)) and columns
                    // (DSWAP(N-K+1, A(I,K), LDA, A(L,K), LDA)).
                    for r in 0..=l {
                        a.swap(r + iu * n, r + l * n);
                    }
                    for c in k..n {
                        a.swap(iu + c * n, l + c * n);
                    }
                }
                noconv = true;
                if l == 0 {
                    return (0, 0, scale);
                }
                l -= 1;
            }
            i -= 1;
        }
    }

    // Search for columns isolating an eigenvalue and push them left.
    noconv = true;
    while noconv {
        noconv = false;
        let mut j = k;
        while j <= l {
            let mut canswap = true;
            for i in k..=l {
                if i != j && a[i + j * n] != 0.0 {
                    canswap = false;
                    break;
                }
            }
            if canswap {
                scale[k] = (j + 1) as f64;
                if j != k {
                    for r in 0..=l {
                        a.swap(r + j * n, r + k * n);
                    }
                    for c in k..n {
                        a.swap(j + c * n, k + c * n);
                    }
                }
                noconv = true;
                k += 1;
            }
            j += 1;
        }
    }

    // Initialize SCALE for the non-permuted submatrix.
    for i in k..=l {
        scale[i] = 1.0;
    }

    // ==== Balance the submatrix in rows K to L. ====
    //
    // SFMIN1 = DLAMCH('S') / DLAMCH('P').
    let sfmin1 = f64::MIN_POSITIVE / f64::EPSILON;
    let sfmax1 = 1.0 / sfmin1;
    let sfmin2 = sfmin1 * SCLFAC;
    let sfmax2 = 1.0 / sfmin2;

    let mut noconv = true;
    while noconv {
        noconv = false;
        for i in k..=l {
            let col: Vec<f64> = (k..=l).map(|r| a[r + i * n]).collect();
            let row: Vec<f64> = (k..=l).map(|c| a[i + c * n]).collect();
            let mut c = crate::blas::dnrm2(&col);
            let mut r = crate::blas::dnrm2(&row);

            // ICA = IDAMAX(L, A(1,I), 1); IRA = IDAMAX(N-K+1, A(I,K), LDA).
            let col_full: Vec<f64> = (0..=l).map(|rr| a[rr + i * n]).collect();
            let ica = crate::blas::idamax(&col_full);
            let mut ca = a[ica + i * n].abs();
            let row_full: Vec<f64> = (k..n).map(|cc| a[i + cc * n]).collect();
            let ira = crate::blas::idamax(&row_full) + k;
            let mut ra = a[i + ira * n].abs();

            // Guard against zero C or R due to underflow.
            if c == 0.0 || r == 0.0 {
                continue;
            }
            if (c + ca + r + ra).is_nan() {
                // DISNAN(C+CA+R+RA): the reference reports an error here. The
                // inputs to dgeev are finite, so this is unreachable in practice;
                // leaving the row unbalanced is the graceful fallback.
                continue;
            }

            let mut g = r / SCLFAC;
            let mut f = 1.0f64;
            let s = c + r;

            while c < g && f.max(c).max(ca) < sfmax2 && r.min(g).min(ra) > sfmin2 {
                f *= SCLFAC;
                c *= SCLFAC;
                ca *= SCLFAC;
                r /= SCLFAC;
                g /= SCLFAC;
                ra /= SCLFAC;
            }

            g = c / SCLFAC;

            while g >= r && r.max(ra) < sfmax2 && f.min(c).min(g).min(ca) > sfmin2 {
                f /= SCLFAC;
                c /= SCLFAC;
                g /= SCLFAC;
                ca /= SCLFAC;
                r *= SCLFAC;
                ra *= SCLFAC;
            }

            // Now balance.
            if c + r >= FACTOR * s {
                continue;
            }
            if f < 1.0 && scale[i] < 1.0 {
                if f * scale[i] <= sfmin1 {
                    continue;
                }
            }
            if f > 1.0 && scale[i] > 1.0 {
                if scale[i] >= sfmax1 / f {
                    continue;
                }
            }
            let g2 = 1.0 / f;
            scale[i] *= f;
            noconv = true;

            // DSCAL(N-K+1, G, A(I,K), LDA) then DSCAL(L, F, A(1,I), 1): scale
            // row I across columns K..N and column I across rows 1..L.
            for cc in k..n {
                a[i + cc * n] *= g2;
            }
            for rr in 0..=l {
                a[rr + i * n] *= f;
            }
        }
    }

    (k, l, scale)
}

/// `dgebak` with `JOB = 'B'`, `SIDE = 'R'` (LAPACK 3.12.1, `SRC/dgebak.f`):
/// undo the balancing that `dgebal` recorded in `scale`.
///
/// The right eigenvectors of the balanced matrix are those of the original
/// matrix with row `i` multiplied by `scale(i)` (for `i` in `ilo..=ihi`) and the
/// isolated rows permuted back.
fn dgebak_right(v: &mut [f64], n: usize, ilo: usize, ihi: usize, scale: &[f64]) {
    if n == 0 {
        return;
    }

    // Backward balance.
    if ilo != ihi {
        for i in ilo..=ihi {
            let s = scale[i];
            for j in 0..n {
                v[i + j * n] *= s;
            }
        }
    }

    // Backward permutation: I = ILO-1 step -1 until 1, then IHI+1 step 1 until N.
    for ii in 0..n {
        let i = if ii < ilo { ilo - ii - 1 } else { ii };
        if i >= ilo && i <= ihi {
            continue;
        }
        let k = scale[i] as usize - 1;
        if k == i {
            continue;
        }
        for j in 0..n {
            v.swap(i + j * n, k + j * n);
        }
    }
}

/// `dgeev` (right eigenvectors only): the eigenvalues and right eigenvectors of
/// a general real square matrix.
///
/// `a` is `n x n` column-major and is destroyed. On exit:
///
/// * the return value is `(wr, wi)`, the eigenvalues, in the same order as the
///   eigenvector columns; a complex conjugate pair is adjacent with `wi > 0`
///   first;
/// * `v_real` and `v_imag` hold the right eigenvectors as their real and
///   imaginary parts: `v(j) = v_real[:, j] + i * v_imag[:, j]`. A real
///   eigenvalue's column has zero imaginary part, and for a pair at `(j, j+1)`
///   the two columns are conjugates, `v(j+1) = conj(v(j))`. This is LAPACK's
///   compact `VR` layout.
///
/// # What this follows
///
/// `dgeev.f` with `JOBVL = 'N'`, `JOBVR = 'V'`, and the `dgeev.f` dispatch to
/// `dtrevc3` (the unblocked `NB = 1` branch, which is what `dgeev`'s
/// `lwork = 4n + 1` selects):
///
/// 1. scale `A` if its largest element is outside `[SMLNUM, BIGNUM]`;
/// 2. `dgebal` (`JOB = 'B'`): permute and scale to improve conditioning;
/// 3. [`dgehrd`]: `A = Q H Q^T`, with `Q` accumulated;
/// 4. `dhseqr_with_schur_vectors` (`COMPZ = 'V'`): `A = Z T Z^T`, accumulating
///    the Schur vectors into `Q`, so that afterwards `Z` holds them;
/// 5. `dtrevc_right` (`DTREVC3` with `SIDE = 'R'`, `HOWMNY = 'B'`):
///    eigenvectors of `T`, back-transformed by `Z`;
/// 6. `dgebak_right` (`JOB = 'B'`, `SIDE = 'R'`): undo the balancing of the
///    eigenvectors;
/// 7. normalize the eigenvectors, then unscale the eigenvalues.
///
/// Balancing is what makes the eigenvalues of a defectively repeated eigenvalue
/// come out exact: the permutation in step 2 isolates such an eigenvalue, so the
/// QR iteration never has to separate it from its neighbors.
///
/// # Errors
///
/// [`LapackError::Convergence`] if the QR iteration in step 4 fails to converge.
pub fn dgeev(v_real: &mut [f64], v_imag: &mut [f64], a: &mut [f64], n: usize) -> Result<(Vec<f64>, Vec<f64>), LapackError> {
    assert_eq!(a.len(), n * n, "dgeev: a must be n*n");
    assert_eq!(v_real.len(), n * n, "dgeev: v_real must be n*n");
    assert_eq!(v_imag.len(), n * n, "dgeev: v_imag must be n*n");

    let mut wr = vec![0.0; n];
    let mut wi = vec![0.0; n];
    if n == 0 {
        return Ok((wr, wi));
    }

    // ---- machine constants (dgeev.f) ----
    let eps = f64::EPSILON;
    let smlnum0 = f64::MIN_POSITIVE;
    let smlnum = smlnum0.sqrt() / eps;
    let bignum = 1.0 / smlnum;

    // ---- scale A if its largest element is outside [SMLNUM, BIGNUM] ----
    let mut anrm = 0.0f64;
    for x in a.iter() {
        anrm = anrm.max(x.abs());
    }
    let mut cscale = 0.0f64;
    let scalea = if anrm > 0.0 && anrm < smlnum {
        cscale = smlnum;
        true
    } else if anrm > bignum {
        cscale = bignum;
        true
    } else {
        false
    };
    if scalea {
        let f = cscale / anrm;
        for x in a.iter_mut() {
            *x *= f;
        }
    }

    // ---- DGEBAL('B'): permute and scale to improve conditioning ----
    let (ilo, ihi, bal_scale) = dgebal(a, n);

    // ---- DGEHRD: A = Q H Q^T ----
    let mut q = vec![0.0; n * n];
    dgehrd(a, n, Some(&mut q));

    // ---- DHSEQR with COMPZ = 'V': A = Z T Z^T (Z accumulated into q) ----
    dhseqr_with_schur_vectors(a, &mut wr, &mut wi, &mut q, n)?;

    // ---- DTREVC3('R','B'): eigenvectors of T, back-transformed by Z ----
    let mut vr = q;
    dtrevc_right(a, &mut vr, n);

    // ---- DGEBAK('B','R'): undo the balancing of the eigenvectors ----
    dgebak_right(&mut vr, n, ilo, ihi, &bal_scale);

    // ---- normalize: unit Euclidean norm, largest component real ----
    // (dgeev.f, the WANTVR branch; `i0` is LAPACK's 1-based I minus one)
    let mut i0 = 0usize;
    while i0 < n {
        if wi[i0] == 0.0 {
            let scl = 1.0 / column_norm2(&vr, n, i0);
            for r in 0..n {
                vr[r + i0 * n] *= scl;
            }
            i0 += 1;
        } else if wi[i0] > 0.0 {
            let scl = 1.0 / f64::hypot(column_norm2(&vr, n, i0), column_norm2(&vr, n, i0 + 1));
            for r in 0..n {
                vr[r + i0 * n] *= scl;
                vr[r + (i0 + 1) * n] *= scl;
            }
            // Rotate the largest component onto the real axis.
            let mut k = 0usize;
            let mut best = -1.0f64;
            for r in 0..n {
                let mag = vr[r + i0 * n] * vr[r + i0 * n] + vr[r + (i0 + 1) * n] * vr[r + (i0 + 1) * n];
                if mag > best {
                    best = mag;
                    k = r;
                }
            }
            let (cs, sn, _r) = dlartg(vr[k + i0 * n], vr[k + (i0 + 1) * n]);
            for r in 0..n {
                let x = vr[r + i0 * n];
                let y = vr[r + (i0 + 1) * n];
                vr[r + i0 * n] = cs * x + sn * y;
                vr[r + (i0 + 1) * n] = cs * y - sn * x;
            }
            vr[k + (i0 + 1) * n] = 0.0;
            i0 += 2;
        } else {
            // The second member of a pair: already written by the first.
            i0 += 1;
        }
    }

    // ---- undo the scaling of A on the eigenvalues ----
    if scalea {
        let f = anrm / cscale;
        for x in wr.iter_mut() {
            *x *= f;
        }
        for x in wi.iter_mut() {
            *x *= f;
        }
    }

    // ---- unpack the compact storage into real and imaginary parts ----
    for i in 0..n {
        if wi[i] == 0.0 {
            for r in 0..n {
                v_real[r + i * n] = vr[r + i * n];
                v_imag[r + i * n] = 0.0;
            }
        } else if wi[i] > 0.0 {
            for r in 0..n {
                let re = vr[r + i * n];
                let im = vr[r + (i + 1) * n];
                v_real[r + i * n] = re;
                v_imag[r + i * n] = im;
                v_real[r + (i + 1) * n] = re;
                v_imag[r + (i + 1) * n] = -im;
            }
        }
    }

    Ok((wr, wi))
}

/// `dtrevc3` (`dtrevc3.f`, LAPACK 3.12.1) with `SIDE = 'L'`, `HOWMNY = 'B'` and
/// the unblocked `NB = 1` back-transformation.
///
/// On entry `t` is the real Schur form (quasi-upper-triangular) and `vl` holds
/// the orthogonal Schur vectors `Z`; on exit `vl` holds the left eigenvectors in
/// LAPACK's compact storage: for a conjugate pair at `(j, j+1)`, columns `j` and
/// `j + 1` are the real and imaginary parts of `u(j)`, with `u(j+1) = conj(u(j))`.
/// The defining equation is `u(j)^T A = lambda(j) u(j)^T`.
///
/// The only `dlaln2` call with `LTRANS = .TRUE.` is the 2x2 block, and `dlaln2`
/// uses `LTRANS` for nothing but the off-diagonal swap in `CR`; passing the two
/// off-diagonals in the opposite argument slots is therefore exactly the
/// transposed system, and no second `dlaln2` is needed.
/// `VCRIT = BIGNUM` is set again after each rescale, and the reference does the
/// same even though the next statement overwrites it (`VCRIT = BIGNUM / VMAX`):
/// the write is dead there too. Kept so a diff against `dtrevc3.f` still lines
/// up; `unused_assignments` is the only lint it trips.
#[allow(unused_assignments)]
fn dtrevc_left(t: &[f64], vl: &mut [f64], n: usize) {
    let unfl = f64::MIN_POSITIVE;
    let ulp = f64::EPSILON;
    let smlnum = unfl * (n as f64 / ulp);
    let bignum = (1.0 - ulp) / smlnum;

    // 1-norm of the strictly upper triangular part of each column of T.
    let mut colnorm = vec![0.0; n];
    for j in 2..=n {
        let mut s = 0.0;
        for i in 1..j {
            s += tget(t, n, i, j).abs();
        }
        colnorm[j - 1] = s;
    }

    // The two columns of `WORK` the reference keeps for the current eigenvector:
    // `w1` is the real part and `w2` the imaginary part.
    let mut w1 = vec![0.0; n];
    let mut w2 = vec![0.0; n];
    let mut y1 = vec![0.0; n];
    let mut y2 = vec![0.0; n];

    // IP = 0: real eigenvalue; 1: the pair starting here was just processed;
    // -1: the second member of a pair, to be skipped.
    let mut ip = 0i32;
    let mut ki = 1usize;
    while ki <= n {
        if ip == 1 {
            ip = -1;
            ki += 1;
            continue;
        }
        if ki == n || tget(t, n, ki + 1, ki) == 0.0 {
            ip = 0;
        } else {
            ip = 1;
        }
        let wr = tget(t, n, ki, ki);
        let wi = if ip != 0 {
            tget(t, n, ki, ki + 1).abs().sqrt() * tget(t, n, ki + 1, ki).abs().sqrt()
        } else {
            0.0
        };
        let smin = (ulp * (wr.abs() + wi.abs())).max(smlnum);

        if ip == 0 {
            // ================== real left eigenvector ==================
            w1.fill(0.0);
            w1[ki - 1] = 1.0;
            for k in (ki + 1)..=n {
                w1[k - 1] = -tget(t, n, ki, k);
            }

            // Solve (T(KI+1:N,KI+1:N) - WR)^T x = scale * work by forward
            // substitution over the quasi-triangular blocks.
            let mut vmax = 1.0;
            let mut vcrit = bignum;
            let mut jnxt = ki + 1;
            let mut j = ki + 1;
            while j <= n {
                if j >= jnxt {
                    let j1 = j;
                    let mut j2 = j;
                    jnxt = j + 1;
                    if j < n && tget(t, n, j + 1, j) != 0.0 {
                        j2 = j + 1;
                        jnxt = j + 2;
                    }
                    if j1 == j2 {
                        if colnorm[j - 1] > vcrit {
                            let rec = 1.0 / vmax;
                            for l in (ki - 1)..n {
                                w1[l] *= rec;
                            }
                            vmax = 1.0;
                            vcrit = bignum;
                        }
                        let mut s = 0.0;
                        for l in (ki + 1)..j {
                            s += tget(t, n, l, j) * w1[l - 1];
                        }
                        w1[j - 1] -= s;
                        let r = dlaln2_real(1, 1, smin, tget(t, n, j, j), 0.0, 0.0, 0.0, w1[j - 1], 0.0, 0.0, 0.0, wr, 0.0);
                        if r.scale != 1.0 {
                            for l in (ki - 1)..n {
                                w1[l] *= r.scale;
                            }
                        }
                        w1[j - 1] = r.x11;
                        vmax = w1[j - 1].abs().max(vmax);
                        vcrit = bignum / vmax;
                    } else {
                        let beta = colnorm[j - 1].max(colnorm[j]);
                        if beta > vcrit {
                            let rec = 1.0 / vmax;
                            for l in (ki - 1)..n {
                                w1[l] *= rec;
                            }
                            vmax = 1.0;
                            vcrit = bignum;
                        }
                        let mut s = 0.0;
                        for l in (ki + 1)..j {
                            s += tget(t, n, l, j) * w1[l - 1];
                        }
                        w1[j - 1] -= s;
                        let mut s2 = 0.0;
                        for l in (ki + 1)..j {
                            s2 += tget(t, n, l, j + 1) * w1[l - 1];
                        }
                        w1[j] -= s2;
                        // LTRANS = .TRUE.: the off-diagonals are swapped.
                        let r = dlaln2_real(
                            2,
                            1,
                            smin,
                            tget(t, n, j, j),
                            tget(t, n, j, j + 1),
                            tget(t, n, j + 1, j),
                            tget(t, n, j + 1, j + 1),
                            w1[j - 1],
                            w1[j],
                            0.0,
                            0.0,
                            wr,
                            0.0,
                        );
                        if r.scale != 1.0 {
                            for l in (ki - 1)..n {
                                w1[l] *= r.scale;
                            }
                        }
                        w1[j - 1] = r.x11;
                        w1[j] = r.x21;
                        vmax = w1[j - 1].abs().max(w1[j].abs()).max(vmax);
                        vcrit = bignum / vmax;
                    }
                }
                j += 1;
            }

            // Back-transform with the Schur vectors: VL(:,KI) := Z * x.
            if ki < n {
                for i in 1..=n {
                    let mut s = w1[ki - 1] * vget(vl, n, i, ki);
                    for l in (ki + 1)..=n {
                        s += w1[l - 1] * vget(vl, n, i, l);
                    }
                    y1[i - 1] = s;
                }
                for i in 1..=n {
                    vset(vl, n, i, ki, y1[i - 1]);
                }
            }
            let ii = crate::blas::idamax(&vl[(ki - 1) * n..ki * n]);
            let remax = 1.0 / vl[ii + (ki - 1) * n].abs();
            for i in 0..n {
                vl[i + (ki - 1) * n] *= remax;
            }
        } else {
            // =============== complex left eigenvector pair ===============
            // The first column of the pair is the real part, the second the
            // imaginary part, with the "seam" component of each set to zero.
            if tget(t, n, ki, ki + 1).abs() >= tget(t, n, ki + 1, ki).abs() {
                w1[ki - 1] = wi / tget(t, n, ki, ki + 1);
                w2[ki] = 1.0;
            } else {
                w1[ki - 1] = 1.0;
                w2[ki] = -wi / tget(t, n, ki + 1, ki);
            }
            w1[ki] = 0.0;
            w2[ki - 1] = 0.0;
            for k in (ki + 2)..=n {
                w1[k - 1] = -w1[ki - 1] * tget(t, n, ki, k);
                w2[k - 1] = -w2[ki] * tget(t, n, ki + 1, k);
            }

            let mut vmax = 1.0;
            let mut vcrit = bignum;
            let mut jnxt = ki + 2;
            let mut j = ki + 2;
            while j <= n {
                if j >= jnxt {
                    let j1 = j;
                    let mut j2 = j;
                    jnxt = j + 1;
                    if j < n && tget(t, n, j + 1, j) != 0.0 {
                        j2 = j + 1;
                        jnxt = j + 2;
                    }
                    if j1 == j2 {
                        if colnorm[j - 1] > vcrit {
                            let rec = 1.0 / vmax;
                            for l in (ki - 1)..n {
                                w1[l] *= rec;
                                w2[l] *= rec;
                            }
                            vmax = 1.0;
                            vcrit = bignum;
                        }
                        let mut s1 = 0.0;
                        for l in (ki + 2)..j {
                            s1 += tget(t, n, l, j) * w1[l - 1];
                        }
                        w1[j - 1] -= s1;
                        let mut s2 = 0.0;
                        for l in (ki + 2)..j {
                            s2 += tget(t, n, l, j) * w2[l - 1];
                        }
                        w2[j - 1] -= s2;
                        let r = dlaln2_real(1, 2, smin, tget(t, n, j, j), 0.0, 0.0, 0.0, w1[j - 1], 0.0, w2[j - 1], 0.0, wr, -wi);
                        if r.scale != 1.0 {
                            for l in (ki - 1)..n {
                                w1[l] *= r.scale;
                                w2[l] *= r.scale;
                            }
                        }
                        w1[j - 1] = r.x11;
                        w2[j - 1] = r.x12;
                        vmax = w1[j - 1].abs().max(w2[j - 1].abs()).max(vmax);
                        vcrit = bignum / vmax;
                    } else {
                        let beta = colnorm[j - 1].max(colnorm[j]);
                        if beta > vcrit {
                            let rec = 1.0 / vmax;
                            for l in (ki - 1)..n {
                                w1[l] *= rec;
                                w2[l] *= rec;
                            }
                            vmax = 1.0;
                            vcrit = bignum;
                        }
                        let mut s1 = 0.0;
                        for l in (ki + 2)..j {
                            s1 += tget(t, n, l, j) * w1[l - 1];
                        }
                        w1[j - 1] -= s1;
                        let mut s2 = 0.0;
                        for l in (ki + 2)..j {
                            s2 += tget(t, n, l, j) * w2[l - 1];
                        }
                        w2[j - 1] -= s2;
                        let mut s3 = 0.0;
                        for l in (ki + 2)..j {
                            s3 += tget(t, n, l, j + 1) * w1[l - 1];
                        }
                        w1[j] -= s3;
                        let mut s4 = 0.0;
                        for l in (ki + 2)..j {
                            s4 += tget(t, n, l, j + 1) * w2[l - 1];
                        }
                        w2[j] -= s4;
                        // LTRANS = .TRUE.: the off-diagonals are swapped.
                        let r = dlaln2_real(
                            2,
                            2,
                            smin,
                            tget(t, n, j, j),
                            tget(t, n, j, j + 1),
                            tget(t, n, j + 1, j),
                            tget(t, n, j + 1, j + 1),
                            w1[j - 1],
                            w1[j],
                            w2[j - 1],
                            w2[j],
                            wr,
                            -wi,
                        );
                        if r.scale != 1.0 {
                            for l in (ki - 1)..n {
                                w1[l] *= r.scale;
                                w2[l] *= r.scale;
                            }
                        }
                        w1[j - 1] = r.x11;
                        w2[j - 1] = r.x12;
                        w1[j] = r.x21;
                        w2[j] = r.x22;
                        vmax = w1[j - 1].abs().max(w2[j - 1].abs()).max(w1[j].abs()).max(w2[j].abs()).max(vmax);
                        vcrit = bignum / vmax;
                    }
                }
                j += 1;
            }

            // Back-transform: the two columns are Q*w1 and Q*w2.
            if ki < n - 1 {
                for i in 1..=n {
                    let mut s1 = w1[ki - 1] * vget(vl, n, i, ki);
                    let mut s2 = w2[ki] * vget(vl, n, i, ki + 1);
                    for l in (ki + 2)..=n {
                        s1 += w1[l - 1] * vget(vl, n, i, l);
                        s2 += w2[l - 1] * vget(vl, n, i, l);
                    }
                    y1[i - 1] = s1;
                    y2[i - 1] = s2;
                }
                for i in 0..n {
                    vl[i + (ki - 1) * n] = y1[i];
                    vl[i + ki * n] = y2[i];
                }
            } else {
                for i in 0..n {
                    vl[i + (ki - 1) * n] *= w1[ki - 1];
                    vl[i + ki * n] *= w2[ki];
                }
            }
            let mut emax = 0.0f64;
            for i in 0..n {
                emax = emax.max(vl[i + (ki - 1) * n].abs() + vl[i + ki * n].abs());
            }
            let remax = 1.0 / emax;
            for i in 0..n {
                vl[i + (ki - 1) * n] *= remax;
                vl[i + ki * n] *= remax;
            }
        }
        ki += 1;
    }
}

/// `dgebak` with `JOB = 'B'`, `SIDE = 'L'`: undo the balancing that `dgebal`
/// recorded in `scale`.
///
/// The left eigenvectors of the balanced matrix are those of the original matrix
/// with row `i` divided by `scale(i)` (for `i` in `ilo..=ihi`), and the isolated
/// rows permuted back.
fn dgebak_left(v: &mut [f64], n: usize, ilo: usize, ihi: usize, scale: &[f64]) {
    if n == 0 {
        return;
    }

    // Backward balance: the division is the only difference from SIDE='R'.
    if ilo != ihi {
        for i in ilo..=ihi {
            let s = 1.0 / scale[i];
            for j in 0..n {
                v[i + j * n] *= s;
            }
        }
    }

    // Backward permutation, identical to SIDE='R'.
    for ii in 0..n {
        let i = if ii < ilo { ilo - ii - 1 } else { ii };
        if i >= ilo && i <= ihi {
            continue;
        }
        let k = scale[i] as usize - 1;
        if k == i {
            continue;
        }
        for j in 0..n {
            v.swap(i + j * n, k + j * n);
        }
    }
}

/// Normalize a set of eigenvector columns the way `dgeev` does: unit Euclidean
/// norm, and for a conjugate pair the largest component is rotated onto the real
/// axis with [`dlartg`].
fn dgeev_normalize_columns(v: &mut [f64], n: usize, wi: &[f64]) {
    let mut k = 0usize;
    while k < n {
        if wi[k] == 0.0 {
            let scl = 1.0 / column_norm2(v, n, k);
            for r in 0..n {
                v[r + k * n] *= scl;
            }
            k += 1;
        } else if wi[k] > 0.0 {
            let scl = 1.0 / f64::hypot(column_norm2(v, n, k), column_norm2(v, n, k + 1));
            for r in 0..n {
                v[r + k * n] *= scl;
                v[r + (k + 1) * n] *= scl;
            }
            let mut kk = 0usize;
            let mut best = -1.0f64;
            for r in 0..n {
                let mag = v[r + k * n] * v[r + k * n] + v[r + (k + 1) * n] * v[r + (k + 1) * n];
                if mag > best {
                    best = mag;
                    kk = r;
                }
            }
            let (cs, sn, _r) = dlartg(v[kk + k * n], v[kk + (k + 1) * n]);
            for r in 0..n {
                let x = v[r + k * n];
                let y = v[r + (k + 1) * n];
                v[r + k * n] = cs * x + sn * y;
                v[r + (k + 1) * n] = cs * y - sn * x;
            }
            v[kk + (k + 1) * n] = 0.0;
            k += 2;
        } else {
            k += 1;
        }
    }
}

/// Unpack LAPACK's compact eigenvector storage into real and imaginary parts.
///
/// For a real eigenvalue the column is real; for a pair at `(k, k + 1)`, `u(k)`
/// has the two columns as its real and imaginary parts and `u(k + 1)` is its
/// conjugate, matching `dgeev`'s documented layout for both `VL` and `VR`.
fn dgeev_unpack(raw: &[f64], wi: &[f64], re: &mut [f64], im: &mut [f64], n: usize) {
    for i in 0..n {
        if wi[i] == 0.0 {
            for r in 0..n {
                re[r + i * n] = raw[r + i * n];
                im[r + i * n] = 0.0;
            }
        } else if wi[i] > 0.0 {
            for r in 0..n {
                let a = raw[r + i * n];
                let b = raw[r + (i + 1) * n];
                re[r + i * n] = a;
                im[r + i * n] = b;
                re[r + (i + 1) * n] = a;
                im[r + (i + 1) * n] = -b;
            }
        }
    }
}

/// `dgeev` with `JOBVL = 'V'` and `JOBVR = 'V'`: the eigenvalues and the left and
/// right eigenvectors of a general real square matrix.
///
/// This is the same pipeline as [`dgeev`] -- scale, `dgebal`, [`dgehrd`],
/// `dhseqr_with_schur_vectors`, the `dtrevc3` back-substitutions,
/// `dgebak_right`/`dgebak_left`, normalization -- with the left eigenvectors
/// taken from the *same* Schur factorization by `dtrevc_left` rather than from
/// a second decomposition of `A^T`.
///
/// `a` is `n x n` column-major and is destroyed. On exit the return value is
/// `(wr, wi)`; `u_real`/`u_imag` and `v_real`/`v_imag` hold the left and right
/// eigenvectors, with `u(j) = u_real[:, j] + i u_imag[:, j]` (and likewise for
/// `v`), and the second member of a pair the conjugate of the first.
///
/// # Errors
///
/// [`LapackError::Convergence`] if the QR iteration fails to converge.
#[allow(clippy::too_many_arguments)]
pub fn dgeev_lr(
    u_real: &mut [f64],
    u_imag: &mut [f64],
    v_real: &mut [f64],
    v_imag: &mut [f64],
    a: &mut [f64],
    n: usize,
) -> Result<(Vec<f64>, Vec<f64>), LapackError> {
    assert_eq!(a.len(), n * n, "dgeev_lr: a must be n*n");
    assert_eq!(u_real.len(), n * n, "dgeev_lr: u_real must be n*n");
    assert_eq!(u_imag.len(), n * n, "dgeev_lr: u_imag must be n*n");
    assert_eq!(v_real.len(), n * n, "dgeev_lr: v_real must be n*n");
    assert_eq!(v_imag.len(), n * n, "dgeev_lr: v_imag must be n*n");

    let mut wr = vec![0.0; n];
    let mut wi = vec![0.0; n];
    if n == 0 {
        return Ok((wr, wi));
    }

    // ---- machine constants (dgeev.f) ----
    let eps = f64::EPSILON;
    let smlnum0 = f64::MIN_POSITIVE;
    let smlnum = smlnum0.sqrt() / eps;
    let bignum = 1.0 / smlnum;

    // ---- scale A if its largest element is outside [SMLNUM, BIGNUM] ----
    let mut anrm = 0.0f64;
    for x in a.iter() {
        anrm = anrm.max(x.abs());
    }
    let mut cscale = 0.0f64;
    let scalea = if anrm > 0.0 && anrm < smlnum {
        cscale = smlnum;
        true
    } else if anrm > bignum {
        cscale = bignum;
        true
    } else {
        false
    };
    if scalea {
        let f = cscale / anrm;
        for x in a.iter_mut() {
            *x *= f;
        }
    }

    // ---- DGEBAL('B'): permute and scale to improve conditioning ----
    let (ilo, ihi, bal_scale) = dgebal(a, n);

    // ---- DGEHRD: A = Q H Q^T ----
    let mut q = vec![0.0; n * n];
    dgehrd(a, n, Some(&mut q));

    // ---- DHSEQR with COMPZ = 'V': A = Z T Z^T ----
    dhseqr_with_schur_vectors(a, &mut wr, &mut wi, &mut q, n)?;

    // ---- DTREVC3('B','B'): eigenvectors of T, back-transformed by Z ----
    let mut vl = q.clone();
    let mut vr = q;
    dtrevc_left(a, &mut vl, n);
    dtrevc_right(a, &mut vr, n);

    // ---- DGEBAK('B','L') and DGEBAK('B','R') ----
    dgebak_left(&mut vl, n, ilo, ihi, &bal_scale);
    dgebak_right(&mut vr, n, ilo, ihi, &bal_scale);

    // ---- normalize: unit Euclidean norm, largest component real ----
    dgeev_normalize_columns(&mut vl, n, &wi);
    dgeev_normalize_columns(&mut vr, n, &wi);

    // ---- undo the scaling of A on the eigenvalues ----
    if scalea {
        let f = anrm / cscale;
        for x in wr.iter_mut() {
            *x *= f;
        }
        for x in wi.iter_mut() {
            *x *= f;
        }
    }

    // ---- unpack the compact storages ----
    dgeev_unpack(&vl, &wi, u_real, u_imag, n);
    dgeev_unpack(&vr, &wi, v_real, v_imag, n);

    Ok((wr, wi))
}

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

#[cfg(test)]
mod tests {
    use super::super::auxiliary::dladiv;
    use super::*;
    use crate::Complex64;

    /// A deterministic, well-conditioned, non-symmetric, non-tridiagonal matrix,
    /// built from a fixed recurrence so the tests need no random generator.
    fn sample(n: usize) -> Vec<f64> {
        let mut a = vec![0.0; n * n];
        for i in 0..n {
            for j in 0..n {
                a[i + j * n] = if i == j {
                    (n as f64) + 1.0 + (i as f64) * 0.5
                } else {
                    (((i * 7 + j * 13) % 11) as f64) - 5.0
                };
            }
        }
        a
    }

    fn norm_max(a: &[f64]) -> f64 {
        a.iter().fold(0.0f64, |m, v| m.max(v.abs()))
    }

    fn mat_mul(a: &[f64], b: &[f64], n: usize) -> Vec<f64> {
        let mut c = vec![0.0; n * n];
        for j in 0..n {
            for i in 0..n {
                let mut s = 0.0;
                for l in 0..n {
                    s += a[i + l * n] * b[l + j * n];
                }
                c[i + j * n] = s;
            }
        }
        c
    }

    /// Copies a matrix into column-major storage. The LAPACK fixture files write
    /// matrices row by row, as Fortran reads them.
    fn column_major(rows: &[f64], n: usize) -> Vec<f64> {
        let mut a = vec![0.0; n * n];
        for i in 0..n {
            for j in 0..n {
                a[i + j * n] = rows[i * n + j];
            }
        }
        a
    }

    /// The 1-norm of a dense column-major matrix: LAPACK's `DLANGE('1', ...)`.
    fn norm_1(a: &[f64], m: usize, n: usize) -> f64 {
        let mut best = 0.0f64;
        for j in 0..n {
            let mut sum = 0.0;
            for i in 0..m {
                sum += a[i + j * m].abs();
            }
            best = best.max(sum);
        }
        best
    }

    fn transpose(a: &[f64], n: usize) -> Vec<f64> {
        let mut t = vec![0.0; n * n];
        for i in 0..n {
            for j in 0..n {
                t[i + j * n] = a[j + i * n];
            }
        }
        t
    }

    #[test]
    fn dlarfg_annihilates_the_tail() {
        // H x = beta e1 with H = I - tau v v^T.
        let x = [3.0, 4.0];
        let (beta, tau, v) = dlarfg(x[0], &x[1..]);
        assert!((beta.abs() - 5.0).abs() < 1e-15, "beta = {beta}");

        let mut hx = x.to_vec();
        let w: f64 = v.iter().zip(x.iter()).map(|(vi, xi)| vi * xi).sum::<f64>() * tau;
        for i in 0..2 {
            hx[i] -= w * v[i];
        }
        assert!((hx[0] - beta).abs() < 1e-14, "first component {} vs {beta}", hx[0]);
        assert!(hx[1].abs() < 1e-14, "tail must be annihilated, got {}", hx[1]);
    }

    #[test]
    fn dlarfg_of_a_zero_tail_is_the_identity() {
        let (beta, tau, _v) = dlarfg(2.5, &[0.0, 0.0]);
        assert_eq!(beta, 2.5);
        assert_eq!(tau, 0.0);
    }

    /// The underflow branch of `dlarfg`: when the tail's norm is *below*
    /// `SAFMIN` the reference rescales by it repeatedly, and that whole path --
    /// 28 lines -- was unexecuted. `1e-300` is not enough to reach it (it is
    /// larger than `SAFMIN`); `1e-310` is subnormal.
    #[test]
    fn dlarfg_rescales_a_subnormal_tail() {
        let tiny = 1e-310f64;
        let (beta, tau, v) = dlarfg(tiny, &[tiny, -tiny]);
        assert!(beta.is_finite() && beta != 0.0, "beta = {beta:e}");
        assert!(tau.is_finite() && tau != 0.0, "tau = {tau:e}");
        assert_eq!(v[0], 1.0);
        // H x = beta e1 with the tail annihilated, to the accuracy available at
        // this magnitude: every operation is subnormal.
        let x = [tiny, tiny, -tiny];
        let mut hx = x;
        let w: f64 = v.iter().zip(x.iter()).map(|(vi, xi)| vi * xi).sum::<f64>() * tau;
        for i in 0..3 {
            hx[i] -= w * v[i];
        }
        assert!((hx[0] - beta).abs() <= 4.0 * f64::MIN_POSITIVE, "first component {} vs {beta:e}", hx[0]);
        assert!(
            hx[1].abs() <= 4.0 * f64::MIN_POSITIVE && hx[2].abs() <= 4.0 * f64::MIN_POSITIVE,
            "tail not annihilated: {:?}",
            &hx[1..]
        );
    }

    #[test]
    fn dlarfg_handles_a_tiny_vector() {
        // Near the underflow threshold the naive hypot underflows to zero unless
        // the vector is rescaled first.
        let tiny = 1e-300;
        let (beta, tau, v) = dlarfg(tiny, &[tiny]);
        assert!(beta != 0.0, "beta underflowed to zero");
        assert!(tau != 0.0, "tau underflowed to zero");

        // H must still be a reflector: H x = beta e1 with the tail annihilated.
        let x = [tiny, tiny];
        let mut hx = x;
        let w: f64 = v.iter().zip(x.iter()).map(|(vi, xi)| vi * xi).sum::<f64>() * tau;
        for i in 0..2 {
            hx[i] -= w * v[i];
        }
        assert!((hx[0] - beta).abs() <= tiny.abs(), "first component {} vs {beta}", hx[0]);
        assert!(hx[1].abs() <= tiny.abs(), "tail must be annihilated, got {}", hx[1]);
    }

    /// The reduction, checked against the definition rather than against an
    /// internal factorization: `H` upper Hessenberg, `Q` orthogonal, and
    /// `A = Q H Q^T`, with the latter two in the form of LAPACK's `dhst01`.
    #[test]
    fn dgehrd_produces_hessenberg_and_an_orthogonal_q() {
        for n in [1usize, 2, 3, 5, 8, 20] {
            let orig = sample(n);
            let mut h = orig.clone();
            let mut q = vec![0.0; n * n];
            dgehrd(&mut h, n, Some(&mut q));
            for j in 0..n {
                for i in (j + 2)..n {
                    assert!(
                        h[i + j * n].abs() < 1e-13,
                        "n={n}: H({i},{j}) = {} is below the subdiagonal",
                        h[i + j * n]
                    );
                }
            }

            // The two residuals of LAPACK's `dhst01`, each divided by the norms
            // and by `n * eps`:
            //
            //   ratio 1 = norm(A - Q H Q^T) / (norm(A) * n * eps)
            //   ratio 2 = norm(I - Q^T Q)   / (n * eps)
            //
            // with the 1-norm, as the reference computes them, and LAPACK's own
            // threshold of 20. This replaces a hand-set absolute tolerance
            // (`1e-12`), which is a different assertion at n=1 and at n=20: the
            // residual of a backward-stable reduction grows with `n`, so the
            // absolute bound was silently tight at the small end and the only
            // thing holding at the large end was the choice of test matrices.
            let qh = mat_mul(&q, &h, n);
            let qhqt = mat_mul(&qh, &transpose(&q, n), n);
            let (mut diff, mut qtq_diff) = (0.0, 0.0);
            for idx in 0..n * n {
                diff += (orig[idx] - qhqt[idx]).abs();
                let (i, j) = (idx % n, idx / n);
                let mut dot = 0.0;
                for l in 0..n {
                    dot += q[l + i * n] * q[l + j * n];
                }
                qtq_diff += (dot - if i == j { 1.0 } else { 0.0 }).abs();
            }
            let eps = f64::EPSILON;
            let ratio1 = diff / (norm_1(&orig, n, n).max(f64::MIN_POSITIVE) * n as f64 * eps);
            let ratio2 = qtq_diff / (n as f64 * eps);
            assert!(ratio1 <= 20.0, "n={n}: dhst01 ratio 1 = {ratio1} exceeds 20");
            assert!(ratio2 <= 20.0, "n={n}: dhst01 ratio 2 = {ratio2} exceeds 20");
        }
    }

    /// The thirteen examples of LAPACK's `dchkbl`: `(n, A, ILO, IHI, balanced A,
    /// scale)`, the matrices row-major as the fixture files write them and
    /// `ILO`/`IHI` 1-based.
    ///
    /// The inputs are those of `TESTING/dbal.in`. The expected outputs are what
    /// LAPACK 3.12.1's own `dgebal` returns for them, not what `dbal.in` states:
    /// the file's expectation for the fifth example is stale -- it lists the
    /// permutation as `(4, 3, 5)` where the current `dgebal` gives `(4, 5, 3)`,
    /// the column phase of the permutation having been restructured after those
    /// examples were written down. `dchkbl` only *reports* its largest error
    /// instead of asserting it, and its own run over this file prints "example
    /// number having largest error = 5", which is how the drift went unnoticed.
    /// The other twelve examples agree with the file exactly.
    const DBAL_CASES: &[(usize, &[f64], usize, usize, &[f64], &[f64])] = &[
        (
            5,
            &[
                1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0,
            ],
            1,
            1,
            &[
                1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0,
            ],
            &[1.0, 2.0, 3.0, 4.0, 5.0],
        ),
        (
            5,
            &[
                1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 0.0, 0.0, 1.0, 2.0, 3.0, 4.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0,
            ],
            1,
            1,
            &[
                5.0, 4.0, 3.0, 2.0, 1.0, 0.0, 4.0, 3.0, 2.0, 1.0, 0.0, 0.0, 3.0, 2.0, 1.0, 0.0, 0.0, 0.0, 2.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0,
            ],
            &[1.0, 2.0, 3.0, 2.0, 1.0],
        ),
        (
            5,
            &[
                1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0,
            ],
            1,
            1,
            &[
                1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0,
            ],
            &[1.0, 2.0, 3.0, 2.0, 1.0],
        ),
        (
            4,
            &[0.0, 2.0, 0.1, 0.0, 2.0, 0.0, 0.0, 0.1, 100.0, 0.0, 0.0, 2.0, 0.0, 100.0, 2.0, 0.0],
            1,
            4,
            &[0.0, 2.0, 3.2, 0.0, 2.0, 0.0, 0.0, 3.2, 3.125, 0.0, 0.0, 2.0, 0.0, 3.125, 2.0, 0.0],
            &[0.0625, 0.0625, 2.0, 2.0],
        ),
        (
            6,
            &[
                2.0, 0.0, 0.0, 0.0, 0.0, 1024.0, 0.0, 0.0, 0.0, 0.0, 0.0, 128.0, 0.0, 2.0, 3000.0, 0.0, 0.0, 2.0, 128.0, 4.0, 0.004, 5.0, 600.0, 8.0,
                0.0, 0.0, 0.0, 0.0, 0.002, 2.0, 8.0, 8192.0, 0.0, 0.0, 0.0, 2.0,
            ],
            4,
            6,
            &[
                5.0, 600.0, 0.004, 1024.0, 0.5, 8.0, 0.0, 0.002, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 3000.0, 0.0, 0.25, 2.0, 0.0, 0.0, 0.0, 2.0, 0.0,
                128.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1024.0, 0.0, 0.0, 0.0, 64.0, 1024.0, 2.0,
            ],
            &[4.0, 5.0, 3.0, 8.0, 0.125, 1.0],
        ),
        (
            5,
            &[
                1.0, 0.0, 0.0, 0.0, 8.0, 0.0, 2.0, 8192.0, 2.0, 4.0, 0.00025, 0.000125, 4.0, 0.0, 64.0, 0.0, 2.0, 1024.0, 4.0, 8.0, 0.0, 8192.0, 0.0,
                0.0, 8.0,
            ],
            1,
            5,
            &[
                1.0, 0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 1024.0, 16.0, 16.0, 0.032, 0.001, 4.0, 0.0, 2048.0, 0.0, 0.25, 16.0, 4.0, 4.0, 0.0, 2048.0, 0.0,
                0.0, 8.0,
            ],
            &[8.0, 0.5, 0.0625, 4.0, 2.0],
        ),
        (
            4,
            &[
                1.0, 1000000.0, 1000000.0, 1000000.0, -2000000.0, 3.0, 2e-06, 3e-06, -3000000.0, 0.0, 1e-06, 2.0, 1000000.0, 0.0, 3e-06, 4000000.0,
            ],
            1,
            4,
            &[
                1.0, 1000000.0, 2000000.0, 1000000.0, -2000000.0, 3.0, 4e-06, 3e-06, -1500000.0, 0.0, 1e-06, 1.0, 1000000.0, 0.0, 6e-06, 4000000.0,
            ],
            &[1.0, 1.0, 2.0, 1.0],
        ),
        (
            4,
            &[
                1.0, 10000.0, 10000.0, 10000.0, -20000.0, 3.0, 0.002, 0.003, 0.0, 2.0, 0.0, -30000.0, 0.0, 0.0, 10000.0, 0.0,
            ],
            1,
            4,
            &[
                1.0, 10000.0, 10000.0, 5000.0, -20000.0, 3.0, 0.002, 0.0015, 0.0, 2.0, 0.0, -15000.0, 0.0, 0.0, 20000.0, 0.0,
            ],
            &[1.0, 1.0, 1.0, 0.5],
        ),
        (
            5,
            &[
                1.0, 512.0, 4096.0, 32768.0, 262144.0, 8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0,
                0.0,
            ],
            1,
            5,
            &[
                1.0, 32.0, 32.0, 32.0, 32.0, 128.0, 0.0, 0.0, 0.0, 0.0, 0.0, 64.0, 0.0, 0.0, 0.0, 0.0, 0.0, 64.0, 0.0, 0.0, 0.0, 0.0, 0.0, 64.0, 0.0,
            ],
            &[256.0, 16.0, 2.0, 0.25, 0.03125],
        ),
        (
            6,
            &[
                1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0,
                0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0,
            ],
            2,
            5,
            &[
                1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0,
                1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
            ],
            &[3.0, 1.0, 1.0, 1.0, 1.0, 4.0],
        ),
        (
            7,
            &[
                6.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 4.0, 0.0, 0.00025, 0.0125, 0.02, 0.125, 1.0, 128.0, 64.0, 0.0, 0.0, -2.0, 16.0, 0.0, 16384.0,
                0.0, 1.0, -400.0, 256.0, -4000.0, -2.0, -256.0, 0.0, 0.0125, 2.0, 2.0, 32.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0, 0.0, 0.004,
                0.125, -0.2, 3.0,
            ],
            2,
            5,
            &[
                64.0, 1.0, 0.5, 0.0, 0.0, 1.0, -2.0, 0.0, 4.0, 0.5, 1.024, 0.8, 0.0, 2.56, 0.0, 2.0, 3.0, 4.096, 2.0, 0.0, -6.4, 0.0, 4.0, -3.90625,
                1.0, -6.25, 0.0, 8.0, 0.0, -4.0, 2.0, 0.8, 2.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            ],
            &[3.0, 0.0078125, 0.03125, 32.0, 0.5, 1.0, 6.0],
        ),
        (
            5,
            &[
                1000.0, 2.0, 3.0, 4.0, 500000.0, 9.0, 0.0, 0.0002, 1.0, 3.0, 0.0, -300.0, 2.0, 1.0, 1.0, 9.0, 0.002, 1.0, 1.0, -1000.0, 6.0, 200.0,
                1.0, 600.0, 3.0,
            ],
            1,
            5,
            &[
                1000.0, 0.03125, 0.375, 0.03125, 1953.125, 576.0, 0.0, 0.0016, 0.5, 0.75, 0.0, -37.5, 2.0, 0.0625, 0.03125, 1152.0, 0.004, 16.0, 1.0,
                -500.0, 1536.0, 800.0, 32.0, 1200.0, 3.0,
            ],
            &[32.0, 0.5, 4.0, 0.25, 0.125],
        ),
        (
            6,
            &[
                1.0, 1e+120, 0.0, 0.0, 0.0, 0.0, 1e-120, 1.0, 1e+120, 0.0, 0.0, 0.0, 0.0, 1e-120, 1.0, 1e+120, 0.0, 0.0, 0.0, 0.0, 1e-120, 1.0,
                1e+120, 0.0, 0.0, 0.0, 0.0, 1e-120, 1.0, 1e+120, 0.0, 0.0, 0.0, 0.0, 1e-120, 1.0,
            ],
            1,
            6,
            &[
                1.0,
                6344.854593289123,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0001576080247855779,
                1.0,
                6344.854593289123,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0001576080247855779,
                1.0,
                3172.4272966445615,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0003152160495711558,
                1.0,
                1586.2136483222807,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0006304320991423117,
                1.0,
                793.1068241611404,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0012608641982846233,
                1.0,
            ],
            &[
                2.4948003869184e+291,
                1.582914569427869e+175,
                1.004336277661869e+59,
                3.1861838222649046e-58,
                5.053968264940244e-175,
                4.008336720017946e-292,
            ],
        ),
    ];

    /// The examples of LAPACK's `dchkbk`, transcribed from `TESTING/dbak.in`:
    /// `(n, ILO, IHI, scale, V, V after dgebak)`. `dchkbk` exercises `SIDE = 'R'`
    /// only, which is the side `dgeev` uses.
    const DBAK_CASES: &[(usize, usize, usize, &[f64], &[f64], &[f64])] = &[
        (
            5,
            1,
            1,
            &[1.0, 2.0, 3.0, 4.0, 5.0],
            &[
                1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
            ],
            &[
                1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
            ],
        ),
        (
            5,
            1,
            1,
            &[1.0, 2.0, 3.0, 2.0, 1.0],
            &[
                1.0, 1.0, 1.0, -0.6667, -0.04167, 0.0, -0.25, -0.6667, 1.0, 0.1667, 0.0, 0.0, 0.2222, -1.0, -0.5, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0,
                0.0, 0.0, -1.0,
            ],
            &[
                0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.2222, -1.0, -0.5, 0.0, -0.25, -0.6667, 1.0, 0.1667, 1.0, 1.0, 1.0,
                -0.6667, -0.04167,
            ],
        ),
        (
            5,
            1,
            1,
            &[1.0, 2.0, 3.0, 2.0, 1.0],
            &[
                1.0, 1.0, 1.0, 1.0, 1.0, 0.0, -6e-18, -6e-18, -6e-18, -6e-18, 0.0, 0.0, 3.6e-35, 3.6e-35, 3.6e-35, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                0.0, 0.0, 0.0,
            ],
            &[
                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.6e-35, 3.6e-35, 3.6e-35, 0.0, -6e-18, -6e-18, -6e-18, -6e-18, 1.0, 1.0,
                1.0, 1.0, 1.0,
            ],
        ),
        (
            6,
            4,
            6,
            &[4.0, 3.0, 5.0, 100.0, 0.1, 1.0],
            &[
                1.0, 1.336e-06, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, -3.001e-11, -3.252e-05, 0.01305, 0.0, 0.0, -0.00833, 8.929e-10, -6.712e-05,
                6.687e-05, 0.0, 0.0, 0.0, -4.455e-06, -0.003355, 0.003345, 0.0, 0.0, 0.0, 4.455e-07, -0.03356, 0.03344, 0.0, 0.0, 0.0, 4.411e-10,
                0.1011, 0.1008,
            ],
            &[
                0.0, 0.0, 0.0, -0.0004455, -0.3355, 0.3345, 0.0, 0.0, 0.0, 4.455e-08, -0.003356, 0.003344, 0.0, 1.0, 0.0, -3.001e-11, -3.252e-05,
                0.01305, 1.0, 1.336e-06, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, -0.00833, 8.929e-10, -6.712e-05, 6.687e-05, 0.0, 0.0, 0.0, 4.411e-10, 0.1011,
                0.1008,
            ],
        ),
        (
            5,
            1,
            5,
            &[100.0, 0.1, 0.01, 1.0, 10.0],
            &[
                0.0001366, -6.829e-05, 0.0001252, 1.0, 1.95e-15, 1.0, 1.0, -2.776e-17, 3.601e-06, -6.073e-18, 0.2736, -0.1363, 0.2503, -3.322e-06,
                -0.002, 0.006909, -0.003443, 0.006196, 0.01666, 1.0, 0.3899, -0.2033, -0.342, -0.001, 6e-15,
            ],
            &[
                0.01366, -0.006829, 0.01252, 100.0, 1.95e-13, 0.1, 0.1, -2.776e-18, 3.601e-07, -6.073e-19, 0.002736, -0.001363, 0.002503, -3.322e-08,
                -2e-05, 0.006909, -0.003443, 0.006196, 0.01666, 1.0, 3.899, -2.033, -3.42, -0.01, 6e-14,
            ],
        ),
        (
            6,
            2,
            5,
            &[3.0, 1.0, 1.0, 1.0, 1.0, 4.0],
            &[
                1.0, 1.0, 2.776e-16, -2.405e-17, 0.0, 1.0, 0.0, 0.75, 1.0, 0.0852, 0.0, -1.52e-17, 0.0, 0.75, -0.8093, 1.0, 0.0, -1.52e-17, 0.0,
                0.75, -0.09533, -0.5426, 1.0, -1.52e-17, 0.0, 0.75, -0.09533, -0.5426, -1.0, -1.52e-17, 0.0, 0.0, 0.0, 0.0, 0.0, 4.559e-17,
            ],
            &[
                0.0, 0.75, -0.8093, 1.0, 0.0, -1.52e-17, 0.0, 0.75, 1.0, 0.0852, 0.0, -1.52e-17, 1.0, 1.0, 2.776e-16, -2.405e-17, 0.0, 1.0, 0.0, 0.0,
                0.0, 0.0, 0.0, 4.559e-17, 0.0, 0.75, -0.09533, -0.5426, -1.0, -1.52e-17, 0.0, 0.75, -0.09533, -0.5426, 1.0, -1.52e-17,
            ],
        ),
        (
            7,
            2,
            5,
            &[3.0, 0.001, 0.01, 10.0, 0.1, 1.0, 6.0],
            &[
                1.0, -0.01105, 0.03794, -0.09378, -0.03481, 0.4465, -0.03602, 0.0, -0.4556, -0.4545, 1.0, 0.4639, -0.6512, 0.4781, 0.0, -0.2734,
                -0.7946, 0.6303, 1.0, -0.6279, 1.0, 0.0, 1.0, -6.939e-18, 0.04259, -0.6495, -0.5581, -0.6452, 0.0, -0.3904, -0.4029, -0.1685,
                -0.9429, 1.0, -0.9371, 0.0, 0.0, 0.0, 0.0, 0.0, -0.2558, 0.0003308, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.001985,
            ],
            &[
                0.0, 0.0, 0.0, 0.0, 0.0, -0.2558, 0.0003308, 0.0, -0.0004556, -0.0004545, 0.001, 0.0004639, -0.0006512, 0.0004781, 1.0, -0.01105,
                0.03794, -0.09378, -0.03481, 0.4465, -0.03602, 0.0, 10.0, -6.939e-17, 0.4259, -6.495, -5.581, -6.452, 0.0, -0.03904, -0.04029,
                -0.01685, -0.09429, 0.1, -0.09371, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.001985, 0.0, -0.002734, -0.007946, 0.006303, 0.01, -0.006279,
                0.01,
            ],
        ),
    ];

    /// LAPACK's `dchkbl` examples: for each matrix, `dgebal` must return the
    /// permutation (`ILO`, `IHI`), the balanced matrix and the scaling vector
    /// that the reference does.
    ///
    /// `dgebal` had no test of its own before this one: it ran only inside
    /// `dgeev`, where a wrong permutation or scale surfaces as a wrong
    /// eigenvector several hundred lines away. Permuting and scaling are exact --
    /// row/column interchanges and powers of two -- so the comparison is exact up
    /// to the rounding of a dot product, hence the `1e-13` relative bound.
    #[test]
    fn dgebal_reproduces_the_reference_examples() {
        for (k, &(n, a_in, ilo_in, ihi_in, a_ref_rows, scale_ref)) in DBAL_CASES.iter().enumerate() {
            let case = k + 1;
            let mut a = column_major(a_in, n);
            let a_ref = column_major(a_ref_rows, n);
            let (ilo, ihi, scale) = dgebal(&mut a, n);
            assert_eq!((ilo + 1, ihi + 1), (ilo_in, ihi_in), "case {case} (n={n}): ILO/IHI disagree with dbal.in");
            let norm = a_ref.iter().fold(0.0f64, |m, v| m.max(v.abs())).max(1.0);
            for j in 0..n {
                for i in 0..n {
                    let (got, want) = (a[i + j * n], a_ref[i + j * n]);
                    assert!(
                        (got - want).abs() <= 1e-13 * norm,
                        "case {case} (n={n}): balanced A({i},{j}) = {got}, the reference gives {want}"
                    );
                }
            }
            for i in 0..n {
                let (got, want) = (scale[i], scale_ref[i]);
                assert!(
                    (got - want).abs() <= 1e-13 * want.abs().max(1.0),
                    "case {case} (n={n}): scale[{i}] = {got}, the reference gives {want}"
                );
            }
        }
    }

    /// The complex twin must reproduce the real one. Balancing is a permutation
    /// and powers of two, so for a real matrix stored as complex the two routines
    /// have to agree entry for entry -- and `zgebal` had no test at all.
    #[test]
    fn zgebal_matches_dgebal_on_the_reference_examples() {
        for (k, &(n, a_in, ilo_in, ihi_in, _, _)) in DBAL_CASES.iter().enumerate() {
            let case = k + 1;
            let a_real = column_major(a_in, n);
            let mut a_complex: Vec<Complex64> = a_real.iter().map(|&v| Complex64::new(v, 0.0)).collect();
            let mut a_work = a_real.clone();
            let (ilo, ihi, scale) = dgebal(&mut a_work, n);
            let (zilo, zihi, zscale) = crate::lapack::zeig::zgebal(&mut a_complex, n);
            assert_eq!((zilo, zihi), (ilo, ihi), "case {case} (n={n}): zgebal disagrees with dgebal on ILO/IHI");
            assert_eq!(
                (zilo + 1, zihi + 1),
                (ilo_in, ihi_in),
                "case {case} (n={n}): both must agree with dbal.in"
            );
            let norm = a_work.iter().fold(0.0f64, |m, v| m.max(v.abs())).max(1.0);
            for j in 0..n {
                for i in 0..n {
                    let want = a_work[i + j * n];
                    let got = a_complex[i + j * n].re;
                    assert!(
                        (got - want).abs() <= 1e-13 * norm,
                        "case {case} (n={n}): zgebal gives {got} at ({i},{j}), dgebal gives {want}"
                    );
                    assert_eq!(a_complex[i + j * n].im, 0.0, "case {case}: a real matrix must stay real");
                }
            }
            assert_eq!(zscale, scale, "case {case} (n={n}): scale vectors differ");
        }
    }

    /// `dgebak_left` (`SIDE = 'L'`) had no test: `dgeev` uses the right one, and
    /// the left one is reached only through `dgeev_lr`.
    ///
    /// The left back-transformation scales *rows* by `1/scale` where the right one
    /// scales columns by `scale`, and then applies the same permutation. The
    /// expected matrix is built here from that definition, and also for a negative
    /// scale factor, which LAPACK allows.
    #[test]
    fn dgebak_left_scales_rows_and_applies_the_permutation() {
        // ILO = 0, IHI = 4 (0-based) for n = 6: inside the block the scale entries
        // are factors, outside they are 1-based permutation indices.
        let n = 6;
        let (ilo, ihi) = (1usize, 4usize);
        let scale = vec![3.0, 0.5, 2.0, 0.25, 8.0, 6.0];
        let mut v = vec![0.0; n * n];
        for i in 0..n {
            for j in 0..n {
                v[i + j * n] = (i as f64) + 0.5 * (j as f64) + 1.0;
            }
        }

        // Reference definition: rows ILO..IHI are divided by SCALE(I), then rows
        // are permuted, I = ILO-1 down to 1 first and I = IHI+1 up to N after.
        let expected = |v: &[f64], scale: &[f64], ilo: usize, ihi: usize| -> Vec<f64> {
            let mut w = v.to_vec();
            // The reference guards the scaling with `ILO != IHI`.
            if ilo != ihi {
                for i in ilo..=ihi {
                    for j in 0..n {
                        w[i + j * n] /= scale[i];
                    }
                }
            }
            for i in (0..ilo).rev() {
                let k = scale[i] as usize - 1;
                for j in 0..n {
                    w.swap(i + j * n, k + j * n);
                }
            }
            for i in (ihi + 1)..n {
                let k = scale[i] as usize - 1;
                for j in 0..n {
                    w.swap(i + j * n, k + j * n);
                }
            }
            w
        };

        for case_scale in [scale.clone(), {
            // A negative factor is legal and must be applied as given, not as
            // its magnitude.
            let mut s = scale.clone();
            s[2] = -2.0;
            s
        }] {
            let want = expected(&v, &case_scale, ilo, ihi);
            let mut got = v.clone();
            dgebak_left(&mut got, n, ilo, ihi, &case_scale);
            for idx in 0..n * n {
                assert!(
                    (got[idx] - want[idx]).abs() <= 1e-13 * want[idx].abs().max(1.0),
                    "scale = {case_scale:?}, index {idx}: dgebak_left gives {}, the definition gives {}",
                    got[idx],
                    want[idx]
                );
            }
        }

        // With ILO == IHI the scaling loop is skipped and only the permutation
        // runs: the reference tests `ILO != IHI` for it.
        let tight = vec![4.0, 3.0, 1.0, 1.0, 1.0, 6.0];
        let mut got = v.clone();
        dgebak_left(&mut got, n, 2, 2, &tight);
        let want = expected(&v, &tight, 2, 2);
        for idx in 0..n * n {
            assert!(
                (got[idx] - want[idx]).abs() <= 1e-13 * want[idx].abs().max(1.0),
                "ILO == IHI, index {idx}: {} != {}",
                got[idx],
                want[idx]
            );
        }
    }

    /// LAPACK's `dchkbk` examples: `dgebak_right` must undo the balancing on the
    /// eigenvectors held in `V`, recovering the vectors in `TESTING/dbak.in`.
    ///
    /// `dgebak_right` and `dgebak_left` likewise had no direct test; both are
    /// applied to every eigenvector `dgeev` returns.
    #[test]
    fn dgebak_undoes_the_reference_balancing() {
        for (k, &(n, ilo, ihi, scale, v_in, v_ref_rows)) in DBAK_CASES.iter().enumerate() {
            let case = k + 1;
            let mut v = column_major(v_in, n);
            let v_ref = column_major(v_ref_rows, n);
            dgebak_right(&mut v, n, ilo - 1, ihi - 1, scale);
            let norm = v_ref.iter().fold(0.0f64, |m, x| m.max(x.abs())).max(1.0);
            for j in 0..n {
                for i in 0..n {
                    let (got, want) = (v[i + j * n], v_ref[i + j * n]);
                    assert!(
                        (got - want).abs() <= 1e-13 * norm,
                        "case {case} (n={n}): V({i},{j}) = {got}, dbak.in has {want}"
                    );
                }
            }
        }
    }

    /// `dlartg`: the plane rotation must satisfy `[c s; -s c] (f, g) = (r, 0)`,
    /// for `f` and `g` in the normal range and in the scaled one.
    #[test]
    fn dlartg_rotates_a_vector_onto_the_first_axis() {
        let cases = [
            (3.0, 4.0),
            (-3.0, 4.0),
            (3.0, -4.0),
            (0.0, 5.0),
            (5.0, 0.0),
            (0.0, -5.0),
            (1e300, 1e300),   // above rtmax: the scaling branch
            (1e-300, 1e-300), // below rtmin: the scaling branch
            (1e300, 1.0),
            (-2.5, 0.0),
        ];
        for (f, g) in cases {
            let (c, s, r) = dlartg(f, g);
            let scale = f.abs().max(g.abs()).max(f64::MIN_POSITIVE);
            assert!((c * c + s * s - 1.0).abs() <= 1e-14, "f = {f:e}, g = {g:e}: rotation not unit ({c},{s})");
            let r1 = c * f + s * g;
            let r2 = -s * f + c * g;
            assert!(
                (r1 - r).abs() <= 1e-13 * scale && r2.abs() <= 1e-13 * scale,
                "f = {f:e}, g = {g:e}: rotation gives ({r1},{r2}) not ({r},0)"
            );
            // `f.hypot(g)` is the reference norm; squaring first underflows here.
            assert!(
                (r.abs() - f.hypot(g)).abs() <= 1e-12 * f.hypot(g) && r.is_finite(),
                "f = {f:e}, g = {g:e}: |r| = {} is not the norm {}",
                r.abs(),
                f.hypot(g)
            );
        }
        // g == 0 and f == 0 are the exact identities.
        assert_eq!(dlartg(2.0, 0.0), (1.0, 0.0, 2.0));
        assert_eq!(dlartg(0.0, -2.0), (0.0, -1.0, 2.0));
    }

    /// `dladiv`: `(a + ib) / (c + id)`, including a divisor whose squares
    /// overflow.
    #[test]
    fn dladiv_divides_without_overflowing() {
        let cases = [
            (1.0, 2.0, 3.0, 4.0),
            (-1.5, 0.25, 2.0, -8.0),
            (1.0, 0.0, 3.0, 4.0),
            (0.0, 1.0, 0.0, 2.0),
            (5.0, 0.0, 0.0, 2.0), // c == 0: the `aa *= 0.5` branch
        ];
        for (a, b, c, d) in cases {
            let (p, q) = dladiv(a, b, c, d);
            let den = c * c + d * d;
            let (want_p, want_q) = ((a * c + b * d) / den, (b * c - a * d) / den);
            let scale = want_p.abs().max(want_q.abs()).max(1.0);
            assert!(
                (p - want_p).abs() <= 1e-13 * scale && (q - want_q).abs() <= 1e-13 * scale,
                "({a},{b})/({c},{d}): ({p},{q}) != ({want_p},{want_q})"
            );
        }
        // Both parts large: the naive formula overflows.
        let (p, q) = dladiv(1e300, 0.0, 1e300, 1e300);
        assert!(p.is_finite() && q.is_finite(), "({p},{q})");
        assert!((p - 0.5).abs() < 1e-14 && (q + 0.5).abs() < 1e-14, "({p},{q})");
    }

    /// `dlaln2`: the four shapes it is called with. Each must solve
    /// `(A - (wr + i wi) I) X = scale * B`, with `scale == 1` when the data is
    /// well conditioned.
    #[test]
    fn dlaln2_solves_the_four_block_shapes() {
        // (a11, b11) with a real eigenvalue: x = b11 / (a11 - wr).
        let r = dlaln2_real(1, 1, 0.0, 4.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 1.0, 0.0);
        assert_eq!(r.scale, 1.0);
        assert!((r.x11 - 2.0 / 3.0).abs() < 1e-15, "x11 = {}", r.x11);

        // Complex 1x1: x = (b11 + i b12) / (a11 - wr - i wi).
        let (a11, b11, b12, wr, wi) = (4.0, 1.0, 2.0, 1.0, 0.5);
        let r = dlaln2_real(1, 2, 0.0, a11, 0.0, 0.0, 0.0, b11, 0.0, b12, 0.0, wr, wi);
        assert_eq!(r.scale, 1.0);
        // (1 + 2i) / (3 - 0.5i) = (1*3 + 2*(-0.5) + i(2*3 - 1*(-0.5)))/9.25
        let den = (a11 - wr) * (a11 - wr) + wi * wi;
        // (b11 + i b12) * ((a11 - wr) + i wi) / den
        let want_re = (b11 * (a11 - wr) - b12 * wi) / den;
        let want_im = (b12 * (a11 - wr) + b11 * wi) / den;
        assert!((r.x11 - want_re).abs() < 1e-14, "x11 = {} want {want_re}", r.x11);
        assert!((r.x12 - want_im).abs() < 1e-14, "x12 = {} want {want_im}", r.x12);

        // Real 2x2: (A - wr I) x = (b11, b21).
        let (a11, a21, a12, a22) = (4.0, 2.0, 1.0, 3.0);
        let (b11, b21, wr) = (1.0, -2.0, 0.5);
        let r = dlaln2_real(2, 1, 0.0, a11, a21, a12, a22, b11, b21, 0.0, 0.0, wr, 0.0);
        assert_eq!(r.scale, 1.0);
        let det = (a11 - wr) * (a22 - wr) - a12 * a21;
        let want_x11 = ((a22 - wr) * b11 - a12 * b21) / det;
        let want_x21 = (-a21 * b11 + (a11 - wr) * b21) / det;
        assert!((r.x11 - want_x11).abs() < 1e-14, "x11 = {}", r.x11);
        assert!((r.x21 - want_x21).abs() < 1e-14, "x21 = {}", r.x21);

        // Complex 2x2: (A - wr I - i wi I) x = (b11 + i b12, b21 + i b22).
        let (wr, wi) = (0.5, 1.5);
        let (b11, b21, b12, b22) = (1.0, -1.0, 0.5, 2.0);
        let r = dlaln2_real(2, 2, 0.0, a11, a21, a12, a22, b11, b21, b12, b22, wr, wi);
        assert_eq!(r.scale, 1.0);
        // Solve with Cramer's rule over the complex field.
        let m11 = Complex64::new(a11 - wr, -wi);
        let m12 = Complex64::new(a12, 0.0);
        let m21 = Complex64::new(a21, 0.0);
        let m22 = Complex64::new(a22 - wr, -wi);
        let rhs1 = Complex64::new(b11, b12);
        let rhs2 = Complex64::new(b21, b22);
        let det = m11 * m22 - m12 * m21;
        let want1 = (rhs1 * m22 - m12 * rhs2) / det;
        let want2 = (m11 * rhs2 - rhs1 * m21) / det;
        let got1 = Complex64::new(r.x11, r.x12);
        let got2 = Complex64::new(r.x21, r.x22);
        assert!((got1 - want1).norm() < 1e-13, "x1 = {got1:?} want {want1:?}");
        assert!((got2 - want2).norm() < 1e-13, "x2 = {got2:?} want {want2:?}");
    }

    /// `dlaln2`'s rescaling: a denominator near the underflow limit and a large
    /// right-hand side must come back with `scale < 1` rather than `inf`, and a
    /// large `SMIN` must clamp the denominator. Both branches leave the defining
    /// relation `(A - wr I) X = scale B` intact, so that is what is asserted.
    #[test]
    fn dlaln2_rescales_instead_of_overflowing() {
        // 1x1 with a denominator 1e-300 and a right-hand side of 1e300: the exact
        // solution would be 1e600, so the routine must return
        // `x = 1e300` with `scale = 1e-300`.
        let (a11, b11) = (1e-300, 1e300);
        let r = dlaln2_real(1, 1, 0.0, a11, 0.0, 0.0, 0.0, b11, 0.0, 0.0, 0.0, 0.0, 0.0);
        assert!(r.scale < 1.0, "scale = {} should be below one", r.scale);
        assert!(r.x11.is_finite(), "x11 = {} overflowed", r.x11);
        let lhs = r.x11 * a11;
        let rhs = r.scale * b11;
        assert!(
            (lhs - rhs).abs() <= 1e-13 * rhs.abs().max(f64::MIN_POSITIVE),
            "({lhs:e}) != scale * b ({rhs:e})"
        );

        // A large SMIN clamps the denominator, so the solution is bounded by
        // b/smin rather than by the true (tiny) denominator.
        let smin = 1e-3;
        let r = dlaln2_real(1, 1, smin, 1e-8, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0);
        assert!(r.x11.abs() <= 1.0 / smin * 1.000001, "x11 = {} exceeds b/smin = {}", r.x11, 1.0 / smin);
        assert!(r.x11.is_finite());
    }

    /// A lower-triangular block must be swapped (`b == 0`), and a block whose
    /// eigenvalues are real but nearly equal takes the equalization branch.
    #[test]
    fn dlanv2_handles_the_triangular_and_nearly_equal_blocks() {
        // [2 0; 3 1]: lower triangular, so the code exchanges a and d.
        let (mut a, mut b, mut c, mut d) = (2.0, 0.0, 3.0, 1.0);
        let (mut cs, mut sn) = (0.0, 0.0);
        let e = dlanv2(&mut a, &mut b, &mut c, &mut d, &mut cs, &mut sn);
        assert_eq!(c, 0.0, "the block must end up upper triangular");
        assert_eq!((a, d), (1.0, 2.0), "the diagonal must be exchanged");
        assert_eq!((cs, sn), (0.0, 1.0), "the rotation must be the exchange");
        assert_eq!((e.rt1i, e.rt2i), (0.0, 0.0));
        assert_eq!((e.rt1r, e.rt2r), (1.0, 2.0));

        // Already upper triangular: the identity rotation.
        let (mut a, mut b, mut c, mut d) = (1.0, 2.0, 0.0, 3.0);
        let (mut cs, mut sn) = (0.0, 0.0);
        let e = dlanv2(&mut a, &mut b, &mut c, &mut d, &mut cs, &mut sn);
        assert_eq!((cs, sn), (1.0, 0.0));
        assert_eq!((e.rt1r, e.rt2r), (1.0, 3.0));

        // Real, nearly equal eigenvalues with b and c of the same sign: the
        // equalization loop runs and the block is reduced to triangular form.
        let (mut a, mut b, mut c, mut d) = (1.0, 1e-9, 2e-9, 1.0);
        let (mut cs, mut sn) = (0.0, 0.0);
        let e = dlanv2(&mut a, &mut b, &mut c, &mut d, &mut cs, &mut sn);
        assert_eq!(e.rt1i, 0.0, "the eigenvalues are real");
        assert_eq!(e.rt2i, 0.0);
        assert!((cs * cs + sn * sn - 1.0).abs() < 1e-14);

        // Extreme magnitudes: the loop rescales instead of overflowing.
        let (mut a, mut b, mut c, mut d) = (1e-200, 1e-200, -1e-200, 1e-200);
        let (mut cs, mut sn) = (0.0, 0.0);
        let e = dlanv2(&mut a, &mut b, &mut c, &mut d, &mut cs, &mut sn);
        assert!((cs * cs + sn * sn - 1.0).abs() < 1e-14);
        assert!(
            e.rt1r.is_finite() && e.rt2r.is_finite() && e.rt1i.is_finite(),
            "eigenvalues must stay finite: {e:?}"
        );
    }

    /// `dgeev`'s eigenvectors are normalized the way LAPACK's `dget22` checks:
    /// unit Euclidean norm, and for a conjugate pair the component of largest
    /// modulus lies on the real axis.
    #[test]
    fn dgeev_returns_normalized_eigenvectors() {
        for n in [2usize, 3, 5, 8] {
            let orig = sample(n);
            let mut a = orig.clone();
            let (mut vr, mut vi) = (vec![0.0; n * n], vec![0.0; n * n]);
            let (wr, wi) = dgeev(&mut vr, &mut vi, &mut a, n).expect("dgeev");
            let mut j = 0;
            while j < n {
                if wi[j] == 0.0 {
                    // Real eigenvector: norm 1, no imaginary part.
                    let norm = (0..n).map(|i| vr[i + j * n] * vr[i + j * n]).sum::<f64>().sqrt();
                    assert!((norm - 1.0).abs() < 1e-13, "n={n}, j={j}: norm = {norm}");
                    for i in 0..n {
                        assert_eq!(vi[i + j * n], 0.0, "n={n}: real column {j} has an imaginary part");
                    }
                    j += 1;
                } else {
                    // A pair: v = vr[:,j] + i vi[:,j] has norm 1, and the
                    // component of largest modulus of v is real.
                    let norm = (0..n)
                        .map(|i| vr[i + j * n] * vr[i + j * n] + vi[i + j * n] * vi[i + j * n])
                        .sum::<f64>()
                        .sqrt();
                    assert!((norm - 1.0).abs() < 1e-13, "n={n}, pair {j}: norm = {norm}");
                    let (mut imax, mut best) = (0, 0.0);
                    for i in 0..n {
                        let m = vr[i + j * n].hypot(vi[i + j * n]);
                        if m > best {
                            best = m;
                            imax = i;
                        }
                    }
                    assert!(
                        vi[imax + j * n].abs() < 1e-12 * best.max(1.0),
                        "n={n}, pair {j}: the largest component (row {imax}) is not real \
                         ({} + {}i)",
                        vr[imax + j * n],
                        vi[imax + j * n]
                    );
                    assert_eq!((wr.len(), wi.len()), (n, n));
                    j += 2;
                }
            }
        }
    }

    #[test]
    fn dlanv2_real_eigenvalues() {
        // [ 1 2 ; 3 4 ] has the real eigenvalues (5 +- sqrt(33))/2.
        let (mut a, mut b, mut c, mut d) = (1.0, 2.0, 3.0, 4.0);
        let (mut cs, mut sn) = (0.0, 0.0);
        let e = dlanv2(&mut a, &mut b, &mut c, &mut d, &mut cs, &mut sn);
        assert_eq!(e.rt1i, 0.0);
        assert_eq!(e.rt2i, 0.0);
        assert!(c.abs() < 1e-14, "the block must be made upper triangular, c = {c}");
        let disc = 33.0f64.sqrt();
        for want in [(5.0 + disc) / 2.0, (5.0 - disc) / 2.0] {
            assert!(
                [e.rt1r, e.rt2r].iter().any(|g| (g - want).abs() < 1e-13),
                "missing eigenvalue {want} in {:?}",
                [e.rt1r, e.rt2r]
            );
        }
    }

    #[test]
    fn dlanv2_complex_pair() {
        // [ 0 1 ; -1 0 ] has eigenvalues +-i.
        let (mut a, mut b, mut c, mut d) = (0.0, 1.0, -1.0, 0.0);
        let (mut cs, mut sn) = (0.0, 0.0);
        let e = dlanv2(&mut a, &mut b, &mut c, &mut d, &mut cs, &mut sn);
        assert!(e.rt1i > 0.0, "expected a complex pair with a positive rt1i, got {e:?}");
        assert!((e.rt1i + e.rt2i).abs() < 1e-14, "the pair must be conjugate");
        assert!(e.rt1r.abs() < 1e-14);
        assert!((e.rt1i - 1.0).abs() < 1e-14, "|imag| should be 1, got {}", e.rt1i);
        // The standard form requires equal diagonal entries and B*C < 0.
        assert!((a - d).abs() < 1e-14, "the diagonal must be equal: {a} vs {d}");
        assert!(b * c < 0.0, "B*C must be negative, got {} * {}", b, c);
    }

    /// The 3x3 case against an independently built `Hr A Hr`.
    ///
    /// For `n = 3` there is exactly one step, so `H` must equal `Hr * A * Hr`
    /// where `Hr` is the 3x3 reflector that is the identity on row/column 0 and
    /// `I - tau v v^T` on rows/columns 1..2. This pins the exact numbers rather
    /// than only the invariants, so a "reduction" that is similar to `A` but not
    /// the reduction LAPACK performs cannot pass.
    #[test]
    fn dgehrd_matches_the_3x3_reference() {
        // Column-major.
        let a0 = [4.0, 2.0, -2.0, -3.0, 4.5, 0.0, -1.0, -5.0, 5.0];
        let n = 3usize;

        // Hr, built directly from dlarfg on the column below the first row.
        let (beta, tau, v) = dlarfg(a0[1], &[a0[2]]);
        assert!((beta + 2.8284271247461903).abs() < 1e-14, "beta = {beta}");
        let mut hr = vec![0.0; n * n];
        for i in 0..n {
            for j in 0..n {
                hr[i + j * n] = if i == 0 && j == 0 {
                    1.0
                } else if i > 0 && j > 0 {
                    (if i == j { 1.0 } else { 0.0 }) - tau * v[i - 1] * v[j - 1]
                } else {
                    0.0
                };
            }
        }

        let mul = |x: &[f64], y: &[f64]| {
            let mut c = vec![0.0; n * n];
            for j in 0..n {
                for i in 0..n {
                    let mut s = 0.0;
                    for l in 0..n {
                        s += x[i + l * n] * y[l + j * n];
                    }
                    c[i + j * n] = s;
                }
            }
            c
        };
        let want = mul(&hr, &mul(&a0, &hr));

        let mut h = a0;
        dgehrd(&mut h, n, None);
        for idx in 0..n * n {
            assert!(
                (h[idx] - want[idx]).abs() < 1e-13,
                "H[{idx}] = {} but Hr A Hr gives {}",
                h[idx],
                want[idx]
            );
        }
        // The published result for this matrix, as a regression pin. Note that
        // H(2,1) -- index 5 -- is the subdiagonal and is legitimately nonzero;
        // only entries *below* the subdiagonal must vanish.
        assert!((h[1] - -2.8284271247461903).abs() < 1e-13, "H(1,0) = {}", h[1]);
        assert_eq!(h[2], 0.0, "H(2,0) must be exactly zero");
        assert!((h[5] - -2.25).abs() < 1e-13, "H(2,1) = {}", h[5]);
    }

    /// `dhseqr` with `COMPZ = 'V'`: the accumulated Schur vectors must satisfy
    /// `A = Z T Z^T` with `Z` orthogonal and `T` quasi-triangular, and must leave
    /// `T` and the eigenvalues bit-for-bit identical to the 4-argument form.
    #[test]
    fn dhseqr_with_schur_vectors_gives_the_schur_decomposition() {
        for n in [1usize, 2, 3, 5, 8, 12] {
            let orig = sample(n);
            let mut h = orig.clone();
            let mut z = vec![0.0; n * n];
            dgehrd(&mut h, n, Some(&mut z));
            let mut wr = vec![0.0; n];
            let mut wi = vec![0.0; n];
            dhseqr_with_schur_vectors(&mut h, &mut wr, &mut wi, &mut z, n).unwrap();

            // Z is orthogonal.
            for j in 0..n {
                for k in 0..n {
                    let s: f64 = (0..n).map(|i| z[i + j * n] * z[i + k * n]).sum();
                    let want = if j == k { 1.0 } else { 0.0 };
                    assert!((s - want).abs() < 1e-12, "n={n}: Z^T Z ({j},{k}) = {s}");
                }
            }

            // A Z = Z T, which is A = Z T Z^T without forming a transpose.
            let az = mat_mul(&orig, &z, n);
            let zt = mat_mul(&z, &h, n);
            let scale = norm_max(&orig).max(1.0);
            for idx in 0..n * n {
                assert!(
                    (az[idx] - zt[idx]).abs() / scale < 1e-12,
                    "n={n}: A Z and Z T differ at {idx}: {} vs {}",
                    az[idx],
                    zt[idx]
                );
            }

            // T is quasi-triangular: nothing below the first subdiagonal.
            for j in 0..n {
                for i in (j + 2)..n {
                    assert!(h[i + j * n].abs() < 1e-13, "n={n}: T({i},{j}) is below the subdiagonal");
                }
            }

            // The eigenvalues are unchanged by recording the Schur vectors.
            let mut h2 = orig.clone();
            dgehrd(&mut h2, n, None);
            let mut wr2 = vec![0.0; n];
            let mut wi2 = vec![0.0; n];
            dhseqr(&mut h2, &mut wr2, &mut wi2, n).unwrap();
            for i in 0..n {
                assert!((wr[i] - wr2[i]).abs() < 1e-13 && (wi[i] - wi2[i]).abs() < 1e-13);
            }
        }
    }

    /// Right eigenvectors of a general real matrix, by residual.
    ///
    /// An eigenvector has no unique value -- it is defined up to scaling, and its
    /// sign and phase are arbitrary -- so comparing against the reference's vectors
    /// would be meaningless. What *is* well defined is the defining equation, so
    /// this checks `A v = lambda v` for each pair, at the tolerance the reduction and
    /// the back-substitution accumulate.
    ///
    /// The eigenvalues are deliberately recomputed here with an independent
    /// `dgehrd` + `dhseqr` rather than taken from `dgeev`, so the pairing of
    /// vectors to eigenvalues is not assumed. `dgebal` leaves the `sample`
    /// family untouched (its dominant diagonal already isolates nothing and
    /// needs no scaling), so the two deflation orders coincide and index `j` of
    /// the returned vectors is still the one paired with index `j` here.
    #[test]
    fn dgeev_eigenvectors_satisfy_their_defining_equation() {
        for n in [1usize, 2, 3, 5, 8, 12] {
            // `sample` here returns column-major data; wrap it so the residual can
            // be formed with element access.
            let data = sample(n);
            let a = data.clone();

            let mut av = data.clone();
            // The working copy must actually hold the matrix. An earlier version
            // left it zeroed, so `dgeev` decomposed the zero matrix, every
            // eigenvalue was zero, and the residual was trivially satisfied --
            // the test passed while testing nothing.
            let mut v_real = vec![0.0; n * n];
            let mut v_imag = vec![0.0; n * n];
            let (wr_vec, wi_vec) = dgeev(&mut v_real, &mut v_imag, &mut av, n).unwrap();

            // Use the eigenvalues `dgeev` itself returns, and its own eigenvector
            // layout. Recovering them from an independent `dgehrd` + `dhseqr`
            // would give a different *order*, because `dgeev` balances first --
            // pairing an independent ordering against these vectors would
            // mismatch columns and fail for the wrong reason.
            let (wr, wi) = (&wr_vec, &wi_vec);

            let scale = a.iter().fold(0.0f64, |m, x| m.max(x.abs())).max(1.0);
            let mut j = 0usize;
            while j < n {
                if wi[j] != 0.0 {
                    // Complex pair: v = (v_real[:,j] + i v_imag[:,j]), lambda = wr + i wi.
                    let (lr, li) = (wr[j], wi[j]);
                    let mut worst = 0.0f64;
                    for i in 0..n {
                        // (A v)_i for the complex vector v.
                        let mut sr = 0.0;
                        let mut si = 0.0;
                        for k in 0..n {
                            let ar = a[i + k * n];
                            let vr = v_real[k + j * n];
                            let vi = v_imag[k + j * n];
                            sr += ar * vr;
                            si += ar * vi;
                        }
                        // A v - lambda v
                        let dr = sr - (lr * v_real[i + j * n] - li * v_imag[i + j * n]);
                        let di = si - (lr * v_imag[i + j * n] + li * v_real[i + j * n]);
                        worst = worst.max(dr.abs().max(di.abs()));
                    }
                    assert!(
                        worst / scale < 1e-10,
                        "n={n}: the complex eigenvector for lambda = {lr}+{li}i has residual {worst} (scale {scale})"
                    );
                    j += 2;
                } else {
                    let l = wr[j];
                    let mut worst = 0.0f64;
                    for i in 0..n {
                        let mut s = 0.0;
                        for k in 0..n {
                            s += a[i + k * n] * v_real[k + j * n];
                        }
                        worst = worst.max((s - l * v_real[i + j * n]).abs());
                    }
                    assert!(
                        worst / scale < 1e-10,
                        "n={n}: the real eigenvector for lambda = {l} has residual {worst} (scale {scale})"
                    );
                    j += 1;
                }
            }
        }
    }

    /// The same defining equation as the acceptance test, over matrices the
    /// `sample` family does not reach: defective, structured and pseudo-random,
    /// including cases that exercise `dlaln2`'s small-pivot perturbation.
    ///
    /// Unlike the acceptance test this uses the eigenvalues `dgeev` returns, so it
    /// also checks that the vectors are paired with the eigenvalues reported for
    /// them.
    #[test]
    fn dgeev_residuals_on_a_battery_of_matrices() {
        let mut state: u64 = 0x243f_6a88_85a3_08d3;
        let mut next = move || {
            state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
            ((state >> 11) as f64) / ((1u64 << 53) as f64) - 0.5
        };

        let mut cases: Vec<(usize, Vec<f64>)> = Vec::new();

        // A defective Jordan block: one eigenvector for a quadruple eigenvalue.
        let mut jordan = vec![0.0; 16];
        for i in 0..4 {
            jordan[i + i * 4] = 2.0;
            if i < 3 {
                jordan[i + (i + 1) * 4] = 1.0;
            }
        }
        cases.push((4, jordan));

        // A pure rotation: the eigenvalues are 1 +- 2i.
        cases.push((2, vec![1.0, -2.0, 2.0, 1.0]));

        // Three decoupled rotations: three conjugate pairs.
        let mut blockrot = vec![0.0; 36];
        for b in 0..3usize {
            let o = 2 * b;
            blockrot[o + o * 6] = 1.0 + b as f64;
            blockrot[(o + 1) + o * 6] = -2.0;
            blockrot[o + (o + 1) * 6] = 2.0;
            blockrot[(o + 1) + (o + 1) * 6] = 1.0 + b as f64;
        }
        cases.push((6, blockrot));

        // Pseudo-random matrices, deterministic so the test is reproducible.
        for n in [3usize, 5, 9, 16] {
            for _ in 0..4 {
                let mut a = vec![0.0; n * n];
                for x in a.iter_mut() {
                    *x = next() * 20.0;
                }
                cases.push((n, a));
            }
        }

        for (n, data) in cases {
            let a = data.clone();
            let mut av = vec![0.0; n * n];
            let mut v_real = vec![0.0; n * n];
            let mut v_imag = vec![0.0; n * n];
            let (wr, wi) = dgeev(&mut av, &mut v_real, &mut v_imag, n).unwrap();

            let scale = a.iter().fold(0.0f64, |m, x| m.max(x.abs())).max(1.0);
            let mut j = 0usize;
            while j < n {
                let (lr, li) = (wr[j], wi[j]);
                let complex = li != 0.0;
                let mut worst = 0.0f64;
                for i in 0..n {
                    let mut sr = 0.0;
                    let mut si = 0.0;
                    for k in 0..n {
                        sr += a[i + k * n] * v_real[k + j * n];
                        si += a[i + k * n] * v_imag[k + j * n];
                    }
                    let dr = sr - (lr * v_real[i + j * n] - li * v_imag[i + j * n]);
                    let di = si - (lr * v_imag[i + j * n] + li * v_real[i + j * n]);
                    worst = worst.max(dr.abs().max(di.abs()));
                }
                assert!(worst / scale < 1e-10, "n={n}: residual {worst} (scale {scale}) for lambda = {lr}{li:+}i");
                j += if complex { 2 } else { 1 };
            }
        }
    }

    /// The left eigenvectors must satisfy `u^T A = lambda u^T`, i.e. for the
    /// complex vector `u = u_real + i u_imag`,
    /// `A^H u = conj(lambda) u`. Nothing in the existing tests checks this, so it
    /// is the acceptance test for [`dgeev_lr`]'s `SIDE = 'L'` branch.
    #[test]
    fn dgeev_lr_left_eigenvectors_satisfy_their_defining_equation() {
        // A deterministic generator, as in the battery test below.
        let mut state: u64 = 0x9e37_79b9_7f4a_7c15;
        let mut next = move || {
            state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
            ((state >> 11) as f64) / ((1u64 << 53) as f64) - 0.5
        };

        let mut cases: Vec<(usize, Vec<f64>)> = Vec::new();
        for n in [1usize, 2, 3, 5, 8] {
            cases.push((n, sample(n)));
        }
        // A pure rotation: the eigenvalues are 1 +- 2i, so the complex-pair
        // branch of the left back-substitution is exercised.
        cases.push((2, vec![1.0, -2.0, 2.0, 1.0]));
        // Three decoupled rotations: three conjugate pairs.
        let mut blockrot = vec![0.0; 36];
        for b in 0..3usize {
            let o = 2 * b;
            blockrot[o + o * 6] = 1.0 + b as f64;
            blockrot[(o + 1) + o * 6] = -2.0;
            blockrot[o + (o + 1) * 6] = 2.0;
            blockrot[(o + 1) + (o + 1) * 6] = 1.0 + b as f64;
        }
        cases.push((6, blockrot));
        for n in [3usize, 7] {
            for _ in 0..4 {
                let mut a = vec![0.0; n * n];
                for x in a.iter_mut() {
                    *x = next() * 20.0;
                }
                cases.push((n, a));
            }
        }

        for (n, data) in cases {
            let a = data.clone();

            // The working copy must hold the same matrix; `a` is kept intact so
            // the residual can be formed against the original.
            let mut av = data.clone();
            let mut u_real = vec![0.0; n * n];
            let mut u_imag = vec![0.0; n * n];
            let mut v_real = vec![0.0; n * n];
            let mut v_imag = vec![0.0; n * n];
            let (wr, wi) = dgeev_lr(&mut u_real, &mut u_imag, &mut v_real, &mut v_imag, &mut av, n).unwrap();

            let scale = a.iter().fold(0.0f64, |m, x| m.max(x.abs())).max(1.0);
            let mut j = 0usize;
            while j < n {
                let (lr, li) = (wr[j], wi[j]);
                let complex = li != 0.0;
                // u^H A = lambda u^H  <=>  A^H u = conj(lambda) u.
                let mut worst = 0.0f64;
                for i in 0..n {
                    let mut sr = 0.0;
                    let mut si = 0.0;
                    for k in 0..n {
                        // A^H u, i.e. sum_k conj(A(k,i)) u(k).
                        sr += a[k + i * n] * u_real[k + j * n];
                        si += a[k + i * n] * u_imag[k + j * n];
                    }
                    // conj(lambda) * u(i) = (lr - i*li)(ur + i*ui)
                    let ur = u_real[i + j * n];
                    let ui = u_imag[i + j * n];
                    let dr = sr - (lr * ur + li * ui);
                    let di = si - (lr * ui - li * ur);
                    worst = worst.max(dr.abs().max(di.abs()));
                }
                assert!(
                    worst / scale < 1e-12,
                    "n={n}: left residual {worst} (scale {scale}) for lambda = {lr}{li:+}i"
                );
                j += if complex { 2 } else { 1 };
            }

            // The right eigenvectors must be unchanged by the extra left work.
            let mut vr = vec![0.0; n * n];
            let mut vi = vec![0.0; n * n];
            let mut aw = data.clone();
            let (wr2, wi2) = dgeev(&mut vr, &mut vi, &mut aw, n).unwrap();
            for i in 0..n {
                assert!(wr[i] == wr2[i] && wi[i] == wi2[i], "n={n}: eigenvalue {i} differs");
                for r in 0..n {
                    assert!(
                        (v_real[r + i * n] - vr[r + i * n]).abs() < 1e-13 && (v_imag[r + i * n] - vi[r + i * n]).abs() < 1e-13,
                        "n={n}: right eigenvector ({r},{i}) differs"
                    );
                }
            }
        }
    }

    /// Non-convergence must be reported, not dropped on the floor.
    ///
    /// A matrix of NaN makes the QR iteration's convergence test
    /// (`|h(i+1,i)| <= threshold`) permanently false, so `dlahqr` exhausts its
    /// iteration limit and returns `INFO > 0`. Before that `INFO` was plumbed
    /// through it was a `debug_assert`, so a release build returned partly
    /// converged eigenvalues as if they had converged. To check this test can
    /// actually fail, put the `debug_assert` back and build in release: the
    /// assertions below then see `Ok`.
    ///
    /// `n >= 4` on purpose: at `n == 2` the matrix is a single 2x2 block, which
    /// `dlanv2` reduces algebraically without iterating, so there is no iteration
    /// limit to exhaust -- reference `dhseqr` reports success there too.
    #[test]
    fn non_convergence_is_reported_rather_than_returned_as_a_result() {
        for n in [4usize, 8] {
            let mut wr = vec![0.0; n];
            let mut wi = vec![0.0; n];

            let mut t = vec![f64::NAN; n * n];
            assert_eq!(
                dhseqr(&mut t, &mut wr, &mut wi, n),
                Err(LapackError::Convergence),
                "dhseqr n={n} must report non-convergence"
            );

            let mut z = vec![0.0; n * n];
            for i in 0..n {
                z[i + i * n] = 1.0;
            }
            let mut t2 = vec![f64::NAN; n * n];
            assert_eq!(
                dhseqr_with_schur_vectors(&mut t2, &mut wr, &mut wi, &mut z, n),
                Err(LapackError::Convergence),
                "dhseqr_with_schur_vectors n={n} must report non-convergence"
            );

            let mut v_real = vec![0.0; n * n];
            let mut v_imag = vec![0.0; n * n];
            let mut a = vec![f64::NAN; n * n];
            assert_eq!(
                dgeev(&mut v_real, &mut v_imag, &mut a, n).map(|_| ()),
                Err(LapackError::Convergence),
                "dgeev n={n} must report non-convergence"
            );

            let mut u_real = vec![0.0; n * n];
            let mut u_imag = vec![0.0; n * n];
            let mut v_real2 = vec![0.0; n * n];
            let mut v_imag2 = vec![0.0; n * n];
            let mut a2 = vec![f64::NAN; n * n];
            assert_eq!(
                dgeev_lr(&mut u_real, &mut u_imag, &mut v_real2, &mut v_imag2, &mut a2, n).map(|_| ()),
                Err(LapackError::Convergence),
                "dgeev_lr n={n} must report non-convergence"
            );
        }
    }
    /// `dgeev_lr` on matrices whose entries span the representable range: balancing,
    /// the Hessenberg reduction and the eigenvector back-substitution all have to
    /// rescale or the products overflow. Uniform extremes and the growth cases (a tiny
    /// diagonal with large off-diagonals, where the triangular solve passes `bignum`)
    /// are both covered, and both sides are checked. An absolute tolerance says
    /// nothing at 1e300, so the residual is the reference's eps-normalized ratio,
    /// `||A v - lambda v|| / (n ||A|| eps) <= 20`, and the same for the left side.
    #[test]
    fn dgeev_handles_matrices_at_the_edges_of_the_range() {
        let n = 4usize;
        let eps = f64::EPSILON;
        let cases = [
            (1e300f64, 1e300f64),
            (1e-300, 1e-300),
            (1e250, 1e250),
            (1e-250, 1e-250),
            (1.0, 1.0),
            (1.0, 1e200),
            (1e-200, 1e200),
            (1e-100, 1e150),
        ];
        for (d, o) in cases {
            let mut a0 = vec![0.0; n * n];
            for i in 0..n {
                a0[i + i * n] = d * (i as f64 + 1.0);
                for j in (i + 1)..n {
                    a0[i + j * n] = o * 0.3 / (j - i) as f64;
                }
            }
            let mut a = a0.clone();
            let mut u_real = vec![0.0; n * n];
            let mut u_imag = vec![0.0; n * n];
            let mut v_real = vec![0.0; n * n];
            let mut v_imag = vec![0.0; n * n];
            let (wr, wi) = dgeev_lr(&mut u_real, &mut u_imag, &mut v_real, &mut v_imag, &mut a, n).expect("dgeev_lr");

            let anorm = a0.iter().fold(0.0f64, |m, x| m.max(x.abs()));
            let bound = n as f64 * anorm * eps;
            for j in 0..n {
                assert!(wr[j].is_finite() && wi[j].is_finite(), "d = {d}, o = {o}: eigenvalue {j} is not finite");
            }
            let mut j = 0usize;
            while j < n {
                let (lr, li) = (wr[j], wi[j]);
                let complex = li != 0.0;
                for (v_re, v_im, left) in [(&v_real, &v_imag, false), (&u_real, &u_imag, true)] {
                    let mut worst = 0.0f64;
                    for i in 0..n {
                        let mut row = 0.0;
                        let mut row_i = 0.0;
                        for k in 0..n {
                            // left eigenvectors satisfy u^T A = lambda u^T
                            let (a_ik, a_ki) = if left { (a0[k + i * n], 0.0) } else { (a0[i + k * n], 0.0) };
                            row += a_ik * v_re[k + j * n] - a_ki * v_im[k + j * n];
                            row_i += a_ik * v_im[k + j * n] + a_ki * v_re[k + j * n];
                        }
                        let dr = row - (lr * v_re[i + j * n] - li * v_im[i + j * n]);
                        let di = row_i - (lr * v_im[i + j * n] + li * v_re[i + j * n]);
                        worst = worst.max(dr.abs().max(di.abs()));
                    }
                    let ratio = worst / bound;
                    assert!(ratio <= 20.0, "d = {d}, o = {o}, j = {j}, left = {left}: ratio = {ratio:e}");
                }
                j += if complex { 2 } else { 1 };
            }
        }
    }
}