StructuredMatrices: Add DPR1EigenSolver (secular-equation eigensolver)

Add DPR1EigenSolver, a direct \mathcal{O}(n^2) eigensolver for real symmetric diagonal-plus-rank-one matrices A = D + \rho z z^T: the standalone version of the secular-equation kernel (xLAED2/3/4) at the heart of the divide-and-conquer symmetric eigensolvers.

New files only, independent of the other open StructuredMatrices merge requests; phase 9 and the final piece of the module build-out.

Scope

The scalar type is limited at compile time to float, double or long double. Low-precision and custom scalar types are rejected rather than quietly given a solver whose error analysis does not hold for them.

What makes it accurate

Eigenvectors come from Gu-Eisenstat weights (SIMAX 15, 1994): the z-vector for which the computed roots are exact secular eigenvalues, with ratio-paired products against overflow. That is what makes the eigenvector matrix numerically orthogonal to \mathcal{O}(n \varepsilon) without any reorthogonalization. Deflation is the two-stage xLAED2 scheme, and the Givens rotations it generates are recorded and replayed per eigenvector in \mathcal{O}(n) rather than accumulated into a dense matrix, so the whole solve stays \mathcal{O}(n^2).

Exponent range

The other half of the work is exponent range. The problem is scaled by an exact power of two into a safe range before the secular solve and rescaled once at output, which is exact, so results are bit-identical for well-scaled inputs, where the scale is one. A ~30k-problem before/after battery — random, wide dynamic range, clustered, and sparse z, at several \rho scales — confirmed that: every eigenvalue on previously-converging inputs is bit-identical.

Algorithm

Each stage is documented in the header.

• \rho < 0 is reduced by negating the matrix (A \to -A = (-D) + (-\rho) z z^T), and the diagonal is sorted. z is prescaled by an exact power of two before stableNorm, and \rho \|z\|^2 is assembled from frexp mantissas plus a single exact ldexp. A z whose true norm exceeds the overflow threshold is therefore still accepted whenever \rho \|z\|^2 is itself representable — for example z = [\text{DBL\_MAX}, \text{DBL\_MAX}] with subnormal \rho.
• The whole problem is scaled by an exact power of two into a safe range. The eigenvalues of (2^{-e} D) + (2^{-e} \rho) z z^T are exactly 2^{-e} times the originals, so the secular solve runs on the scaled data and the spectrum is rescaled once at output. Full-exponent poles such as d = \pm 0.75 \cdot \text{DBL\_MAX} with \rho = 10^{300} no longer overflow internally.
• Two-stage deflation with the xLAED2 threshold 8 \varepsilon \max(\|d\|_\infty, \rho \|z\|^2), evaluated on the scaled data. Negligible z entries deflate outright; nearly equal poles are combined by Givens rotations whose dropped coupling sits below the backward-error budget.
• Secular roots by bisection in pole-shifted coordinates. The shift is the nearer pole, chosen by the sign of f at the interval midpoint, and the last root is bracketed entirely in shifted coordinates on (0, \rho \|\zeta\|^2]. Every distance \lambda - d_i is an exact data difference plus a small offset, so there is no cancellation. Bisection stops when the bracket endpoints become adjacent floats, which is magnitude-independent; the iteration cap is a pure backstop sized to the worst theoretical bracket (\text{expRange} + 2 \cdot \text{digits} + 32) and reports NoConvergence if it is ever reached.
• Eigenvectors from Gu-Eisenstat weights, as described above.
• Non-finite inputs are rejected up front. At extreme exponents, representability of the spectrum is certified from the original inputs by evaluating the secular equation at the largest finite value in double-word arithmetic:

\rho \sum_i \frac{z_i^2}{M - d_i} \le 1

A proven overflow reports InvalidInput; an unresolved boundary case reports NoConvergence. This avoids both the false-overflow and the false-success classifications that normalized-root rounding introduces.

Benchmark

contrib/benchmarks/StructuredMatrices/bench_structured_dpr1.cpp compares DPR1EigenSolver against dense SelfAdjointEigenSolver on the materialized matrix, both full and eigenvalues-only.

Appendix: test coverage (structured_dpr1, 5 subtests)

Eigenvalues are checked against SelfAdjointEigenSolver with the bound taken relative to the data scale \|D\| + |\rho| \|z\|^2. That is the secular formulation's backward-error scale and exceeds \|A\| under cancellation, so it is the honest bound. Residual and \mathcal{O}(n \varepsilon) orthogonality are checked alongside.

• Tight clusters and exact duplicates, exercising the Givens path.
• The Gu-Eisenstat litmus case: close-but-undeflated poles with gaps sweeping \varepsilon to \sqrt{\varepsilon}, where naive eigenvectors measurably lose orthogonality.
• Sparse and zero z; both signs of \rho and \rho = 0; rank-one A; unsorted input; wide dynamic range; n = 128 for growth.
• Float huge-z handled by the exact prescaling. The regression constant is finite now — the previous 4e38f was already infinity, which made the check vacuous — and a genuinely overflowing \rho\|z\|^2 is reported as InvalidInput.
• NaN and Inf rejection; the all-deflated huge-diagonal case; agreement between the eigenvalues-only and full paths.
• The review regressions with the reviewers' exact inputs: full-range poles, z = [\text{DBL\_MAX}, \text{DBL\_MAX}] with \rho = 10^{-310}, the near-deflation-boundary bisection case d = [0, 1], z = [16 \varepsilon, 1], the 1x1 overflow contract, maximum-range double-rounding cases for both signs of \rho, and rejection of invalid eigensolver options.

🤖 Generated with Claude Code