// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2026 Rasmus Munk Larsen <rmlarsen@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0

#include "main.h"
#include "fp_control.h"
#include "tridiag_test_matrices.h"
#include <cfenv>
#include <limits>
#include <Eigen/Eigenvalues>

template <typename Packet>
EIGEN_DONT_INLINE Packet packet_product_runtime(const Packet& lhs, const Packet& rhs) {
  return internal::pmul(lhs, rhs);
}

// Some SIMD units flush subnormal operands and results to zero regardless of any library-side care
// (ARMv7 NEON is hard-wired FTZ while its scalar VFP unit honors subnormals). Sections that feed
// genuinely subnormal data through vectorized kernels cannot even observe their inputs on such
// hardware, so they are gated on this runtime probe of the packet path. The volatile load prevents
// constant folding, while the EIGEN_DONT_INLINE wrapper prevents Clang from scalarizing the packet
// product and applying scalar IEEE semantics.
template <typename RealScalar>
bool packet_path_flushes_subnormals() {
  typedef typename internal::packet_traits<RealScalar>::type Packet;
  volatile RealScalar vtiny = (std::numeric_limits<RealScalar>::min)();
  const RealScalar tiny = vtiny;
  const RealScalar half_tiny =
      internal::pfirst(packet_product_runtime(internal::pset1<Packet>(tiny), internal::pset1<Packet>(RealScalar(0.5))));
  return numext::is_exactly_zero(half_tiny);
}

// Dense (n x n) symmetric tridiagonal matrix from its diagonal d and sub-diagonal e, for
// residual / reconstruction checks against the solver output.
template <typename DiagType, typename SubType>
Matrix<typename DiagType::Scalar, Dynamic, Dynamic> dense_symmetric_tridiag(const DiagType& d, const SubType& e) {
  typedef typename DiagType::Scalar RealScalar;
  const Index n = d.size();
  Matrix<RealScalar, Dynamic, Dynamic> T = Matrix<RealScalar, Dynamic, Dynamic>::Zero(n, n);
  T.diagonal() = d;
  if (n > 1) {
    T.diagonal(-1) = e;
    T.diagonal(1) = e;
  }
  return T;
}

// Test TridiagonalEigenSolver (SIMD Sturm-sequence spectral bisection) on the full
// structured-tridiagonal catalog: compare against the implicit-QR path of
// SelfAdjointEigenSolver::computeFromTridiagonal(), exercise index/value range subset
// selection, and verify absolute accuracy against matrices with known spectra.
template <typename RealScalar>
void tridiagonal_eigensolver_bisection() {
  typedef Matrix<RealScalar, Dynamic, Dynamic> MatrixType;
  typedef Matrix<RealScalar, Dynamic, 1> VectorType;
  const RealScalar eps = NumTraits<RealScalar>::epsilon();
  const RealScalar tiny = (std::numeric_limits<RealScalar>::min)();

  // (a) Compare bisection against the QR path over every structured matrix, and exercise
  // index- and value-range subset selection.
  test::for_all_symmetric_tridiag_test_matrices<RealScalar>([&](const VectorType& diag, const VectorType& offdiag) {
    const Index n = diag.size();

    SelfAdjointEigenSolver<MatrixType> qr;
    qr.computeFromTridiagonal(diag, offdiag, EigenvaluesOnly);
    if (qr.info() != Success) return;  // skip pathological inputs the QR path itself rejects
    const VectorType w = qr.eigenvalues();

    TridiagonalEigenSolver<RealScalar> bz;
    bz.computeEigenvalues(diag, offdiag);
    VERIFY_IS_EQUAL(bz.info(), Success);
    VERIFY_IS_EQUAL(bz.eigenvalues().size(), n);
    for (Index i = 1; i < n; ++i) VERIFY(bz.eigenvalues()(i) >= bz.eigenvalues()(i - 1));

    // Bisection has an absolute accuracy floor of ~pivmin (the smallest safe pivot),
    // so the tolerance carries an absolute term in addition to the eps*||T|| term.
    const RealScalar radius = w.cwiseAbs().maxCoeff();
    const RealScalar tol = RealScalar(64) * RealScalar(n) * (eps * radius + tiny);
    VERIFY((bz.eigenvalues() - w).cwiseAbs().maxCoeff() <= tol);

    // Index-range subsets: the k smallest and the k largest.
    if (n >= 4) {
      const Index k = n / 3;
      TridiagonalEigenSolver<RealScalar> lo, hi;
      lo.computeEigenvalues(diag, offdiag, EigenvalueRange::indices(0, k));
      VERIFY_IS_EQUAL(lo.eigenvalues().size(), k);
      VERIFY((lo.eigenvalues() - w.head(k)).cwiseAbs().maxCoeff() <= tol);
      hi.computeEigenvalues(diag, offdiag, EigenvalueRange::indices(n - k, n));
      VERIFY_IS_EQUAL(hi.eigenvalues().size(), k);
      VERIFY((hi.eigenvalues() - w.tail(k)).cwiseAbs().maxCoeff() <= tol);
    }

    // Value-range subset [vl, vu) about the spectrum center. The count is compared strictly
    // only when no eigenvalue sits within tol of a boundary (otherwise it is legitimately
    // ambiguous which side the boundary eigenvalue falls on).
    if (n >= 2 && radius > tiny) {
      const RealScalar vl = -RealScalar(0.3) * radius, vu = RealScalar(0.3) * radius;
      TridiagonalEigenSolver<RealScalar> bv;
      bv.computeEigenvalues(diag, offdiag,
                            EigenvalueRange::values(static_cast<long double>(vl), static_cast<long double>(vu)));
      for (Index i = 0; i < bv.eigenvalues().size(); ++i)
        VERIFY(bv.eigenvalues()(i) >= vl - tol && bv.eigenvalues()(i) < vu + tol);
      const bool clean = (w.array() - vl).abs().minCoeff() > tol && (w.array() - vu).abs().minCoeff() > tol;
      if (clean) {
        const Index expected = (w.array() >= vl && w.array() < vu).count();
        VERIFY_IS_EQUAL(bv.eigenvalues().size(), expected);
        if (expected > 0) {
          Index first = 0;
          while (first < n && w(first) < vl) ++first;
          VERIFY((bv.eigenvalues() - w.segment(first, expected)).cwiseAbs().maxCoeff() <= tol);
        }
      }
    }
  });

  // (a2) Half-open value-range convention, fixed-size inputs, and input validation.
  {
    // values(vl, vu) selects the eigenvalues in the lower-closed, upper-open interval [vl, vu) (like
    // indices()): an eigenvalue exactly equal to vl is included, one exactly equal to vu is not. Use a
    // diagonal matrix (zero off-diagonal) with the exact integer spectrum {1, 2, 3} so the endpoints
    // land precisely on eigenvalues.
    VectorType d(3), e(2);
    d << RealScalar(1), RealScalar(2), RealScalar(3);
    e.setZero();
    TridiagonalEigenSolver<RealScalar> bmid;
    bmid.computeEigenvalues(d, e, EigenvalueRange::values(1.0L, 3.0L));
    VERIFY_IS_EQUAL(bmid.eigenvalues().size(), Index(2));  // [1, 3) = {1, 2}: lower end closed, upper end open
    VERIFY_IS_APPROX(bmid.eigenvalues()(0), RealScalar(1));
    VERIFY_IS_APPROX(bmid.eigenvalues()(1), RealScalar(2));
    TridiagonalEigenSolver<RealScalar> bopen;
    bopen.computeEigenvalues(d, e, EigenvalueRange::values(0.0L, 1.0L));
    VERIFY_IS_EQUAL(bopen.eigenvalues().size(), Index(0));  // [0, 1) excludes the eigenvalue at the open end 1
    TridiagonalEigenSolver<RealScalar> blow;
    blow.computeEigenvalues(d, e, EigenvalueRange::values(1.0L, 2.0L));
    VERIFY_IS_EQUAL(blow.eigenvalues().size(), Index(1));  // [1, 2) = {1}: includes the eigenvalue at the closed end 1
    VERIFY_IS_APPROX(blow.eigenvalues()(0), RealScalar(1));

    // The half-open endpoint semantics must survive rounding in the Sturm recurrence, not just the
    // exact-count diagonal case above: [[-3, 2], [2, -3]] has spectrum {-5, -1} exactly, but the
    // recurrence rounds through 2/3-type values, so counts taken exactly at an endpoint eigenvalue
    // are a coin flip (float used to drop -5 entirely, double used to keep -1). The endpoint
    // resolution against the converged eigenvalues must return exactly {-5}.
    VectorType dcross(2), ecross(1);
    dcross << RealScalar(-3), RealScalar(-3);
    ecross << RealScalar(2);
    TridiagonalEigenSolver<RealScalar> bcross;
    bcross.computeEigenvalues(dcross, ecross, EigenvalueRange::values(-5.0L, -1.0L));
    VERIFY_IS_EQUAL(bcross.eigenvalues().size(), Index(1));
    VERIFY_IS_APPROX(bcross.eigenvalues()(0), RealScalar(-5));

    // Infinite endpoints mean unbounded on that side; the endpoint-tolerance arithmetic must not
    // poison the result (inf - inf = NaN used to return an empty spectrum). A finite long double
    // endpoint that narrows to infinity in RealScalar (1e300 for float) behaves the same way.
    TridiagonalEigenSolver<RealScalar> binf;
    binf.computeEigenvalues(d, e,
                            EigenvalueRange::values(-std::numeric_limits<long double>::infinity(),
                                                    std::numeric_limits<long double>::infinity()));
    VERIFY_IS_EQUAL(binf.eigenvalues().size(), Index(3));
    binf.computeEigenvalues(d, e, EigenvalueRange::values(1.5L, std::numeric_limits<long double>::infinity()));
    VERIFY_IS_EQUAL(binf.eigenvalues().size(), Index(2));
    VERIFY_IS_APPROX(binf.eigenvalues()(0), RealScalar(2));
    binf.computeEigenvalues(d, e, EigenvalueRange::values(1.5L, 1e300L));
    VERIFY_IS_EQUAL(binf.eigenvalues().size(), Index(2));

    // Fixed-size input vectors work, including subset selection (the solver's own storage is
    // dynamic, so a subset shorter than the input is not a problem).
    Matrix<RealScalar, 3, 1> fd;
    fd << RealScalar(1), RealScalar(2), RealScalar(3);
    Matrix<RealScalar, 2, 1> fe;
    fe.setZero();
    TridiagonalEigenSolver<RealScalar> fz;
    fz.computeEigenvalues(fd, fe);
    VERIFY_IS_EQUAL(fz.eigenvalues().size(), Index(3));
    fz.computeEigenvalues(fd, fe, EigenvalueRange::indices(0, 1));
    VERIFY_IS_EQUAL(fz.eigenvalues().size(), Index(1));
    VERIFY_IS_APPROX(fz.eigenvalues()(0), RealScalar(1));

    // Querying an uninitialized solver is a usage error.
    TridiagonalEigenSolver<RealScalar> uninit;
    VERIFY_RAISES_ASSERT(uninit.eigenvalues());
    VERIFY_RAISES_ASSERT(uninit.info());

    // A non-finite input (NaN or Inf, diagonal or sub-diagonal) is rejected up front and
    // reported via info().
    VectorType dbad = d;
    dbad(1) = std::numeric_limits<RealScalar>::quiet_NaN();
    TridiagonalEigenSolver<RealScalar> bnan;
    bnan.computeEigenvalues(dbad, e);
    VERIFY_IS_EQUAL(bnan.info(), NoConvergence);
    dbad(1) = std::numeric_limits<RealScalar>::infinity();
    bnan.computeEigenvalues(dbad, e);
    VERIFY_IS_EQUAL(bnan.info(), NoConvergence);
    VectorType ebad = e;
    ebad(0) = std::numeric_limits<RealScalar>::quiet_NaN();
    bnan.computeEigenvalues(d, ebad);
    VERIFY_IS_EQUAL(bnan.info(), NoConvergence);
  }

  // (b) Absolute-accuracy checks against matrices with known spectra (independent of QR).
  // A subnormal-magnitude matrix exercises the normalization: dividing entries by the (subnormal) max
  // magnitude keeps the scaled matrix O(1), whereas multiplying by 1/scale would overflow to infinity
  // and return wrong eigenvalues. Compare against the same matrix in a normal range, scaled back down.
  // Skipped on flush-to-zero packet hardware, where the subnormal inputs read as zero (see the probe).
  if (!packet_path_flushes_subnormals<RealScalar>()) {
    const Index nn = 7;
    VectorType bd(nn), be(nn - 1);
    test::tridiag_1_2_1(bd, be);
    const RealScalar sub_scale = (std::numeric_limits<RealScalar>::min)() / RealScalar(64);  // subnormal
    TridiagonalEigenSolver<RealScalar> bref, bsub;
    bref.computeEigenvalues(bd, be);
    bsub.computeEigenvalues(VectorType(bd * sub_scale), VectorType(be * sub_scale));
    VERIFY_IS_EQUAL(bsub.info(), Success);
    // The subnormal eigenvalues, scaled back to O(1), agree with the reference to the subnormal
    // granularity at this magnitude (denorm_min / sub_scale) plus the usual bisection error.
    const RealScalar radius = bref.eigenvalues().cwiseAbs().maxCoeff();
    const RealScalar gran = (std::numeric_limits<RealScalar>::denorm_min)() / sub_scale;
    VERIFY((bsub.eigenvalues() / sub_scale - bref.eigenvalues()).cwiseAbs().maxCoeff() <=
           RealScalar(16) * (RealScalar(nn) * eps + gran) * radius);
  }
  const double pi = 3.14159265358979323846;
  for (Index n : {7, 16, 33, 64}) {
    VectorType d(n), e(n - 1), wexact(n);

    // 1-2-1 Toeplitz: lambda_k = 2 - 2 cos(k*pi/(n+1)), ascending for k = 1..n.
    test::tridiag_1_2_1(d, e);
    TridiagonalEigenSolver<RealScalar> b121;
    b121.computeEigenvalues(d, e);
    for (Index k = 0; k < n; ++k) wexact(k) = RealScalar(2.0 - 2.0 * std::cos(double(k + 1) * pi / double(n + 1)));
    VERIFY((b121.eigenvalues() - wexact).cwiseAbs().maxCoeff() <= RealScalar(64) * RealScalar(n) * eps * RealScalar(4));

    // Clement / Kac: integer spectrum -(n-1), -(n-3), ..., (n-1).
    test::tridiag_clement(d, e);
    TridiagonalEigenSolver<RealScalar> bcl;
    bcl.computeEigenvalues(d, e);
    for (Index k = 0; k < n; ++k) wexact(k) = RealScalar(-(n - 1) + 2 * k);
    VERIFY((bcl.eigenvalues() - wexact).cwiseAbs().maxCoeff() <= RealScalar(64) * RealScalar(n) * eps * RealScalar(n));

    // Spectra exactly symmetric about 0: signed Wilkinson, Hermite, Legendre, [1,0,1] Toeplitz.
    auto check_symmetric = [&](const VectorType& dd, const VectorType& ee) {
      TridiagonalEigenSolver<RealScalar> bs;
      bs.computeEigenvalues(dd, ee);
      const VectorType s = bs.eigenvalues();
      const RealScalar rad = s.cwiseAbs().maxCoeff();
      const RealScalar tol = RealScalar(64) * RealScalar(n) * (eps * rad + tiny);
      for (Index i = 0; i < n; ++i) VERIFY(numext::abs(s(i) + s(n - 1 - i)) <= tol);
    };
    test::tridiag_wilkinson_signed(d, e);
    check_symmetric(d, e);
    test::tridiag_hermite(d, e);
    check_symmetric(d, e);
    test::tridiag_legendre(d, e);
    check_symmetric(d, e);
    test::tridiag_toeplitz(d, e);  // default a = 1, b = 0
    check_symmetric(d, e);
  }
}

// Test the inverse-iteration eigenvector stage (LAPACK xSTEIN analog): staged
// (computeEigenvalues() followed by computeEigenvectors()), one-call compute(), and the direct
// computeEigenvectors(diag, subdiag, eigenvalues) API.
template <typename RealScalar>
void tridiagonal_eigensolver_eigenvectors() {
  typedef Matrix<RealScalar, Dynamic, Dynamic> MatrixType;
  typedef Matrix<RealScalar, Dynamic, 1> VectorType;
  const RealScalar eps = NumTraits<RealScalar>::epsilon();
  const RealScalar tiny = (std::numeric_limits<RealScalar>::min)();
  const double pi = 3.14159265358979323846;

  // (a) Eigenvectors over the structured catalog: staged and direct APIs.
  test::for_all_symmetric_tridiag_test_matrices<RealScalar>([&](const VectorType& diag, const VectorType& offdiag) {
    const Index n = diag.size();

    // Dense form of T for residual / reconstruction checks.
    const MatrixType T = dense_symmetric_tridiag(diag, offdiag);
    RealScalar scale = T.cwiseAbs().maxCoeff();
    if (!(numext::isfinite)(scale)) return;  // skip non-finite inputs, like the eigenvalue path
    if (numext::is_exactly_zero(scale)) scale = RealScalar(1);

    TridiagonalEigenSolver<RealScalar> es;
    es.computeEigenvalues(diag, offdiag);
    if (es.info() != Success) return;
    es.computeEigenvectors();  // staged inverse-iteration pass
    VERIFY_IS_EQUAL(es.info(), Success);
    const VectorType w = es.eigenvalues();
    const MatrixType V = es.eigenvectors();
    VERIFY_IS_EQUAL(V.rows(), n);
    VERIFY_IS_EQUAL(V.cols(), n);

    const RealScalar tol = RealScalar(128) * RealScalar(n) * (eps * scale + tiny);
    const RealScalar otol = RealScalar(128) * RealScalar(n) * eps;
    // Per-column residual; stableNorm() since the catalog includes near-overflow entries.
    for (Index i = 0; i < n; ++i) VERIFY((T * V.col(i) - w(i) * V.col(i)).stableNorm() <= tol);
    // Orthonormality and full-spectrum reconstruction (the latter in scaled coordinates so the
    // extreme-magnitude catalog entries cannot overflow the products).
    VERIFY((V.transpose() * V - MatrixType::Identity(n, n)).cwiseAbs().maxCoeff() <= otol);
    const RealScalar recon =
        ((V * (w / scale).asDiagonal() * V.transpose()) - (T / scale)).cwiseAbs().maxCoeff() * scale;
    VERIFY(recon <= tol);

    // The one-call compute() reproduces the staged result exactly.
    TridiagonalEigenSolver<RealScalar> onecall;
    onecall.compute(diag, offdiag);
    VERIFY_IS_EQUAL(onecall.info(), Success);
    VERIFY_IS_EQUAL((onecall.eigenvalues() - w).cwiseAbs().maxCoeff(), RealScalar(0));
    VERIFY_IS_EQUAL((onecall.eigenvectors() - V).cwiseAbs().maxCoeff(), RealScalar(0));

    // The direct API on the same eigenvalues reproduces the staged result exactly, unless T splits into
    // disconnected blocks (xSTEBZ's criterion, here with a 2x margin): the staged pass knows its eigenvalues
    // are only good to the bisection tolerance and may re-bisect a small block, while supplied eigenvalues
    // are taken as exact.
    TridiagonalEigenSolver<RealScalar> dir;
    dir.computeEigenvectors(diag, offdiag, w);
    VERIFY_IS_EQUAL(dir.eigenvectors().rows(), n);
    VERIFY_IS_EQUAL(dir.eigenvectors().cols(), n);
    const bool may_split =
        n > 1 && (offdiag.array().abs() <=
                  RealScalar(2) * eps * (diag.head(n - 1).array().abs().sqrt() * diag.tail(n - 1).array().abs().sqrt()))
                     .any();
    if (!may_split) VERIFY_IS_EQUAL((dir.eigenvectors() - V).cwiseAbs().maxCoeff(), RealScalar(0));
    VERIFY_IS_EQUAL(dir.info(), Success);
    VERIFY_IS_EQUAL(dir.eigenvalues(), w);
    const MatrixType& supplied_vectors = dir.eigenvectors();
    for (Index i = 0; i < n; ++i)
      VERIFY((T * supplied_vectors.col(i) - w(i) * supplied_vectors.col(i)).stableNorm() <= tol);
    VERIFY((supplied_vectors.transpose() * supplied_vectors - MatrixType::Identity(n, n)).cwiseAbs().maxCoeff() <=
           otol);

    // Index-subset eigenvectors: the bisection range selects a band, inverse iteration produces just
    // those columns; check residual and that they are orthonormal among themselves.
    if (n >= 4) {
      const Index il = n / 4, iu = n - n / 4;
      TridiagonalEigenSolver<RealScalar> sub;
      sub.compute(diag, offdiag, ComputeEigenvectors, EigenvalueRange::indices(il, iu));
      const MatrixType Vs = sub.eigenvectors();
      const VectorType ws = sub.eigenvalues();
      VERIFY_IS_EQUAL(Vs.rows(), n);
      VERIFY_IS_EQUAL(Vs.cols(), iu - il);
      for (Index k = 0; k < iu - il; ++k) VERIFY((T * Vs.col(k) - ws(k) * Vs.col(k)).stableNorm() <= tol);
      VERIFY((Vs.transpose() * Vs - MatrixType::Identity(iu - il, iu - il)).cwiseAbs().maxCoeff() <= otol);
    }
  });

  // (b) Closed-form 1-2-1 Toeplitz eigenvectors: v_k(j) = sin(j*(n-k)*pi/(n+1)) for the k-th
  // (ascending) eigenvalue. An analytic, solver-independent check of the inverse-iteration vectors.
  for (Index n : {5, 16, 33, 64}) {
    VectorType d(n), e(n - 1);
    test::tridiag_1_2_1(d, e);
    TridiagonalEigenSolver<RealScalar> es(d, e);
    const MatrixType V = es.eigenvectors();
    for (Index k = 0; k < n; ++k) {
      VectorType exact(n);
      for (Index j = 0; j < n; ++j) exact(j) = RealScalar(std::sin(double(j + 1) * double(n - k) * pi / double(n + 1)));
      exact.normalize();
      const RealScalar err = (std::min)((V.col(k) - exact).norm(), (V.col(k) + exact).norm());
      VERIFY(err <= RealScalar(256) * RealScalar(n) * eps);
    }
  }

  // (c) Edge cases: fixed-size inputs, eigenvalues-only requests, subnormal-magnitude matrices,
  // and convergence reporting.
  {
    // Fixed-size input vectors work for the full eigendecomposition and for subsets (the solver's
    // own storage is dynamic).
    Matrix<RealScalar, 4, 1> fdiag = Matrix<RealScalar, 4, 1>::Random();
    Matrix<RealScalar, 3, 1> fsub = Matrix<RealScalar, 3, 1>::Random();
    TridiagonalEigenSolver<RealScalar> fz;
    fz.compute(fdiag, fsub);
    VERIFY_IS_EQUAL(fz.info(), Success);
    VERIFY_IS_EQUAL(fz.eigenvectors().cols(), Index(4));
    {
      const MatrixType Tf = dense_symmetric_tridiag(fdiag, fsub);
      const MatrixType Vf = fz.eigenvectors();
      VERIFY((Vf.transpose() * Vf - MatrixType::Identity(4, 4)).cwiseAbs().maxCoeff() <= RealScalar(64) * eps);
      VERIFY((Tf * Vf - Vf * fz.eigenvalues().asDiagonal()).cwiseAbs().maxCoeff() <=
             RealScalar(64) * eps * Tf.cwiseAbs().maxCoeff());
    }

    // An EigenvaluesOnly compute() does not produce eigenvectors; querying them is a usage error,
    // and the staged computeEigenvectors() supplies them afterwards.
    TridiagonalEigenSolver<RealScalar> staged;
    staged.compute(fdiag, fsub, EigenvaluesOnly);
    VERIFY_IS_EQUAL(staged.info(), Success);
    VERIFY_RAISES_ASSERT(staged.eigenvectors());
    staged.computeEigenvectors();
    VERIFY_IS_EQUAL(staged.info(), Success);
    VERIFY_IS_EQUAL((staged.eigenvectors() - fz.eigenvectors()).cwiseAbs().maxCoeff(), RealScalar(0));

    // A matrix whose entries are all subnormal must still yield genuine unit-norm, orthonormal
    // eigenvectors rather than all-zero columns. Normalization must avoid forming 1/scale when
    // that reciprocal overflows; subnormal blocks require two finite scaling steps. The eigenvalues are
    // taken from the same matrix in a normal magnitude range (then scaled back down), so this exercises
    // the eigenvector normalization in isolation from the eigenvalue solver.
    // Skipped on flush-to-zero packet hardware, where the subnormal inputs read as zero (see the probe).
    if (!packet_path_flushes_subnormals<RealScalar>()) {
      const Index n = 6;
      VectorType base_d(n), base_e(n - 1);
      test::tridiag_1_2_1(base_d, base_e);
      // sub_scale is subnormal and small enough that 1/sub_scale overflows to infinity.
      const RealScalar sub_scale = (std::numeric_limits<RealScalar>::min)() / RealScalar(64);
      const VectorType d = base_d * sub_scale, e = base_e * sub_scale;

      TridiagonalEigenSolver<RealScalar> ref;
      ref.computeEigenvalues(base_d, base_e);
      const VectorType w_sub = ref.eigenvalues() * sub_scale;

      TridiagonalEigenSolver<RealScalar> se;
      se.computeEigenvectors(d, e, w_sub);
      VERIFY_IS_EQUAL(se.info(), Success);
      const MatrixType V = se.eigenvectors();
      for (Index i = 0; i < n; ++i)
        VERIFY(numext::abs(V.col(i).norm() - RealScalar(1)) <= RealScalar(64) * eps);  // unit norm, not all-zero
      VERIFY((V.transpose() * V - MatrixType::Identity(n, n)).cwiseAbs().maxCoeff() <=
             RealScalar(64) * RealScalar(n) * eps);
      // Residual in normal-range coordinates (the eigenvectors are invariant under the uniform scale).
      const MatrixType Tn = dense_symmetric_tridiag(base_d, base_e);
      for (Index i = 0; i < n; ++i)
        VERIFY((Tn * V.col(i) - ref.eigenvalues()(i) * V.col(i)).norm() <= RealScalar(256) * RealScalar(n) * eps);
    }

    // Numerically disconnected blocks at wildly different scales: each 1x1 block's eigenvector is
    // the corresponding unit vector, exactly (components outside the owning block are exact zeros).
    // Processing the blocks with one global perturbation floor (eps * ||T||) would swamp the pivots
    // of the small blocks and return arbitrary mixtures of the first four vectors instead.
    {
      VectorType d(5), e = VectorType::Zero(4);
      d << RealScalar(1), RealScalar(2), RealScalar(3), RealScalar(4), RealScalar(1e16);
      TridiagonalEigenSolver<RealScalar> blocks;
      blocks.computeEigenvectors(d, e, d);  // the diagonal is the exact (ascending) spectrum
      VERIFY_IS_EQUAL(blocks.info(), Success);
      VERIFY_IS_EQUAL((blocks.eigenvectors() - MatrixType::Identity(5, 5)).cwiseAbs().maxCoeff(), RealScalar(0));

      // Any single requested eigenvalue must select its own block. The localization windows are
      // expressed at each block's own scale; a window at the global scale (~ eps * 1e16 here) spans
      // all four small blocks at once and used to hand w = (3) to the first block.
      for (Index j = 0; j < 5; ++j) {
        VectorType wone(1);
        wone << d(j);
        TridiagonalEigenSolver<RealScalar> one;
        one.computeEigenvectors(d, e, wone);
        VERIFY_IS_EQUAL(one.info(), Success);
        VERIFY_IS_EQUAL((one.eigenvectors().col(0) - MatrixType::Identity(5, 5).col(j)).cwiseAbs().maxCoeff(),
                        RealScalar(0));
      }
    }

    // Small-block shifts need local accuracy even above the old global clustering cutoff.
    for (RealScalar c : {RealScalar(0.125), RealScalar(2), RealScalar(3.9), RealScalar(4), RealScalar(4.1),
                         RealScalar(5), RealScalar(10), RealScalar(100), RealScalar(1000)}) {
      VectorType d(3), e(2);
      d << RealScalar(1), RealScalar(0), RealScalar(0);
      e << RealScalar(0), c * eps;
      for (const auto& range : {EigenvalueRange::all(), EigenvalueRange::indices(0, 2), EigenvalueRange::indices(1, 2),
                                EigenvalueRange::values(-0.5L, 0.5L)}) {
        TridiagonalEigenSolver<RealScalar> blocks;
        blocks.computeEigenvalues(d, e, range);
        const VectorType values = blocks.eigenvalues();
        blocks.computeEigenvectors();
        VERIFY_IS_EQUAL(blocks.info(), Success);
        VERIFY_IS_EQUAL(blocks.eigenvalues(), values);
        const MatrixType V = blocks.eigenvectors();
        VERIFY(V.allFinite());
        VERIFY((V.transpose() * V - MatrixType::Identity(V.cols(), V.cols())).cwiseAbs().maxCoeff() <=
               RealScalar(32) * eps);
        VERIFY((dense_symmetric_tridiag(d, e) * V - V * values.asDiagonal()).cwiseAbs().maxCoeff() <=
               RealScalar(32) * eps);
        // Independent local eigenvectors are (1,+/-1)/sqrt(2), irrespective of global shift error.
        for (Index j = 0; j < V.cols(); ++j) {
          if (numext::abs(values(j)) < RealScalar(0.5))
            VERIFY(numext::abs(numext::abs(V(1, j)) - numext::abs(V(2, j))) <= RealScalar(32) * eps);
        }
      }
    }

    // Accurate supplied shifts must not inherit the unrelated block's scale, even after solver reuse.
    {
      const Index nb = 16;
      VectorType d = VectorType::Constant(nb + 1, RealScalar(2));
      VectorType e = VectorType::Ones(nb), w(nb);
      d(0) = numext::ldexp(RealScalar(1), 60);
      e(0) = RealScalar(0);
      for (Index k = 0; k < nb; ++k) {
        const RealScalar angle = RealScalar(EIGEN_PI) * RealScalar(k + 1) / RealScalar(2 * (nb + 1));
        w(k) = RealScalar(4) * numext::abs2(numext::sin(angle));
      }
      TridiagonalEigenSolver<RealScalar> supplied, local;
      supplied.computeEigenvalues(d, e);
      supplied.computeEigenvectors(d, e, w);
      local.computeEigenvectors(d.tail(nb).eval(), e.tail(nb - 1).eval(), w);
      VERIFY_IS_EQUAL(supplied.info(), Success);
      VERIFY_IS_EQUAL(local.info(), Success);
      VERIFY_IS_EQUAL(supplied.eigenvalues(), w);
      VERIFY_IS_EQUAL(supplied.eigenvectors().bottomRows(nb), local.eigenvectors());
      VERIFY_IS_EQUAL(supplied.eigenvectors().row(0).squaredNorm(), RealScalar(0));
    }

    // A zero block's normalization placeholder must not change clustering under power-of-two scaling.
    {
      VectorType d(4), e(3), w(4);
      d << RealScalar(0), RealScalar(2), RealScalar(2), RealScalar(2);
      e << RealScalar(0), RealScalar(1), RealScalar(1);
      const RealScalar root2 = numext::sqrt(RealScalar(2));
      w << RealScalar(0), RealScalar(2) - root2, RealScalar(2), RealScalar(2) + root2;
      TridiagonalEigenSolver<RealScalar> reference, scaled;
      reference.computeEigenvectors(d, e, w);
      const RealScalar scale = numext::ldexp(RealScalar(1), -60);
      scaled.computeEigenvectors((d * scale).eval(), (e * scale).eval(), (w * scale).eval());
      VERIFY_IS_EQUAL(reference.info(), Success);
      VERIFY_IS_EQUAL(scaled.info(), Success);
      VERIFY_IS_EQUAL((scaled.eigenvectors() - reference.eigenvectors()).cwiseAbs().maxCoeff(), RealScalar(0));
    }

    // A huge disconnected block must not force splitting of a strongly connected small block: a
    // splitting safety floor expressed at the global scale used to break the trailing 2x2 into
    // singletons, returning coordinate vectors with O(1) local residual.
    {
      const RealScalar big = (std::numeric_limits<RealScalar>::max)() / RealScalar(1e10);
      VectorType d(3), e(2), w(3);
      d << big, RealScalar(1), RealScalar(2);
      e << RealScalar(0), RealScalar(0.5);
      // Exact spectrum: big and (3 -+ sqrt(2))/2 from the trailing 2x2 block.
      w << RealScalar((3.0 - std::sqrt(2.0)) / 2.0), RealScalar((3.0 + std::sqrt(2.0)) / 2.0), big;
      TridiagonalEigenSolver<RealScalar> conn;
      conn.computeEigenvectors(d, e, w);
      VERIFY_IS_EQUAL(conn.info(), Success);
      const MatrixType Vc = conn.eigenvectors();
      const MatrixType Tc = dense_symmetric_tridiag(d, e);
      for (Index j = 0; j < 3; ++j)
        VERIFY((Tc * Vc.col(j) - w(j) * Vc.col(j)).stableNorm() <=
               RealScalar(64) * eps * numext::maxi(RealScalar(1), numext::abs(w(j))));
      // The small-block vectors have no support on the big block, exactly.
      VERIFY_IS_EQUAL(numext::abs(Vc(0, 0)), RealScalar(0));
      VERIFY_IS_EQUAL(numext::abs(Vc(0, 1)), RealScalar(0));
      VERIFY(numext::abs(Vc.col(0).dot(Vc.col(1))) <= RealScalar(8) * eps);
    }

    // Non-finite input to the direct eigenvector path is rejected up front and reported via info(),
    // mirroring computeEigenvalues(): NaN data or shifts must not produce NaN vectors with Success.
    {
      VectorType d = VectorType::Ones(3), e = VectorType::Zero(2), w = VectorType::Ones(3);
      TridiagonalEigenSolver<RealScalar> bad;
      VectorType dnan = d;
      dnan(1) = std::numeric_limits<RealScalar>::quiet_NaN();
      bad.computeEigenvectors(dnan, e, w);
      VERIFY_IS_EQUAL(bad.info(), NoConvergence);
      VectorType wnan = w;
      wnan(2) = std::numeric_limits<RealScalar>::quiet_NaN();
      bad.computeEigenvectors(d, e, wnan);
      VERIFY_IS_EQUAL(bad.info(), NoConvergence);
    }

    // Convergence is reported through info() rather than hard-coded to Success: inverse iteration
    // counts the eigenvectors that fail to converge within its step limit (cf. LAPACK xSTEIN), and a
    // well-conditioned spectrum converges fully, so that count is zero and info() is Success.
    {
      const Index n = 32;
      VectorType d = VectorType::Random(n), e = VectorType::Random(n - 1);
      TridiagonalEigenSolver<RealScalar> es2(d, e);
      VERIFY_IS_EQUAL(es2.info(), Success);
    }
  }
}

// Scalar types narrower than float (half, bfloat16) are computed internally in float (see
// TridiagonalEigenSolver::ComputeScalar): Sturm counts and bisection targets stay exact past the
// narrow type's ~1/eps integer ceiling, the start-vector RNG cannot overflow to infinity, and the
// coincident-shift perturbation (~0.08|x| at bfloat16 precision) cannot jump into a neighbouring
// eigenspace. Each block below used to fail in the corresponding way.
template <typename RealScalar>
void tridiagonal_eigensolver_narrow() {
  typedef Matrix<RealScalar, Dynamic, 1> VectorType;
  typedef Matrix<RealScalar, Dynamic, Dynamic> MatrixType;
  const RealScalar eps = NumTraits<RealScalar>::epsilon();

  // (a) 2x2 [2, 1; 1, 2]: exact eigenpairs (1, 3) with vectors (1, -1)/sqrt(2), (1, 1)/sqrt(2).
  // With the RNG overflow this returned an all-NaN eigenvector matrix with info() == Success.
  {
    VectorType d(2), e(1);
    d << RealScalar(2), RealScalar(2);
    e << RealScalar(1);
    TridiagonalEigenSolver<RealScalar> es(d, e);
    VERIFY_IS_EQUAL(es.info(), Success);
    VERIFY(es.eigenvectors().allFinite());
    VERIFY_IS_APPROX(es.eigenvalues()(0), RealScalar(1));
    VERIFY_IS_APPROX(es.eigenvalues()(1), RealScalar(3));
    const MatrixType V = es.eigenvectors();
    VERIFY((V.transpose() * V - MatrixType::Identity(2, 2)).cwiseAbs().maxCoeff() <= RealScalar(8) * eps);
    const MatrixType T = dense_symmetric_tridiag(d, e);
    VERIFY((T * V - V * es.eigenvalues().asDiagonal()).cwiseAbs().maxCoeff() <= RealScalar(16) * eps);
  }

  // (b) Identity of dimension 257 (> 1/eps for bfloat16): every eigenvalue is exactly 1. With the
  // Sturm counts held in the narrow type this returned eigenvalues as large as 5.19.
  {
    const Index n = 257;
    const VectorType d = VectorType::Constant(n, RealScalar(1));
    const VectorType e = VectorType::Zero(n - 1);
    TridiagonalEigenSolver<RealScalar> es;
    es.computeEigenvalues(d, e);
    VERIFY_IS_EQUAL(es.info(), Success);
    VERIFY_IS_EQUAL(es.eigenvalues().size(), n);
    for (Index i = 0; i < n; ++i) VERIFY_IS_EQUAL(es.eigenvalues()(i), RealScalar(1));
  }

  // (c) diag (1, 1, 1 + 2^-4), zero off-diagonal: the eigenvectors for the two eigenvalues at 1
  // must span exactly the first two coordinates. An un-capped shift perturbation at this precision
  // used to select the third eigenspace for the second vector.
  {
    VectorType d(3), e = VectorType::Zero(2), w(2);
    d << RealScalar(1), RealScalar(1), RealScalar(1.0625f);
    w << RealScalar(1), RealScalar(1);
    TridiagonalEigenSolver<RealScalar> es;
    es.computeEigenvectors(d, e, w);
    VERIFY_IS_EQUAL(es.info(), Success);
    const MatrixType V = es.eigenvectors();
    VERIFY_IS_EQUAL(numext::abs(V(2, 0)), RealScalar(0));
    VERIFY_IS_EQUAL(numext::abs(V(2, 1)), RealScalar(0));
    VERIFY((V.transpose() * V - MatrixType::Identity(2, 2)).cwiseAbs().maxCoeff() <= RealScalar(8) * eps);
  }

  // (d) 1-2-1 Toeplitz of dimension 64, staged full solve: the eigenvector pass must shift by the
  // retained unrounded float eigenvalues. Shifts rounded to RealScalar collapse neighbouring
  // eigenvalues and used to return near-parallel columns (Gram error ~0.4 in bfloat16).
  {
    const Index n = 64;
    const VectorType d = VectorType::Constant(n, RealScalar(2.0f));
    const VectorType e = VectorType::Constant(n - 1, RealScalar(1.0f));
    TridiagonalEigenSolver<RealScalar> es(d, e);
    VERIFY_IS_EQUAL(es.info(), Success);
    const MatrixType V = es.eigenvectors();
    VERIFY((V.transpose() * V - MatrixType::Identity(n, n)).cwiseAbs().maxCoeff() <= RealScalar(16) * eps);
  }
}

template <typename = void>
void tridiagonal_eigensolver_power_of_two_scaling() {
  // An input sensitive to normalization rounding; bisection promises normwise accuracy, not correct rounding.
  Vector3f diag;
  diag << numext::bit_cast<float>(numext::uint32_t(0xf055994f)), numext::bit_cast<float>(numext::uint32_t(0x6fb2331a)),
      numext::bit_cast<float>(numext::uint32_t(0x70d02374));
  Vector2f subdiag;
  subdiag << numext::bit_cast<float>(numext::uint32_t(0x70b84b8a)),
      numext::bit_cast<float>(numext::uint32_t(0x6fc495b1));
  Vector3f expected;
  expected << numext::bit_cast<float>(numext::uint32_t(0xf0e811b7)),
      numext::bit_cast<float>(numext::uint32_t(0x708f34f5)), numext::bit_cast<float>(numext::uint32_t(0x70eac055));

  TridiagonalEigenSolver<float> eigenvalues;
  eigenvalues.computeEigenvalues(diag, subdiag);
  VERIFY_IS_EQUAL(eigenvalues.info(), Success);
  VERIFY(eigenvalues.eigenvalues().allFinite());
  const Matrix3d bisectionMatrix = dense_symmetric_tridiag(diag, subdiag).cast<double>();
  const double bisectionNorm = bisectionMatrix.cwiseAbs().rowwise().sum().maxCoeff();
  // Allow O(n*eps*||T||_inf) Sturm-count displacement, the final bracket width and reference rounding.
  // Eight covers the 2.1*n*eps bracket expansion on both sides, with room for the latter two terms.
  const double eigenvalueTolerance = 8.0 * double(diag.size()) * double(NumTraits<float>::epsilon()) * bisectionNorm;
  VERIFY((eigenvalues.eigenvalues().cast<double>() - expected.cast<double>()).cwiseAbs().maxCoeff() <=
         eigenvalueTolerance);

  // Exercise normalization-sensitive inverse iteration with rounded supplied eigenvalues.
  diag << numext::bit_cast<float>(numext::uint32_t(0x44d75b8b)), numext::bit_cast<float>(numext::uint32_t(0x45b9403d)),
      numext::bit_cast<float>(numext::uint32_t(0x4554211f));
  subdiag << numext::bit_cast<float>(numext::uint32_t(0x45609c76)),
      numext::bit_cast<float>(numext::uint32_t(0xc590ed04));
  expected << numext::bit_cast<float>(numext::uint32_t(0xc4e47d7d)),
      numext::bit_cast<float>(numext::uint32_t(0x45173159)), numext::bit_cast<float>(numext::uint32_t(0x46235731));

  TridiagonalEigenSolver<float> eigenvectors;
  eigenvectors.computeEigenvectors(diag, subdiag, expected);
  VERIFY_IS_EQUAL(eigenvectors.info(), Success);
  const Matrix3d matrix = dense_symmetric_tridiag(diag, subdiag).cast<double>();
  const Matrix3d vectors = eigenvectors.eigenvectors().cast<double>();
  VERIFY(vectors.allFinite());
  const double residual = (matrix * vectors - vectors * expected.cast<double>().asDiagonal()).cwiseAbs().maxCoeff();
  // Match the inverse-iteration suite's 128*n*eps allowance for perturbed pivots, repeated solves and normalization.
  const double tolerance = 128.0 * double(diag.size()) * double(NumTraits<float>::epsilon());
  VERIFY(residual <= tolerance * matrix.cwiseAbs().maxCoeff());
  VERIFY((vectors.transpose() * vectors - Matrix3d::Identity()).cwiseAbs().maxCoeff() <= tolerance);
}

template <typename Scalar>
void tridiagonal_eigensolver_scaling_units() {
  using Vector = Matrix<Scalar, Dynamic, 1>;
  using MatrixType = Matrix<Scalar, Dynamic, Dynamic>;
  Vector diag(4), subdiag(3);
  diag << Scalar(2), Scalar(2), Scalar(10), Scalar(10);
  subdiag << Scalar(1), Scalar(0), Scalar(2);
  Vector expected(4);
  expected << Scalar(1), Scalar(3), Scalar(8), Scalar(12);
  const MatrixType matrix = dense_symmetric_tridiag(diag, subdiag);
  const Scalar tolerance = Scalar(32) * NumTraits<Scalar>::epsilon();
  for (int exponent : {-60, 0, 60}) {
    const Scalar scale = numext::ldexp(Scalar(1), exponent);
    const Vector d = diag * scale, e = subdiag * scale;
    TridiagonalEigenSolver<Scalar> solver(d, e);
    VERIFY_IS_EQUAL(solver.info(), Success);
    const Vector values = solver.eigenvalues() / scale;
    VERIFY(values.allFinite());
    VERIFY((values - expected).cwiseAbs().maxCoeff() <= tolerance * Scalar(12));
    const MatrixType vectors = solver.eigenvectors();
    VERIFY(vectors.allFinite());
    VERIFY((matrix * vectors - vectors * values.asDiagonal()).norm() <= tolerance * matrix.norm());
    VERIFY((vectors.transpose() * vectors - MatrixType::Identity(4, 4)).norm() <= tolerance);
    solver.computeEigenvalues(
        d, e, EigenvalueRange::values(static_cast<long double>(scale) * 2, static_cast<long double>(scale) * 9));
    VERIFY_IS_EQUAL(solver.eigenvalues().size(), 2);
    VERIFY((solver.eigenvalues() / scale - expected.segment(1, 2)).cwiseAbs().maxCoeff() <= tolerance * Scalar(12));

    // The internal absolute tolerance is in input units, not normalized units.
    Vector approximate;
    const Scalar abs_tol = scale / Scalar(128);
    internal::tridiagonal_bisection(d, e, EigenvalueRange::all(), abs_tol, approximate);
    VERIFY((approximate / scale - expected).cwiseAbs().maxCoeff() <= Scalar(1) / Scalar(128) + tolerance);
  }
}

template <typename Scalar>
void tridiagonal_eigensolver_subnormal_staged() {
  if (packet_path_flushes_subnormals<Scalar>()) return;
  using Vector = Matrix<Scalar, 4, 1>;
  using Subdiag = Matrix<Scalar, 3, 1>;
  using MatrixType = Matrix<Scalar, 4, 4>;
  const Vector diag = Vector::Constant(Scalar(2));
  const Subdiag subdiag = Subdiag::Ones();
  const MatrixType matrix = dense_symmetric_tridiag(diag, subdiag);
  for (int divisor : {64, 1024, 4096}) {
    const Scalar scale = (std::numeric_limits<Scalar>::min)() / Scalar(divisor);
    const Vector d = diag * scale;
    const Subdiag e = subdiag * scale;
    TridiagonalEigenSolver<Scalar> solver;
    solver.computeEigenvalues(d, e);
    VERIFY_IS_EQUAL(solver.info(), Success);
    solver.computeEigenvectors();
    VERIFY_IS_EQUAL(solver.info(), Success);
    const Vector values = solver.eigenvalues() / scale;
    const MatrixType vectors = solver.eigenvectors();
    VERIFY(values.allFinite());
    VERIFY(vectors.allFinite());
    // Restored eigenvalues have absolute quantization denorm_min/scale. This matrix has
    // norm <= 4 and minimum eigenvalue gap 1, so both residual and orthogonality see this error.
    const Scalar quantization = std::numeric_limits<Scalar>::denorm_min() / scale;
    const Scalar tolerance = Scalar(128 * 4) * NumTraits<Scalar>::epsilon() + Scalar(8 * 4) * quantization;
    VERIFY((matrix * vectors - vectors * values.asDiagonal()).norm() <= tolerance);
    VERIFY((vectors.transpose() * vectors - MatrixType::Identity()).norm() <= tolerance);
  }
}

// diag(big, tridiag(1, 2, 1)) in float: the staged eigenvalues are accurate to eps * big, which passes the small
// block's Sturm check yet left inverse iteration with |v_i' v_j| ~ (eps * big / gap)^3, O(1) for big ~ 1e4 and
// nb = 1000. The staged pass re-bisects the block at its own scale, so its vectors are those of the block solved alone.
// (Double needs ||T|| / ||T_b|| ~ 1e7 before the shift error matters, past where the Sturm check already refines.)
template <typename Scalar>
void tridiagonal_eigensolver_staged_small_block() {
  using VectorType = Matrix<Scalar, Dynamic, 1>;
  using MatrixType = Matrix<Scalar, Dynamic, Dynamic>;
  const Index nb = 1000;
  VectorType d = VectorType::Constant(nb + 1, Scalar(2)), e = VectorType::Ones(nb);
  d(0) = Scalar(1e4);
  e(0) = Scalar(0);
  TridiagonalEigenSolver<Scalar> local, staged;
  local.computeEigenvalues(d.tail(nb).eval(), e.tail(nb - 1).eval());
  local.computeEigenvectors();
  staged.computeEigenvalues(d, e);
  staged.computeEigenvectors();
  VERIFY_IS_EQUAL(local.info(), Success);
  VERIFY_IS_EQUAL(staged.info(), Success);
  // The block's eigenvalues all lie below big, so they are the first nb columns. Bitwise equality also relies on
  // the Rayleigh-Ritz refinement leaving both results alone, which it does while their residuals are at eps.
  const MatrixType block = staged.eigenvectors().bottomLeftCorner(nb, nb);
  VERIFY_IS_EQUAL(block, local.eigenvectors());
  VERIFY_IS_EQUAL(staged.eigenvectors().topLeftCorner(1, nb).squaredNorm(), Scalar(0));
  // Vectors in different xSTEIN clusters are orthogonal to about eps * ||T_b|| / ortol, with ortol = 1e-3 ||T_b||.
  const Scalar cross_cluster_bound = Scalar(1e3) * NumTraits<Scalar>::epsilon();
  VERIFY((block.transpose() * block - MatrixType::Identity(nb, nb)).cwiseAbs().maxCoeff() <= cross_cluster_bound);
}

// Which of `excepts` f() raises, run with the floating-point environment held; -1 if it cannot be held.
template <typename Func>
int fp_exceptions_raised_by(int excepts, const Func& f) {
  std::fenv_t environment;
  if (std::feholdexcept(&environment) != 0) {
    std::cout << "SKIP: floating-point exception check: feholdexcept failed.\n";
    f();
    return -1;
  }
  f();
  const int raised = std::fetestexcept(excepts);
  VERIFY_IS_EQUAL(std::fesetenv(&environment), 0);
  return raised;
}

// The gate deciding whether to re-bisect a block's shifts compares their error bound with the block's gap without
// forming the ratio of the two, which finite input can make overflow when cubed, or 0/0. FE_DIVBYZERO is not checked:
// an exact shift's vanishing LU pivot has an infinite reciprocal by design.
template <typename Scalar>
void tridiagonal_eigensolver_block_gate_fp_exceptions() {
#if defined(FE_OVERFLOW) && defined(FE_INVALID)
  using VectorType = Matrix<Scalar, Dynamic, 1>;
  // diag(sqrt(highest), tridiag(1, 2, 1)_4), staged: (shift error / gap)^3 overflows. The block's vectors are still
  // those of the block solved alone.
  const Index nb = 4;
  VectorType d = VectorType::Constant(nb + 1, Scalar(2)), e = VectorType::Ones(nb);
  d(0) = numext::sqrt(NumTraits<Scalar>::highest());
  e(0) = Scalar(0);
  TridiagonalEigenSolver<Scalar> local, staged;
  local.computeEigenvalues(d.tail(nb).eval(), e.tail(nb - 1).eval());
  local.computeEigenvectors();
  staged.computeEigenvalues(d, e);
  int raised = fp_exceptions_raised_by(FE_OVERFLOW | FE_INVALID, [&] { staged.computeEigenvectors(); });
  if (raised >= 0) VERIFY_IS_EQUAL(raised, 0);
  VERIFY_IS_EQUAL(staged.info(), Success);
  VERIFY_IS_EQUAL(staged.eigenvectors().bottomLeftCorner(nb, nb), local.eigenvectors());

  // diag(1, [0 a; a 0]) with a = 8 denorm_min and its exact eigenvalues supplied: the block's gap underflows to zero
  // and supplied shifts carry a zero error bound, so the ratio was 0/0. (Hardware that flushes subnormal inputs splits
  // the block instead.)
  const Scalar a = Scalar(8) * std::numeric_limits<Scalar>::denorm_min();
  VectorType ds(3), es(2), w(3);
  ds << Scalar(1), Scalar(0), Scalar(0);
  es << Scalar(0), a;
  w << -a, a, Scalar(1);
  TridiagonalEigenSolver<Scalar> supplied;
  raised = fp_exceptions_raised_by(FE_OVERFLOW | FE_INVALID, [&] { supplied.computeEigenvectors(ds, es, w); });
  if (raised >= 0) VERIFY_IS_EQUAL(raised, 0);
  VERIFY_IS_EQUAL(supplied.info(), Success);
#else
  std::cout << "SKIP: block gate exception check: FE_OVERFLOW or FE_INVALID is unavailable.\n";
#endif
}

// d and e against their power-of-two scaled copies, computed beforehand, in every flush-to-zero mode: eigenvalues to
// the bisection error plus the quantization of the range they are stored in, eigenvectors (invariant under the
// scaling) against the scaled residual.
template <typename RealScalar>
void tridiagonal_check_scaled_pair(const Matrix<RealScalar, Dynamic, 1>& d, const Matrix<RealScalar, Dynamic, 1>& e,
                                   int k) {
  using MatrixType = Matrix<RealScalar, Dynamic, Dynamic>;
  using VectorType = Matrix<RealScalar, Dynamic, 1>;
  const RealScalar eps = NumTraits<RealScalar>::epsilon();
  constexpr int digits = std::numeric_limits<RealScalar>::digits;
  const Index n = d.size();
  const VectorType ds = d.unaryExpr(internal::scale_by_exponent_op<RealScalar>(k));
  const VectorType es = e.unaryExpr(internal::scale_by_exponent_op<RealScalar>(k));
  const MatrixType Ts = dense_symmetric_tridiag(ds, es);
  TridiagonalEigenSolver<RealScalar> ref;
  ref.compute(ds, es);
  VERIFY_IS_EQUAL(ref.info(), Success);

  forEachFlushToZeroMode([&](FlushToZeroMode) {
    TridiagonalEigenSolver<RealScalar> sub;
    sub.compute(d, e);
    VERIFY_IS_EQUAL(sub.info(), Success);
    const VectorType up = sub.eigenvalues().unaryExpr(internal::scale_by_exponent_op<RealScalar>(k));
    const RealScalar radius = ref.eigenvalues().cwiseAbs().maxCoeff();
    // denorm_min 2^k: the spacing of the subnormal range the eigenvalues were stored in, scaled up with them.
    const RealScalar granularity =
        numext::ldexp(RealScalar(1), std::numeric_limits<RealScalar>::min_exponent - digits + k);
    // Named before the checks: MSVC does not capture a local that a lambda first uses after a `.template` call in the
    // same expression (C2065).
    const RealScalar valueBound = RealScalar(16) * (RealScalar(n) * eps * radius + granularity);
    const RealScalar residualBound = RealScalar(64) * RealScalar(n) * eps * radius;
    VERIFY((up - ref.eigenvalues()).cwiseAbs().template maxCoeff<PropagateNaN>() <= valueBound);
    const MatrixType V = sub.eigenvectors();
    VERIFY_IS_UNITARY(V);
    VERIFY((Ts * V - V * ref.eigenvalues().asDiagonal()).cwiseAbs().template maxCoeff<PropagateNaN>() <= residualBound);
  });
}

// All-subnormal d and e, assembled from significands without floating-point arithmetic, which a SIMD unit that
// flushes subnormal inputs (ARMv7 NEON, Arm FZ, DAZ) reads as zero in the maximum that selects the normalization.
template <typename RealScalar>
void tridiagonal_eigensolver_flushed_subnormal(Index n) {
  using VectorType = Matrix<RealScalar, Dynamic, 1>;
  using Binary = internal::binary_floating_point_traits<RealScalar>;
  using Bits = typename Binary::Bits;
  constexpr int digits = std::numeric_limits<RealScalar>::digits;
  // Distinct diagonal significands of digits - 2 bits, off-diagonal ones of alternating sign.
  const Bits top = Bits(1) << (digits - 3);
  VectorType d(n), e(n - 1);
  for (Index i = 0; i < n; ++i) d(i) = numext::bit_cast<RealScalar>(top + Bits(i) * (top / Bits(16)));
  for (Index i = 0; i < n - 1; ++i) {
    const Bits sign = (i % 2 == 0) ? Bits(0) : Binary::kSignBit;
    e(i) = numext::bit_cast<RealScalar>(sign | (top / Bits(2) + Bits(i) * (top / Bits(32))));
  }
  tridiagonal_check_scaled_pair(d, e, 60);
}

// A normal maximum below min / eps with subnormal couplings: |e| flushed to zero would split the matrix into 1 x 1
// blocks with eigenvectors e_i, although the true ones are (1, +-1) / sqrt(2) for equal d.
template <typename RealScalar>
void tridiagonal_eigensolver_flushed_subnormal_coupling() {
  using VectorType = Matrix<RealScalar, Dynamic, 1>;
  constexpr int minExponent = std::numeric_limits<RealScalar>::min_exponent;
  VectorType d(4), e(3);
  d.setConstant(numext::ldexp(RealScalar(1), minExponent + 5));
  for (Index i = 0; i < 3; ++i) e(i) = numext::ldexp(RealScalar(1) + RealScalar(i) / RealScalar(4), minExponent - 4);
  tridiagonal_check_scaled_pair(d, e, 40);
}

// The staged and the supplied-eigenvalue eigenvector passes, out of line so that none of their floating-point
// operations moves across the caller's exception-flag test.
template <typename Scalar>
EIGEN_DONT_INLINE void tridiagonal_staged_and_supplied_eigenvectors(TridiagonalEigenSolver<Scalar>& staged,
                                                                    TridiagonalEigenSolver<Scalar>& supplied,
                                                                    const Matrix<Scalar, Dynamic, 1>& d,
                                                                    const Matrix<Scalar, Dynamic, 1>& e,
                                                                    const Matrix<Scalar, Dynamic, 1>& w) {
  staged.computeEigenvectors();
  supplied.computeEigenvectors(d, e, w);
}

// Assigning eigenvalues to disconnected blocks takes Sturm counts at shifts in each block's normalized units. For a
// shift far outside a small block, or within its tolerance of the overflow threshold, the count is 0 or the block's
// size, decided without forming the normalized shift, which overflows. FE_DIVBYZERO is not checked: an exact shift's
// vanishing LU pivot has an infinite reciprocal by design.
template <typename Scalar>
void tridiagonal_eigensolver_block_assignment_overflow() {
  using VectorType = Matrix<Scalar, Dynamic, 1>;
  using MatrixType = Matrix<Scalar, Dynamic, Dynamic>;
  const Scalar highest = NumTraits<Scalar>::highest(), eps = NumTraits<Scalar>::epsilon();
  const Index nb = 3;
  // T = A (+) s * tridiag(1, 2, 1)_3, with A first and last, through the staged and supplied-eigenvalue paths.
  auto check = [&](const VectorType& ad, const VectorType& ae, Scalar s) {
    const Index na = ad.size(), n = na + nb;
    for (const bool a_first : {true, false}) {
      const Index a0 = a_first ? 0 : nb, b0 = a_first ? na : 0;
      VectorType d(n), e = VectorType::Zero(n - 1);
      d.segment(a0, na) = ad;
      e.segment(a0, na - 1) = ae;
      d.segment(b0, nb).setConstant(Scalar(2) * s);
      e.segment(b0, nb - 1).setConstant(s);
      TridiagonalEigenSolver<Scalar> staged, supplied;
      staged.computeEigenvalues(d, e);
      VERIFY_IS_EQUAL(staged.info(), Success);
      const VectorType w = staged.eigenvalues();
#if defined(FE_OVERFLOW) && defined(FE_INVALID)
      const int raised = fp_exceptions_raised_by(
          FE_OVERFLOW | FE_INVALID, [&] { tridiagonal_staged_and_supplied_eigenvectors(staged, supplied, d, e, w); });
      if (raised >= 0) VERIFY_IS_EQUAL(raised, 0);
#else
      std::cout << "SKIP: block assignment overflow flag check: FE_OVERFLOW or FE_INVALID is unavailable.\n";
      tridiagonal_staged_and_supplied_eigenvectors(staged, supplied, d, e, w);
#endif
      for (const TridiagonalEigenSolver<Scalar>* es : {&staged, &supplied}) {
        VERIFY_IS_EQUAL(es->info(), Success);
        const MatrixType& V = es->eigenvectors();
        VERIFY((V.transpose() * V - MatrixType::Identity(n, n)).cwiseAbs().maxCoeff() <= Scalar(128) * Scalar(n) * eps);
        // Every vector lies in exactly one block, and A holds na of them.
        Index in_a = 0;
        for (Index j = 0; j < n; ++j) {
          const bool zero_a = (V.col(j).segment(a0, na).array() == Scalar(0)).all();
          const bool zero_b = (V.col(j).segment(b0, nb).array() == Scalar(0)).all();
          VERIFY(zero_a != zero_b);
          in_a += zero_b ? 1 : 0;
        }
        VERIFY_IS_EQUAL(in_a, na);
      }
    }
  };
  VectorType ad(1), ae(0);
  // A's shifts in units of a block at sqrt(min) overflow.
  ad << highest / Scalar(4);
  check(ad, ae, numext::sqrt((std::numeric_limits<Scalar>::min)()));
  // The largest shift plus its tolerance overflows.
  ad << highest;
  check(ad, ae, Scalar(1));
  // The same in A's own refinement check, which needs na > 1: with c = highest eps, A's largest eigenvalue is about
  // highest - 2c, within its tolerance of highest, while its row sums stay finite.
  const Scalar c = highest * eps;
  ad.resize(2);
  ae.resize(1);
  ad << highest - Scalar(2) * c, Scalar(0);
  ae << c;
  check(ad, ae, Scalar(1));
}

EIGEN_DECLARE_TEST(tridiagonal_eigensolver) {
  CALL_SUBTEST_1(tridiagonal_eigensolver_subnormal_staged<double>());
  CALL_SUBTEST_2(tridiagonal_eigensolver_subnormal_staged<float>());
  CALL_SUBTEST_2(tridiagonal_eigensolver_staged_small_block<float>());
  CALL_SUBTEST_1(tridiagonal_eigensolver_block_gate_fp_exceptions<double>());
  CALL_SUBTEST_2(tridiagonal_eigensolver_block_gate_fp_exceptions<float>());
  CALL_SUBTEST_2(tridiagonal_eigensolver_power_of_two_scaling<>());
  CALL_SUBTEST_2(tridiagonal_eigensolver_scaling_units<float>());
  CALL_SUBTEST_1(tridiagonal_eigensolver_scaling_units<double>());
  for (Index n : {Index(1), Index(2), Index(9)}) {
    EIGEN_UNUSED_VARIABLE(n);  // parts 3-6 compile no subtest that takes n
    CALL_SUBTEST_1(tridiagonal_eigensolver_flushed_subnormal<double>(n));
    CALL_SUBTEST_2(tridiagonal_eigensolver_flushed_subnormal<float>(n));
  }
  CALL_SUBTEST_1(tridiagonal_eigensolver_flushed_subnormal_coupling<double>());
  CALL_SUBTEST_2(tridiagonal_eigensolver_flushed_subnormal_coupling<float>());
  CALL_SUBTEST_1(tridiagonal_eigensolver_block_assignment_overflow<double>());
  CALL_SUBTEST_2(tridiagonal_eigensolver_block_assignment_overflow<float>());
  for (int i = 0; i < g_repeat; i++) {
    CALL_SUBTEST_1(tridiagonal_eigensolver_bisection<double>());
    CALL_SUBTEST_2(tridiagonal_eigensolver_bisection<float>());
    CALL_SUBTEST_3(tridiagonal_eigensolver_eigenvectors<double>());
    CALL_SUBTEST_4(tridiagonal_eigensolver_eigenvectors<float>());
    CALL_SUBTEST_5(tridiagonal_eigensolver_narrow<Eigen::half>());
    CALL_SUBTEST_6(tridiagonal_eigensolver_narrow<Eigen::bfloat16>());
  }
}