Fix three bugs in SelfAdjointEigenSolver and improve test coverage

Bug fixes:

1. computeDirect 3x3: eigenvalues from computeRoots() are theoretically
   sorted via the trigonometric formula, but floating-point rounding
   (especially in float) can break the ordering. Add a 3-element sorting
   network after computeRoots() to guarantee sorted output.

2. compute(): Inf input silently returns Success with garbage results,
   unlike NaN which correctly returns NoConvergence. Add an isfinite()
   check on the scaling factor (which is maxCoeff of the matrix) to
   detect Inf/NaN early and return NoConvergence.

3. computeFromTridiagonal(): lacks the equilibration scaling that
   compute() applies, causing overflow and NoConvergence for tridiagonal
   matrices with large entries. Add the same scale-to-[-1,1] pattern.

New tests:

- Eigenvalue sorting verification (both iterative and direct solvers)
- Repeated/degenerate eigenvalues (all equal, multiplicity n-1, two
  clusters, nearly repeated separated by O(epsilon))
- Extreme eigenvalue ranges (high condition number spanning many orders
  of magnitude, near-underflow, near-overflow, mixed positive/negative,
  rank-deficient with zero eigenvalue)
- computeFromTridiagonal with large and tiny values
- Diagonal matrices (eigenvalues must match sorted diagonal)
- operatorInverseSqrt accuracy (sqrtA*invSqrtA=I, invSqrtA*A*invSqrtA=I,
  symmetry)
- RowMajor storage for computeDirect (2x2, 3x3, float, double) and
  iterative solver (dynamic RowMajor)
- Inf input detection
- Tridiagonal structure verification (off-tridiagonal entries are zero)
- Direct solver stress tests: 3x3 (near-planar covariance, triple
  eigenvalue, double eigenvalue, large off-diagonal, nearly singular)
  and 2x2 (equal eigenvalues, tiny off-diagonal, huge diagonal ratio,
  anti-diagonal dominant, negative entries)
- Tightened unitary tolerance from fixed 32*test_precision to
  4*n*epsilon (scales with matrix size)
- Fixed typo: "expponential" -> "exponential"

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

@@ -451,6 +451,13 @@ EIGEN_DEVICE_FUNC SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<Mat
   // map the matrix coefficients to [-1:1] to avoid over- and underflow.
   mat = matrix.template triangularView<Lower>();
   RealScalar scale = mat.cwiseAbs().maxCoeff();
+  if (!(numext::isfinite)(scale)) {
+    // Input contains Inf or NaN.
+    m_info = NoConvergence;
+    m_isInitialized = true;
+    m_eigenvectorsOk = false;
+    return *this;
+  }
   if (numext::is_exactly_zero(scale)) scale = RealScalar(1);
   mat.template triangularView<Lower>() /= scale;
   m_subdiag.resize(n - 1);
@@ -470,16 +477,36 @@ EIGEN_DEVICE_FUNC SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<Mat
 template <typename MatrixType>
 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::computeFromTridiagonal(
     const RealVectorType& diag, const SubDiagonalType& subdiag, int options) {
-  // TODO : Add an option to scale the values beforehand
   bool computeEigenvectors = (options & ComputeEigenvectors) == ComputeEigenvectors;
 
   m_eivalues = diag;
   m_subdiag = subdiag;
+
+  // Scale the tridiagonal matrix to [-1:1] to avoid over- and underflow,
+  // just like compute() does for the full matrix.
+  RealScalar scale = m_eivalues.cwiseAbs().maxCoeff();
+  if (m_subdiag.size() > 0) scale = numext::maxi(scale, m_subdiag.cwiseAbs().maxCoeff());
+  if (!(numext::isfinite)(scale)) {
+    m_info = NoConvergence;
+    m_isInitialized = true;
+    m_eigenvectorsOk = false;
+    return *this;
+  }
+  if (numext::is_exactly_zero(scale)) scale = RealScalar(1);
+  const bool needsScaling = scale != RealScalar(1);
+  if (needsScaling) {
+    m_eivalues /= scale;
+    m_subdiag /= scale;
+  }
+
   if (computeEigenvectors) {
     m_eivec.setIdentity(diag.size(), diag.size());
   }
   m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
 
+  // Scale back the eigenvalues.
+  if (needsScaling) m_eivalues *= scale;
+
   m_isInitialized = true;
   m_eigenvectorsOk = computeEigenvectors;
   return *this;
@@ -662,6 +689,28 @@ struct direct_selfadjoint_eigenvalues<SolverType, 3, false> {
     // compute the eigenvalues
     computeRoots(scaledMat, eivals);
 
+    // computeRoots produces theoretically sorted roots, but floating-point
+    // rounding in the trigonometric formulas can break the ordering.
+    // Enforce sorting with a branchless min/max network (3 elements).
+    {
+      Scalar tmp;
+      if (eivals(0) > eivals(1)) {
+        tmp = eivals(0);
+        eivals(0) = eivals(1);
+        eivals(1) = tmp;
+      }
+      if (eivals(1) > eivals(2)) {
+        tmp = eivals(1);
+        eivals(1) = eivals(2);
+        eivals(2) = tmp;
+      }
+      if (eivals(0) > eivals(1)) {
+        tmp = eivals(0);
+        eivals(0) = eivals(1);
+        eivals(1) = tmp;
+      }
+    }
+
     // compute the eigenvectors
     if (computeEigenvectors) {
       if ((eivals(2) - eivals(0)) <= Eigen::NumTraits<Scalar>::epsilon()) {