TridiagonalInverseIteration.h

// 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
 
#ifndef EIGEN_TRIDIAGONAL_INVERSE_ITERATION_H
#define EIGEN_TRIDIAGONAL_INVERSE_ITERATION_H
 
#include "./SelfAdjointEigenSolver.h"
// For tridiagonal_sturm_count_below(), used to assign eigenvalues to disconnected blocks.
#include "./TridiagonalBisection.h"
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
namespace internal {

struct inverse_iteration_rng {
  numext::uint64_t state;
 
  explicit inverse_iteration_rng(numext::uint64_t seed) : state(seed) {}
 
  template <typename RealScalar>
  RealScalar next() {
    // SplitMix64: a strong avalanche so that consecutive per-column seeds yield well-separated
    // (uncorrelated) start vectors -- consecutive LCG seeds would produce near-parallel starts that
    // collapse under the cluster reorthogonalization.
    state += 0x9E3779B97F4A7C15ULL;
    const numext::uint64_t z = splitmix64_mix(state);
    // Top 32 bits give a uniform integer in [0, 2^32); map to [0, 1) and then to (-1, 1).
    // RealScalar must be float or wider: a narrower scalar overflows RealScalar(hi) to infinity
    // (inf * 0 = NaN start vectors), which is why TridiagonalEigenSolver computes scalars narrower
    // than float in float (see its ComputeScalar) and none of this code ever runs on them.
    const numext::uint32_t hi = numext::uint32_t(z >> 32);
    const RealScalar u = RealScalar(hi) * (RealScalar(1) / RealScalar(4294967296.0));
    return RealScalar(2) * u - RealScalar(1);
  }
};

template <typename RealScalar>
void tridiagonal_lagtf(RealScalar* d, RealScalar* du, RealScalar* dl, RealScalar* du2, Index* piv, RealScalar lambda,
                       Index n) {
  d[0] -= lambda;
  piv[n - 1] = 0;
  if (n == 1) return;
 
  RealScalar scale1 = numext::abs(d[0]) + numext::abs(du[0]);
  for (Index k = 0; k < n - 1; ++k) {
    d[k + 1] -= lambda;
    RealScalar scale2 = numext::abs(dl[k]) + numext::abs(d[k + 1]);
    if (k < n - 2) scale2 += numext::abs(du[k + 1]);
    // Scaled magnitudes of the two candidate pivots; pick the larger (partial pivoting).
    const RealScalar piv1 = numext::is_exactly_zero(d[k]) ? RealScalar(0) : numext::abs(d[k]) / scale1;
    if (numext::is_exactly_zero(dl[k])) {
      piv[k] = 0;
      scale1 = scale2;
      if (k < n - 2) du2[k] = RealScalar(0);
    } else {
      const RealScalar piv2 = numext::abs(dl[k]) / scale2;
      if (piv2 <= piv1) {
        // No interchange: eliminate using the diagonal pivot d[k] (non-zero, since piv2 > 0).
        piv[k] = 0;
        scale1 = scale2;
        dl[k] = dl[k] / d[k];
        d[k + 1] -= dl[k] * du[k];
        if (k < n - 2) du2[k] = RealScalar(0);
      } else {
        // Interchange rows k and k+1; the sub-diagonal entry becomes the pivot.
        piv[k] = 1;
        const RealScalar mult = d[k] / dl[k];
        d[k] = dl[k];
        const RealScalar temp = d[k + 1];
        d[k + 1] = du[k] - mult * temp;
        if (k < n - 2) {
          du2[k] = du[k + 1];
          du[k + 1] = -mult * du2[k];
        }
        du[k] = temp;
        dl[k] = mult;
      }
    }
  }
}

template <typename RealScalar>
void tridiagonal_lagts(const RealScalar* d, const RealScalar* rcp, const RealScalar* du, const RealScalar* dl,
                       const RealScalar* du2, const Index* piv, RealScalar* b, Index n) {
  const RealScalar eps = NumTraits<RealScalar>::epsilon();
  const RealScalar sfmin = (std::numeric_limits<RealScalar>::min)();
  const RealScalar bignum = RealScalar(1) / sfmin;
 
  // Perturbation floor: eps times the largest magnitude entry of U.
  using ArrayBuf = Map<const Array<RealScalar, Dynamic, 1>>;
  RealScalar tol = ArrayBuf(d, n).abs().maxCoeff();
  if (n > 1) tol = numext::maxi(tol, ArrayBuf(du, n - 1).abs().maxCoeff());
  if (n > 2) tol = numext::maxi(tol, ArrayBuf(du2, n - 2).abs().maxCoeff());
  tol *= eps;
  if (numext::is_exactly_zero(tol)) tol = eps;
 
  // Forward substitution: apply P then L^{-1}. Written branchlessly in the (data-dependent, roughly
  // 50/50) pivot-interchange flag so a misprediction cannot stall the serial recurrence; both arms
  // are evaluated and selected, reproducing the two branches bit for bit.
  for (Index k = 1; k < n; ++k) {
    const RealScalar dlk = dl[k - 1];
    const RealScalar bkm1 = b[k - 1];
    const RealScalar bk = b[k];
    const bool interchange = piv[k - 1] != 0;
    b[k - 1] = interchange ? bk : bkm1;
    b[k] = interchange ? (bkm1 - dlk * bk) : (bk - dlk * bkm1);
  }
 
  // Back substitution: solve U x = b. Away from the singularity the pivot is safe and the quotient is
  // a reciprocal multiply (rcp[k] = 1/d[k], exact and precomputed once for all the solves); only a
  // dangerously small pivot -- the rare near-singular element of the deliberately singular system --
  // falls back to the LAPACK xLAGTS perturbation/rescale with a true division. The guard is written so
  // the fast path is taken exactly when the original would divide by the unperturbed pivot, and is
  // virtually always predicted taken.
  for (Index k = n - 1; k >= 0; --k) {
    RealScalar temp = b[k];
    if (k <= n - 2) temp -= du[k] * b[k + 1];
    if (k <= n - 3) temp -= du2[k] * b[k + 2];
    const RealScalar ak = d[k];
    const RealScalar absak = numext::abs(ak);
    if (EIGEN_PREDICT_TRUE(absak >= RealScalar(1) || (absak >= sfmin && numext::abs(temp) <= absak * bignum))) {
      b[k] = temp * rcp[k];
    } else {
      // Tiny pivot: perturb (or rescale) exactly as LAPACK xLAGTS so the quotient cannot overflow.
      RealScalar a = ak;
      RealScalar t = temp;
      RealScalar pert = (a >= RealScalar(0)) ? tol : -tol;
      while (true) {
        const RealScalar aa = numext::abs(a);
        if (aa < RealScalar(1)) {
          if (aa < sfmin) {
            if (numext::is_exactly_zero(aa) || numext::abs(t) * sfmin > aa) {
              a += pert;
              pert *= RealScalar(2);
              continue;
            } else {
              t *= bignum;
              a *= bignum;
            }
          } else if (numext::abs(t) > aa * bignum) {
            a += pert;
            pert *= RealScalar(2);
            continue;
          }
        }
        break;
      }
      b[k] = t / a;
    }
  }
}

template <typename RealScalar, typename EivecType>
Index tridiagonal_inverse_iteration_block(const RealScalar* sdiag, const RealScalar* ssub, const RealScalar* xj_scaled,
                                          const Index* clstart, Index n, RealScalar onenrm, RealScalar dtpcrt,
                                          int maxits, int extra, EivecType& eivecs, Index j_lo, Index j_hi) {
  using RealVectorType = Matrix<RealScalar, Dynamic, 1>;
  const RealScalar eps = NumTraits<RealScalar>::epsilon();
 
  // Work arrays for the LU factors of T - xj*I (reused across the block's columns) and the iterate.
  // lu_rcp holds the exact reciprocal of the U diagonal, computed once per factorization and reused by
  // the inverse-iteration back-solves. Being locals, each thread gets its own copy.
  RealVectorType lu_d(n), lu_dl(n), lu_du(n), lu_du2(n), lu_rcp(n), b(n);
  Matrix<Index, Dynamic, 1> piv(n);
 
  Index nonconv = 0;  // columns that exhausted maxits without satisfying the growth test
  for (Index j = j_lo; j < j_hi; ++j) {
    const RealScalar xj = xj_scaled[j];
    const Index gpind = clstart[j];  // first column of j's cluster (>= j_lo by block alignment)
 
    // Factor T - xj*I = P L U once; the inverse-iteration loop reuses it.
    lu_d = Map<const RealVectorType>(sdiag, n);
    lu_du.head(n - 1) = Map<const RealVectorType>(ssub, n - 1);
    lu_dl.head(n - 1) = Map<const RealVectorType>(ssub, n - 1);
    tridiagonal_lagtf<RealScalar>(lu_d.data(), lu_du.data(), lu_dl.data(), lu_du2.data(), piv.data(), xj, n);
    // Exact reciprocals of the U diagonal (pdiv, not the approximate preciprocal of cwiseInverse());
    // a vanishing pivot yields inf, which the back-solve fast path never consumes.
    lu_rcp.array() = RealScalar(1) / lu_d.array();
 
    // Deterministic pseudo-random start vector (seeded per column for thread independence).
    inverse_iteration_rng rng(numext::uint64_t(j) + 1);
    for (Index i = 0; i < n; ++i) b[i] = rng.template next<RealScalar>();
 
    int nrmchk = 0;
    bool converged = false;
    for (int its = 0; its < maxits; ++its) {
      // Scale the right-hand side so the near-singular solve neither overflows nor underflows.
      const RealScalar bmax = b.cwiseAbs().maxCoeff();
      if (numext::is_exactly_zero(bmax)) break;  // degenerate iterate; cannot grow -> not converged
      const RealScalar target = RealScalar(n) * onenrm * numext::maxi(eps, numext::abs(lu_d[n - 1]));
      const RealScalar scl = target / bmax;
      if (EIGEN_PREDICT_TRUE(scl >= (std::numeric_limits<RealScalar>::min)())) {
        b *= scl;
      } else {
        // scl underflowed (huge iterate from a shift within ~pivmin of an exact eigenvalue). A
        // subnormal or flushed-to-zero scale would zero the iterate on flush-to-zero hardware
        // (ARMv7 NEON), so scale in two in-range steps instead.
        b /= bmax;
        b *= target;
      }
      tridiagonal_lagts<RealScalar>(lu_d.data(), lu_rcp.data(), lu_du.data(), lu_dl.data(), lu_du2.data(), piv.data(),
                                    b.data(), n);
 
      // Modified Gram-Schmidt against the already-accepted eigenvectors of this cluster.
      for (Index i = gpind; i < j; ++i) b -= b.dot(eivecs.col(i)) * eivecs.col(i);
 
      const RealScalar nrm = b.cwiseAbs().maxCoeff();
      if (nrm < dtpcrt) continue;          // not yet grown into the eigenvector; iterate
      if (++nrmchk < extra + 1) continue;  // a couple of safety iterations past the threshold
      converged = true;
      break;
    }
    // LAPACK xSTEIN reports vectors that never satisfied the growth test; mirror that. The column is
    // still written (best effort), but the caller surfaces the count as ComputationInfo::NoConvergence.
    if (!converged) ++nonconv;
 
    // Normalize to unit 2-norm with a deterministic sign (largest-magnitude entry positive). The
    // iterate can carry huge entries (~1/eps times the start) on a nearly singular solve, so scale
    // before taking the norm -- otherwise squaredNorm() would overflow and zero out the vector.
    Index jmax = 0;
    const RealScalar binf = b.cwiseAbs().maxCoeff(&jmax);
    safe_scaling<RealScalar>::scale_to(b, b, binf);
    const RealScalar nrm2 = b.norm();
    RealScalar scl = numext::is_exactly_zero(nrm2) ? RealScalar(1) : RealScalar(1) / nrm2;
    if (b[jmax] < RealScalar(0)) scl = -scl;
    eivecs.col(j) = b * scl;
  }
  return nonconv;
}

template <typename DiagType, typename SubdiagType, typename EivalType, typename EivecType>
Index tridiagonal_inverse_iteration_connected(const DiagType& diag, const SubdiagType& subdiag, const EivalType& eivals,
                                              EivecType& eivecs) {
  using RealScalar = typename DiagType::Scalar;
  EIGEN_STATIC_ASSERT(NumTraits<RealScalar>::IsInteger == 0 && NumTraits<RealScalar>::IsComplex == 0,
                      THIS_FUNCTION_IS_NOT_FOR_INTEGER_OR_COMPLEX_TYPES)
 
  const Index n = diag.size();
  const Index m = eivals.size();
  if (n == 0 || m == 0) return 0;
  if (n == 1) {
    eivecs.setOnes();
    return 0;
  }
 
  const RealScalar eps = NumTraits<RealScalar>::epsilon();
 
  // Normalize T (and the shifts) to O(1) so the deliberately near-singular factor/solve cannot
  // overflow or underflow; eigenvectors are invariant under this uniform scaling.
  // The rescan recovers an all-subnormal input that a flushing SIMD unit reads as zero.
  const RealScalar maxCoeff = max_preserving_subnormals(
      safe_scaling<RealScalar>::recover_flushed_max_coeff(diag, diag.cwiseAbs().maxCoeff()),
      safe_scaling<RealScalar>::recover_flushed_max_coeff(subdiag, subdiag.cwiseAbs().maxCoeff()));
  Matrix<RealScalar, Dynamic, 1> sdiag(n), ssub(n - 1);
  const auto factors = safe_scaling<RealScalar>::scale_to(sdiag, diag, maxCoeff);
  safe_scaling<RealScalar>::scale_to(ssub, subdiag, maxCoeff, factors);
  Matrix<RealScalar, Dynamic, 1> xj_scaled(m);
  safe_scaling<RealScalar>::scale_to(xj_scaled, eivals, maxCoeff, factors);
  if (maxCoeff > RealScalar(0) && maxCoeff < (std::numeric_limits<RealScalar>::min)()) {
    // The first scale is clamped to normal range. Finish normalization in a second finite step
    // so the absolute pivot floor and growth threshold still see an O(1) matrix.
    const RealScalar scaledMax = numext::maxi(sdiag.cwiseAbs().maxCoeff(), ssub.cwiseAbs().maxCoeff());
    const auto remaining = safe_scaling<RealScalar>::scale_in_place(sdiag, scaledMax);
    safe_scaling<RealScalar>::scale_in_place(ssub, scaledMax, remaining);
    safe_scaling<RealScalar>::scale_in_place(xj_scaled, scaledMax, remaining);
  }
 
  // Infinity norm of the scaled T: max_i (|e_{i-1}| + |d_i| + |e_i|), missing boundary off-diagonals zero.
  RealScalar onenrm = numext::abs(sdiag[0]) + numext::abs(ssub[0]);
  onenrm = numext::maxi(onenrm, numext::abs(sdiag[n - 1]) + numext::abs(ssub[n - 2]));
  if (n > 2)
    onenrm = numext::maxi(onenrm, (ssub.head(n - 2).array().abs() + sdiag.segment(1, n - 2).array().abs() +
                                   ssub.segment(1, n - 2).array().abs())
                                      .maxCoeff());
  if (numext::is_exactly_zero(onenrm)) onenrm = RealScalar(1);  // T == 0: any orthonormal basis works
 
  // Cluster threshold and convergence threshold (LAPACK xSTEIN constants).
  const RealScalar ortol = RealScalar(1e-3) * onenrm;
  const RealScalar dtpcrt = numext::sqrt(RealScalar(0.1) / RealScalar(n));
  const int maxits = 5;
  const int extra = 2;
 
  // Pre-pass (sequential, O(m)): the perturbed, scaled shift xj_scaled[j] and the cluster start
  // clstart[j] (index of the first eigenvalue of j's cluster). Reproduced exactly as the serial sweep
  // would, so the parallel column blocks below match the serial result bit for bit. A shift sitting on
  // top of its predecessor is nudged up by pertol so the two factorizations stay distinct; a gap larger
  // than ortol starts a new cluster (within which the eigenvectors are reorthogonalized). The
  // perturbation is at most ~pertol << ortol, so it never moves a column across a cluster boundary.
  Matrix<Index, Dynamic, 1> clstart(m);
  {
    Index gpind = 0;
    RealScalar xjm = RealScalar(0);  // previous (possibly perturbed) shift
    for (Index j = 0; j < m; ++j) {
      RealScalar xj = xj_scaled[j];
      if (j > 0) {
        // The xSTEIN nudge 10*eps*|xj| separates coincident shifts so their factorizations differ.
        // Capped at a fraction of the cluster threshold: at low precision (bfloat16: 10*eps ~ 0.08)
        // the un-capped nudge can push a shift across a cluster boundary into a neighbouring
        // eigenspace. Shifts the cap leaves coincident still yield an orthonormal basis, via the
        // per-cluster Gram-Schmidt on independent random starts.
        const RealScalar pertol = numext::mini(RealScalar(10) * numext::abs(eps * xj), RealScalar(0.25) * ortol);
        if (xj - xjm < pertol) xj = xjm + pertol;
      }
      if (j == 0 || numext::abs(xj - xjm) > ortol) gpind = j;
      xj_scaled[j] = xj;
      clstart[j] = gpind;
      xjm = xj;
    }
  }
 
  const RealScalar* sdiag_p = sdiag.data();
  const RealScalar* ssub_p = ssub.data();
  const RealScalar* xj_p = xj_scaled.data();
  const Index* cl_p = clstart.data();
 
  // Distinct clusters are independent, so the columns split across threads with a single fork/join:
  // thread t owns a contiguous slice of columns snapped to cluster boundaries (so no cluster straddles
  // two threads) and writes only those columns of eivecs. No communication until the join, and the
  // result is bit-identical to the serial path for any thread count. Each thread tallies its own
  // non-converged columns; the reduction sums them into the value returned to the caller.
  Index nonconv = 0;
#if defined(EIGEN_HAS_OPENMP)
  int nthreads = 1;
  // Don't nest inside an existing parallel region, and only fork when there is enough work to amortize
  // the thread overhead. Each column costs ~one O(n) factorization plus a few O(n) back-solves.
  if (omp_get_num_threads() == 1) {
    // Work ~ m columns * n rows (a factorization plus a few O(n) back-solves each). kMinTaskSize is the
    // minimum such work to give each thread before adding one: m*n / kMinTaskSize threads, capped at
    // the pool size. The value is measured -- at m*n = 2*kMinTaskSize (the two-thread point) inverse
    // iteration is already ~2x faster than serial, with the gain growing to the core count for larger n.
    const double work = double(m) * double(n);
    const double kMinTaskSize = 2048.0;
    const Index work_threads = Index(work / kMinTaskSize);
    nthreads = int(numext::maxi(Index(1), numext::mini(work_threads, Index(Eigen::nbThreads()))));
  }
  if (nthreads > 1) {
#pragma omp parallel num_threads(nthreads) reduction(+ : nonconv)
    {
      const Index nt = omp_get_num_threads();
      const Index tid = omp_get_thread_num();
      // Balanced split snapped up to the next cluster boundary (clstart[j]==j marks a cluster start).
      // Adjacent threads snap the same raw index to the same boundary, so the blocks tile [0, m)
      // exactly; a thread whose whole share falls inside one cluster simply gets an empty block.
      Index lo = tid * m / nt;
      Index hi = (tid + 1) * m / nt;
      while (lo < m && cl_p[lo] != lo) ++lo;
      while (hi < m && cl_p[hi] != hi) ++hi;
      nonconv += tridiagonal_inverse_iteration_block<RealScalar>(sdiag_p, ssub_p, xj_p, cl_p, n, onenrm, dtpcrt, maxits,
                                                                 extra, eivecs, lo, hi);
    }
  } else
#endif
  {
    nonconv = tridiagonal_inverse_iteration_block<RealScalar>(sdiag_p, ssub_p, xj_p, cl_p, n, onenrm, dtpcrt, maxits,
                                                              extra, eivecs, 0, m);
  }
  return nonconv;
}

template <typename DiagType, typename SubdiagType, typename EivalType, typename EivecType>
Index tridiagonal_inverse_iteration(const DiagType& diag, const SubdiagType& subdiag, const EivalType& eivals,
                                    EivecType& eivecs) {
  using RealScalar = typename DiagType::Scalar;
  using VectorType = Matrix<RealScalar, Dynamic, 1>;
  const Index n = diag.size();
  const Index m = eivals.size();
  if (n == 0 || m == 0) return 0;
 
  // A matrix whose largest entry is below the recovery threshold min / eps can hold significant subnormal
  // couplings, which FTZ/DAZ hardware reads as zero in the comparisons below (ARMv7 NEON also in the packet maxima).
  // Scale it into the normal range first, exactly, through integer significands: the eigenvectors are invariant
  // under the scaling, and every representable entry becomes normal. The test reads the exponent from the
  // representation, since the maximum itself may be subnormal and a floating-point comparison on it is what FTZ/DAZ
  // breaks; a zero or non-finite maximum has no exponent in range and is left alone.
  {
    RealScalar maxCoeff = safe_scaling<RealScalar>::recover_flushed_max_coeff(diag, diag.cwiseAbs().maxCoeff());
    if (n >= 2) {
      maxCoeff = max_preserving_subnormals(
          maxCoeff, safe_scaling<RealScalar>::recover_flushed_max_coeff(subdiag, subdiag.cwiseAbs().maxCoeff()));
    }
    // 2^(e - 1) <= maxCoeff < 2^e, so maxCoeff < 2^recovery iff e <= recovery; e == 0 for a zero maximum and one
    // above the largest finite exponent for infinities and NaN.
    const int e = frexp_exponent_preserving_subnormals(maxCoeff);
    if (e != 0 && e <= safe_scaling<RealScalar>::subnormal_recovery_exponent()) {
      VectorType sdiag(n), ssub(n - 1), seivals(m);
      const auto factors = safe_scaling<RealScalar>::scale_to(sdiag, diag, maxCoeff);
      safe_scaling<RealScalar>::scale_to(ssub, subdiag, maxCoeff, factors);
      safe_scaling<RealScalar>::scale_to(seivals, eivals, maxCoeff, factors);
      return tridiagonal_inverse_iteration(sdiag, ssub, seivals, eivecs);
    }
  }
 
  // Split at negligible couplings (cf. xSTEBZ): |e_k| <= eps * sqrt(|d_k|) * sqrt(|d_k+1|). The
  // geometric-mean form is scale-invariant on its own (both sides scale linearly) and needs no
  // additive floor: a floor expressed at any single scale falsely splits strongly connected blocks
  // living far below that scale. Taking square roots before multiplying keeps every intermediate in
  // range, so the comparison stays exact down to subnormal entries; a zero coupling always splits.
  const RealScalar eps = NumTraits<RealScalar>::epsilon();
  Matrix<Index, Dynamic, 1> bstart(n + 1);  // block b spans rows [bstart(b), bstart(b+1))
  Index nblocks = 1;
  bstart(0) = 0;
  if (n > 1) {
    const Array<bool, Dynamic, 1> split = (subdiag.array().abs() <= eps * (diag.head(n - 1).array().abs().sqrt() *
                                                                           diag.tail(n - 1).array().abs().sqrt()));
    for (Index k = 0; k + 1 < n; ++k)
      if (split(k)) bstart(nblocks++) = k + 1;
  }
  bstart(nblocks) = n;
 
  if (nblocks == 1) return tridiagonal_inverse_iteration_connected(diag, subdiag, eivals, eivecs);
 
  // Per-block normalized data (block-local scale) for the assignment Sturm counts, and the
  // per-block localization tolerance btol: the count's rounding-displacement bound at the block's
  // own scale. The subset-edge windows below are widened by btol so that a requested value equal to
  // a block eigenvalue (up to that block's rounding) is found, without a window at any other
  // block's scale swallowing well-separated foreign eigenvalues -- a tolerance expressed at the
  // global scale would hand a small block's eigenvalue to whichever block comes first.
  const RealScalar safemin = numext::maxi(RealScalar(1) / NumTraits<RealScalar>::highest(),
                                          (RealScalar(1) + eps) * (std::numeric_limits<RealScalar>::min)());
  VectorType alpha_all(n), beta_sq_all(n), bscale(nblocks), bpivmin(nblocks), btol(nblocks);
  RealScalar gscale = RealScalar(0);
  for (Index b = 0; b < nblocks; ++b) {
    const Index b0 = bstart(b), nb = bstart(b + 1) - b0;
    RealScalar s = diag.segment(b0, nb).cwiseAbs().maxCoeff();
    if (nb > 1) s = numext::maxi(s, subdiag.segment(b0, nb - 1).cwiseAbs().maxCoeff());
    if (numext::is_exactly_zero(s)) s = RealScalar(1);
    gscale = numext::maxi(gscale, s);
    auto alpha = alpha_all.segment(b0, nb);
    const auto factors = safe_scaling<RealScalar>::scale_to(alpha, diag.segment(b0, nb), s);
    bscale(b) = factors.scale;
    RealScalar max_bsq = RealScalar(0);
    if (nb > 1) {
      auto beta = beta_sq_all.segment(b0, nb - 1);
      safe_scaling<RealScalar>::scale_to(beta, subdiag.segment(b0, nb - 1), s, factors);
      beta = beta.array().square();
      max_bsq = beta.maxCoeff();
    }
    bpivmin(b) = safemin * numext::maxi(max_bsq, RealScalar(1));
    // In original units, block row sums are bounded by 3*s, independently of the chosen scaling factor.
    btol(b) = RealScalar(2.1) * (RealScalar(3) * RealScalar(nb) * eps + RealScalar(4) * safemin) * s;
  }
 
  // Assign each requested eigenvalue to a block by capacity. Groups are runs of exactly-equal
  // requested values; a block's capacity for a group is its Sturm count over the group's value
  // interval, bounded by the midpoints to the neighbouring distinct values and, at the subset's
  // edges, widened by a tolerance. The edge tolerance is tiered: the per-block btol first, so that
  // exactly supplied eigenvalues can only be claimed by the block that owns them, and -- when the
  // edge group still cannot cover its copies -- the matrix-scale gtol, which admits requested
  // values whose error is at the scale of the whole matrix (e.g. the staged bisection's output for
  // a block much smaller than the matrix norm). The intervals partition the requested span and
  // adjacent intervals share their boundary evaluation, so each block eigenvalue is counted exactly
  // once even at count knife-edges. A copy its own interval cannot supply is carried into the next
  // interval; copies left at the end pair up with the blocks holding leftover span capacity (a
  // miscount that let a block absorb a foreign copy freed exactly one such slot elsewhere).
  const RealScalar gtol = RealScalar(2.1) * (RealScalar(3) * RealScalar(n) * eps + RealScalar(4) * safemin) * gscale;
  Matrix<Index, Dynamic, 1> blockof(m), assigned(nblocks), caps(nblocks), below_prev(nblocks), below_cur(nblocks),
      below_edge(nblocks), carry(m), carry_next(m);
  assigned.setZero();
  Matrix<Index, Dynamic, 1> local_index(m);
  Array<bool, Dynamic, 1> claimed = Array<bool, Dynamic, 1>::Constant(n, false);
  // Sturm count over block b of its eigenvalues strictly below the unnormalized shift x (normalized
  // by the block's own scale, as the block's alpha/beta data is).
  auto count_below = [&](Index b, RealScalar x) -> Index {
    const Index b0 = bstart(b), nb = bstart(b + 1) - b0;
    return tridiagonal_sturm_count_below<RealScalar>(alpha_all.data() + b0, beta_sq_all.data() + b0, nb, bpivmin(b),
                                                     x / bscale(b));
  };
  for (Index b = 0; b < nblocks; ++b) below_prev(b) = count_below(b, eivals[0] - btol(b));
  below_edge = below_prev;
  Index ncarry = 0;
  Index g_begin = 0;
  while (g_begin < m) {
    Index g_end = g_begin + 1;
    while (g_end < m && !(eivals[g_begin] < eivals[g_end])) ++g_end;
    // Upper boundary: midpoint to the next distinct value, or the per-block widened edge at the end.
    const bool first = (g_begin == 0);
    const bool last = (g_end >= m);
    const RealScalar mid_bound =
        last ? RealScalar(0) : RealScalar(0.5) * eivals[g_end - 1] + RealScalar(0.5) * eivals[g_end];
    bool widened = false;
    while (true) {
      Index total = 0;
      for (Index b = 0; b < nblocks; ++b) {
        const RealScalar bound = last ? eivals[m - 1] + (widened ? gtol : btol(b)) : mid_bound;
        below_cur(b) = count_below(b, bound);
        caps(b) = numext::maxi(Index(0), below_cur(b) - below_prev(b));
        total += caps(b);
      }
      // Retry an edge group once with the matrix-scale tolerance if the block-scale windows left it
      // short; interior deficits are handled by the carry chain instead.
      if (widened || !(first || last) || total >= g_end - g_begin + (last ? ncarry : 0)) break;
      widened = true;
      if (first) {
        for (Index b = 0; b < nblocks; ++b) below_prev(b) = count_below(b, eivals[0] - gtol);
        below_edge = below_prev;
      }
    }
    // The group's own copies first, then the carried ones; whatever finds no capacity carries on.
    Index ncarry_next = 0;
    for (Index pass = 0; pass < 2; ++pass) {
      const Index count = (pass == 0) ? g_end - g_begin : ncarry;
      for (Index c = 0; c < count; ++c) {
        const Index j = (pass == 0) ? g_begin + c : carry(c);
        Index chosen = -1;
        for (Index b = 0; b < nblocks && chosen < 0; ++b)
          if (caps(b) > 0) chosen = b;
        if (chosen < 0) {
          carry_next(ncarry_next++) = j;
          continue;
        }
        const Index index = below_cur(chosen) - caps(chosen);
        local_index(j) = index;
        claimed(bstart(chosen) + index) = true;
        --caps(chosen);
        ++assigned(chosen);
        blockof(j) = chosen;
      }
    }
    carry.swap(carry_next);
    ncarry = ncarry_next;
    below_prev = below_cur;
    g_begin = g_end;
  }
  // Leftover copies: first blocks with leftover span capacity, then any block with spare rows.
  for (Index c = 0; c < ncarry; ++c) {
    const Index j = carry(c);
    Index chosen = -1;
    for (Index b = 0; b < nblocks && chosen < 0; ++b)
      if (below_cur(b) - below_edge(b) - assigned(b) > 0) chosen = b;
    if (chosen < 0)
      for (Index b = 0; b < nblocks && chosen < 0; ++b)
        if (bstart(b + 1) - bstart(b) - assigned(b) > 0) chosen = b;
    if (chosen < 0) chosen = 0;
    {
      const Index b0 = bstart(chosen), nb = bstart(chosen + 1) - b0;
      Index index = numext::mini(below_edge(chosen), nb - 1);
      while (index < nb && claimed(b0 + index)) ++index;
      if (index == nb) {
        index = 0;
        while (index < nb && claimed(b0 + index)) ++index;
      }
      if (index == nb) {
        eivecs.setZero();
        return m;
      }
      local_index(j) = index;
      claimed(b0 + index) = true;
    }
    ++assigned(chosen);
    blockof(j) = chosen;
  }
 
  // Run each block independently and scatter its columns; rows outside the block stay exactly zero.
  eivecs.setZero();
  Index nonconv = 0;
  VectorType wloc, refined, endpoints;
  Matrix<numext::int64_t, Dynamic, 1> counts;
  Matrix<Index, Dynamic, 1> indices;
  Matrix<RealScalar, Dynamic, Dynamic> vloc;
  Matrix<Index, Dynamic, 1> colmap(m);
  for (Index b = 0; b < nblocks; ++b) {
    const Index b0 = bstart(b), nb = bstart(b + 1) - b0;
    Index mb = 0;
    for (Index j = 0; j < m; ++j)
      if (blockof(j) == b) colmap(mb++) = j;
    if (mb == 0) continue;
    wloc.resize(mb);
    for (Index k = 0; k < mb; ++k) wloc(k) = eivals[colmap(k)];
    vloc.resize(nb, mb);
    const VectorType bdiag = diag.segment(b0, nb);
    const VectorType bsub = subdiag.segment(b0, nb > 1 ? nb - 1 : 0);
    if (nb > 1) {
      // Retain the capacity assignment, including multiplicities across blocks. As in LAPACK
      // xSTEBZ, use block-local bisection accuracy rather than a global reorthogonalization floor.
      indices.resize(mb);
      for (Index k = 0; k < mb; ++k) indices(k) = local_index(colmap(k));
      std::sort(indices.data(), indices.data() + mb);
      endpoints.resize(2 * mb);
      counts.resize(2 * mb);
      endpoints.head(mb) = (wloc.array() - btol(b)) / bscale(b);
      endpoints.tail(mb) = (wloc.array() + btol(b)) / bscale(b);
      tridiagonal_sturm_counts(alpha_all.data() + b0, beta_sq_all.data() + b0, nb, bpivmin(b), endpoints.data(),
                               counts.data(), 2 * mb);
      bool needs_refinement = false;
      // A supplied shift is locally resolved if its rounding-sized interval contains the
      // assigned eigenvalue index. This also detects repeated coarse shifts claiming one root.
      for (Index k = 0; k < mb && !needs_refinement; ++k) {
        needs_refinement = counts(k) > indices(k) || counts(mb + k) <= indices(k);
      }
      if (needs_refinement) {
        for (Index first = 0; first < mb;) {
          Index last = first + 1;
          while (last < mb && indices(last) == indices(last - 1) + 1) ++last;
          tridiagonal_bisection(bdiag, bsub, EigenvalueRange::indices(indices(first), indices(last - 1) + 1),
                                RealScalar(0), refined);
          wloc.segment(first, last - first) = refined;
          first = last;
        }
      }
    }
    nonconv += tridiagonal_inverse_iteration_connected(bdiag, bsub, wloc, vloc);
    for (Index k = 0; k < mb; ++k) eivecs.col(colmap(k)).segment(b0, nb) = vloc.col(k);
  }
  return nonconv;
}

template <typename DiagType, typename SubdiagType, typename EivalType, typename EivecType>
void tridiagonal_rayleigh_ritz_refine(const DiagType& diag, const SubdiagType& subdiag, const EivalType& eivals,
                                      EivecType& eivecs) {
  using RealScalar = typename DiagType::Scalar;
  using DenseType = Matrix<RealScalar, Dynamic, Dynamic>;
  const Index n = diag.size();
  const Index m = eivals.size();
  if (n < 2 || m < 2) return;
 
  // Inf-norm of T and the cluster threshold, matching tridiagonal_inverse_iteration()'s grouping.
  RealScalar onenrm = numext::abs(diag[0]) + numext::abs(subdiag[0]);
  onenrm = numext::maxi(onenrm, numext::abs(diag[n - 1]) + numext::abs(subdiag[n - 2]));
  for (Index i = 1; i < n - 1; ++i)
    onenrm = numext::maxi(onenrm, numext::abs(subdiag[i - 1]) + numext::abs(diag[i]) + numext::abs(subdiag[i]));
  const RealScalar ortol = RealScalar(1e-3) * onenrm;
  if (!(ortol > RealScalar(0))) return;  // T == 0: nothing to refine
  const RealScalar inv_onenrm = RealScalar(1) / onenrm;
  const RealScalar refine_threshold = RealScalar(16) * NumTraits<RealScalar>::epsilon();
 
  SelfAdjointEigenSolver<DenseType> block_solver;
  // Per-cluster scratch, hoisted out of the loop and reused across clusters: each assignment below
  // resizes only when a cluster is larger than any seen so far, so the common many-small-clusters case
  // allocates once rather than per cluster. Bsym is the distinct target for the symmetrization, which
  // also lets B + B.transpose() avoid an aliasing temporary.
  DenseType Vc, TVc, B, Bsym;
  Index s = 0;
  while (s < m) {
    Index t = s;
    while (t + 1 < m && (eivals[t + 1] - eivals[t]) < ortol) ++t;
    const Index c = t - s + 1;
    if (c > 1) {
      // Project T onto the cluster basis. T V_c via the tridiagonal structure:
      //   (T V_c)(i,:) = e_{i-1} V_c(i-1,:) + d_i V_c(i,:) + e_i V_c(i+1,:).
      Vc = eivecs.middleCols(s, c);  // copy so the write-back below cannot alias
      TVc = diag.asDiagonal() * Vc;
      TVc.topRows(n - 1) += subdiag.asDiagonal() * Vc.bottomRows(n - 1);
      TVc.bottomRows(n - 1) += subdiag.asDiagonal() * Vc.topRows(n - 1);
      // Worst per-vector residual in the cluster (scaled by 1/||T|| so it cannot overflow). Skip the
      // refinement when the cluster is already at machine precision -- refining it would only add the
      // subspace error to vectors that are individually more accurate.
      const RealScalar cluster_resid =
          ((TVc - Vc * eivals.segment(s, c).asDiagonal()) * inv_onenrm).colwise().norm().maxCoeff();
      if (cluster_resid > refine_threshold) {
        B.noalias() = Vc.transpose() * TVc;
        Bsym = B + B.transpose();  // enforce exact symmetry (distinct target avoids an aliasing temp)
        Bsym *= RealScalar(0.5);
        block_solver.compute(Bsym, ComputeEigenvectors);
        eivecs.middleCols(s, c).noalias() = Vc * block_solver.eigenvectors();
      }
    }
    s = t + 1;
  }
}
 
}  // namespace internal
}  // namespace Eigen
 
#endif  // EIGEN_TRIDIAGONAL_INVERSE_ITERATION_H