rtx-fea: repair the eigensolver, and stop the suite lying about the rest

Lifts the 27 `#[ignore]` markers on rtx-cfd and rtx-fea. 21 of them fail;
6 were stale, marking components that have since been implemented. The
suite now reports the truth, which means it is red.

The eigensolver had three independent defects, each individually fatal.
Found by writing closed-form tests first and confirming red:

  - The generalized reduction formed M^-1 K and ran Lanczos on it.
    M^-1 K has the right eigenvalues but is not symmetric even when K
    and M both are, and Lanczos assumes symmetry -- so it returned a
    wrong answer rather than an inaccurate one. On a 2-DOF spring-mass
    chain with M = diag(2,1) it gave 1.633 against an exact root of
    1 - sqrt(2)/2 ~= 0.293. Replaced with the Cholesky reduction
    B = L^-1 (K - sigma M) L^-T.

  - Output was unsorted. nalgebra's symmetric_eigen gives no ordering
    guarantee and none was imposed; modal analysis names modes by index,
    so the ordering is part of the contract.

  - Eigenvectors could not be transformed back out of the Krylov basis.
    The Lanczos block was (n x num_iter) and the tridiagonal
    eigenvectors (min(num_iter, k) x k); whenever those differed the
    multiply panicked on a dimension mismatch -- that is, on every
    problem with more DOFs than requested modes, which is every real
    modal analysis.

Lanczos now runs shift-invert by default. Plain Lanczos converges to the
eigenvalues of largest magnitude and modal analysis wants the lowest, so
without it the solver returns the modes nobody asked for. Also switched
to full reorthogonalization, twice per step, so converged eigenvalues do
not reappear as ghosts indistinguishable from genuine repeated roots.

ModalResults computed f = sqrt(lambda / 2pi) instead of
sqrt(lambda) / 2pi. The two agree only at lambda = 2pi, so a smoke test
asserting a positive frequency would never separate them. A
`#[cfg(disabled)]` module in the same file asserted the correct formula
-- the module was disabled rather than the bug fixed. That module is
removed; tests/eigenvalue_closed_form.rs supersedes it with every
expected value derived analytically.

Corrected a fixture rather than loosening its tolerance:
implementation_tests expected the smallest eigenvalue of
tridiag(-1, 4, -1) at order 3 to be 4 - 2 sqrt(2) ~= 1.172. The
eigenvalues of tridiag(c, a, c) are a + 2c cos(k pi / (n+1)), so the
true value is 4 - sqrt(2) ~= 2.586. The test had been quarantined for
failing to match an expectation that was never right.

rtx-fsi is untouched and stays 26/26.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>

crates/specialized/rtx-fea/src/solvers/eigenvalue.rs

@@ -29,134 +29,218 @@ impl EigenvalueSolver {
         self
     }
 
-    /// Solve generalized eigenvalue problem K*φ = λ*M*φ using Lanczos algorithm
+    /// Solve the generalized eigenvalue problem `K φ = λ M φ`.
+    ///
+    /// Returns the `num_eigenvalues` eigenvalues closest to the shift
+    /// (lowest, when no shift is set) in **ascending order**, together with
+    /// their eigenvectors as columns in the original coordinate space,
+    /// normalized to unit Euclidean length. Modal analysis names modes by
+    /// index, so the ordering is part of the contract.
+    ///
+    /// # Method
+    ///
+    /// Lanczos requires a **symmetric** operator. The obvious reduction
+    /// `M⁻¹K` is not symmetric even when `K` and `M` both are, so it cannot
+    /// be used here: it has the right eigenvalues but drives the Lanczos
+    /// recurrence to a wrong answer rather than an inaccurate one. Instead
+    /// `M` is factored as `M = L Lᵀ` and the problem reduced to the
+    /// symmetric standard form
+    ///
+    /// ```text
+    ///     B = L⁻¹ (K - σM) L⁻ᵀ,    B y = (λ - σ) y,    φ = L⁻ᵀ y
+    /// ```
+    ///
+    /// Lanczos then runs on `A = B⁻¹`, not on `B`. Plain Lanczos converges to
+    /// the eigenvalues of **largest** magnitude first, and modal analysis
+    /// wants the **lowest** — so shift-invert is not an optional refinement
+    /// here, it is what makes the answer the one that was asked for. The
+    /// eigenvalues of `A` map back as `λ = σ + 1/μ`.
+    ///
+    /// # Errors
+    ///
+    /// `M` must be symmetric positive definite — true of any consistent or
+    /// lumped mass matrix with positive densities. `K - σM` must be
+    /// non-singular; with the default shift of zero that means `K` must be
+    /// non-singular, which requires the structure to be constrained against
+    /// rigid-body motion. An unconstrained structure has zero-frequency
+    /// rigid-body modes and needs an explicit negative shift via
+    /// [`Self::with_shift`].
     pub fn solve(
         &self,
         stiffness: &SparseMatrix,
         mass: &SparseMatrix,
     ) -> FeaResult<(DVector<f64>, DMatrix<f64>)> {
         let n = stiffness.nrows();
+        if self.num_eigenvalues == 0 || self.num_eigenvalues > n {
+            return Err(SolverError::SolveError {
+                reason: format!(
+                    "requested {} eigenvalues from a system of order {n}",
+                    self.num_eigenvalues
+                ),
+            }
+            .into());
+        }
 
-        // Convert to dense for Lanczos implementation
-        let k_dense: nalgebra::DMatrix<f64> = stiffness.to_dense();
-        let m_dense: nalgebra::DMatrix<f64> = mass.to_dense();
+        let k_dense: DMatrix<f64> = stiffness.to_dense();
+        let m_dense: DMatrix<f64> = mass.to_dense();
 
-        // Apply shift-invert if specified
-        let (a_matrix, shift_value) = if let Some(sigma) = self.shift {
-            // (K - σM)⁻¹ M for shift-invert
-            let shifted: nalgebra::DMatrix<f64> = &k_dense - sigma * &m_dense;
-            let shifted_inv = shifted
-                .try_inverse()
-                .ok_or_else(|| SolverError::SolveError {
-                    reason: "Shifted matrix is singular".to_string(),
-                })?;
-            (shifted_inv * &m_dense, sigma)
-        } else {
-            // Standard eigenvalue problem: M⁻¹K
-            let m_inv = m_dense
-                .try_inverse()
-                .ok_or_else(|| SolverError::SolveError {
-                    reason: "Mass matrix is singular".to_string(),
-                })?;
-            (m_inv * &k_dense, 0.0)
-        };
+        // M = L Lᵀ. Cholesky fails exactly when M is not positive definite,
+        // which is the condition that makes the reduction below valid.
+        let chol = m_dense
+            .clone()
+            .cholesky()
+            .ok_or_else(|| SolverError::SolveError {
+                reason: "mass matrix is not symmetric positive definite".to_string(),
+            })?;
+        let l = chol.l();
 
-        // Lanczos algorithm
-        let num_iter = (self.num_eigenvalues * 2).min(n);
-        let (eigenvalues, eigenvectors) = self.lanczos(&a_matrix, num_iter)?;
+        let sigma = self.shift.unwrap_or(0.0);
+        let shifted = &k_dense - sigma * &m_dense;
 
-        // Apply shift correction if needed
-        let corrected_eigenvalues = if self.shift.is_some() {
-            eigenvalues.map(|lambda| {
-                if lambda.abs() > 1e-12 {
-                    shift_value + 1.0 / lambda
-                } else {
-                    lambda
-                }
-            })
-        } else {
-            eigenvalues
-        };
+        // B = L⁻¹ S L⁻ᵀ via two triangular solves: X = L⁻¹ S, then
+        // B = (L⁻¹ Xᵀ)ᵀ. No inverse of L is formed.
+        let x = solve_lower(&l, &shifted)?;
+        let b = solve_lower(&l, &x.transpose())?.transpose();
 
-        Ok((corrected_eigenvalues, eigenvectors))
+        // A = B⁻¹. Symmetrized to shed the asymmetry that the inverse of a
+        // numerically-symmetric matrix picks up at rounding level; Lanczos is
+        // sensitive to it.
+        let a = b.try_inverse().ok_or_else(|| SolverError::SolveError {
+            reason: if self.shift.is_some() {
+                "shifted matrix K - σM is singular".to_string()
+            } else {
+                "stiffness matrix is singular — an unconstrained structure has \
+                 rigid-body modes and needs an explicit shift"
+                    .to_string()
+            },
+        })?;
+        let a = 0.5 * (&a + &a.transpose());
+
+        // A generous Krylov subspace: the extra vectors are what separate
+        // converged Ritz values from spurious ones, and the systems reaching
+        // this solver are dense anyway.
+        let num_iter = (2 * self.num_eigenvalues + 20).min(n);
+        let (mu, y) = self.lanczos(&a, num_iter)?;
+
+        // μ = 1/(λ - σ), so λ = σ + 1/μ. Ritz values that have not separated
+        // from zero correspond to eigenvalues infinitely far from the shift;
+        // they are not converged modes and are dropped.
+        let mut pairs: Vec<(f64, usize)> = mu
+            .iter()
+            .enumerate()
+            .filter(|(_, m)| m.abs() > 1e-12)
+            .map(|(i, m)| (sigma + 1.0 / m, i))
+            .collect();
+
+        if pairs.len() < self.num_eigenvalues {
+            return Err(SolverError::SolveError {
+                reason: format!(
+                    "only {} of {} requested modes converged",
+                    pairs.len(),
+                    self.num_eigenvalues
+                ),
+            }
+            .into());
+        }
+
+        // Ascending in λ: with shift-invert these are the modes nearest σ.
+        pairs.sort_by(|a, b| a.0.total_cmp(&b.0));
+        pairs.truncate(self.num_eigenvalues);
+
+        let eigenvalues = DVector::from_iterator(pairs.len(), pairs.iter().map(|&(l, _)| l));
+
+        // φ = L⁻ᵀ y takes the eigenvectors back out of the Cholesky factor's
+        // coordinates into the original space.
+        let selected = DMatrix::from_fn(y.nrows(), pairs.len(), |i, j| y[(i, pairs[j].1)]);
+        let mut eigenvectors = solve_upper(&l.transpose(), &selected)?;
+        for mut col in eigenvectors.column_iter_mut() {
+            let norm = col.norm();
+            if norm > 1e-300 {
+                col /= norm;
+            }
+        }
+
+        Ok((eigenvalues, eigenvectors))
     }
 
-    /// Lanczos algorithm implementation
+    /// Lanczos iteration on a symmetric operator.
+    ///
+    /// Returns **all** Ritz pairs the Krylov subspace produced — every
+    /// eigenvalue of the projected tridiagonal matrix, with eigenvectors
+    /// already lifted back into the operator's space. Selecting among them is
+    /// the caller's job, because which ones are wanted depends on whether a
+    /// shift-invert transformation is in play.
+    ///
+    /// Full reorthogonalization against every previous Lanczos vector is used
+    /// rather than the three-term recurrence alone. Without it the vectors
+    /// lose orthogonality in floating point and the iteration returns
+    /// duplicated ("ghost") copies of converged eigenvalues, which is
+    /// indistinguishable from a genuine repeated root — exactly the failure a
+    /// modal analysis must not have.
     fn lanczos(
         &self,
         a: &DMatrix<f64>,
         num_iter: usize,
     ) -> FeaResult<(DVector<f64>, DMatrix<f64>)> {
         let n = a.nrows();
+        debug_assert!(num_iter <= n);
 
-        // Initialize with random vector
+        // A fixed, non-symmetric start vector. It must not be orthogonal to
+        // any eigenvector of interest; a constant vector can be, by symmetry
+        // of the mode shapes, so the entries are varied. Fixed rather than
+        // random keeps the result reproducible run to run.
         let mut v = DVector::from_fn(n, |i, _| ((i as f64 + 1.0) * 0.1).sin());
-        v.normalize_mut();
+        let norm = v.norm();
+        if norm < 1e-300 {
+            return Err(SolverError::SolveError {
+                reason: "Lanczos start vector vanished".to_string(),
+            }
+            .into());
+        }
+        v /= norm;
 
-        // Lanczos vectors and tridiagonal matrix coefficients
-        let mut lanczos_vectors = DMatrix::zeros(n, num_iter);
-        let mut alpha = DVector::zeros(num_iter);
-        let mut beta = DVector::zeros(num_iter);
+        let mut basis = DMatrix::zeros(n, num_iter);
+        let mut alpha = vec![0.0; num_iter];
+        let mut beta = vec![0.0; num_iter];
 
-        lanczos_vectors.column_mut(0).copy_from(&v);
+        // Number of iterations actually completed — the subspace can be
+        // exhausted early when the start vector lies in an invariant
+        // subspace, and the projected matrix must then be truncated to match.
+        let mut m = num_iter;
+
+        basis.column_mut(0).copy_from(&v);
         let mut w = a * &v;
         alpha[0] = v.dot(&w);
         w -= alpha[0] * &v;
 
         for j in 1..num_iter {
             beta[j - 1] = w.norm();
-
             if beta[j - 1] < 1e-12 {
-                // Early termination if Krylov subspace is exhausted
+                m = j;
                 break;
             }
 
-            let v_prev = v.clone();
             v = &w / beta[j - 1];
-            lanczos_vectors.column_mut(j).copy_from(&v);
+            basis.column_mut(j).copy_from(&v);
 
             w = a * &v;
             alpha[j] = v.dot(&w);
-            w -= alpha[j] * &v + beta[j - 1] * &v_prev;
+            w -= alpha[j] * &v;
 
-            // Re-orthogonalization for numerical stability
-            for k in 0..=j {
-                let vk = lanczos_vectors.column(k);
-                let h = w.dot(&vk);
-                w -= h * vk;
+            // Full reorthogonalization, twice. One pass is not always enough
+            // when the loss of orthogonality is severe.
+            for _ in 0..2 {
+                for k in 0..=j {
+                    let vk = basis.column(k);
+                    let h = w.dot(&vk);
+                    w -= h * vk;
+                }
             }
         }
 
-        // Solve tridiagonal eigenvalue problem
-        let (tri_eigenvalues, tri_eigenvectors) =
-            self.solve_tridiagonal(&alpha, &beta, self.num_eigenvalues)?;
-
-        // Transform eigenvectors back to original space
-        let eigenvectors = &lanczos_vectors * &tri_eigenvectors;
-
-        // Normalize eigenvectors
-        let mut normalized_eigenvectors = DMatrix::zeros(n, self.num_eigenvalues);
-        for i in 0..self.num_eigenvalues {
-            let mut col = eigenvectors.column(i).clone_owned();
-            col.normalize_mut();
-            normalized_eigenvectors.column_mut(i).copy_from(&col);
-        }
-
-        Ok((tri_eigenvalues, normalized_eigenvectors))
-    }
-
-    /// Solve tridiagonal eigenvalue problem using QR algorithm
-    fn solve_tridiagonal(
-        &self,
-        alpha: &DVector<f64>,
-        beta: &DVector<f64>,
-        num_eigenvalues: usize,
-    ) -> FeaResult<(DVector<f64>, DMatrix<f64>)> {
-        let n = alpha.len().min(num_eigenvalues);
-
-        // Build tridiagonal matrix
-        let mut t = DMatrix::zeros(n, n);
-        for i in 0..n {
+        // Project onto the Krylov subspace actually built.
+        let mut t = DMatrix::zeros(m, m);
+        for i in 0..m {
             t[(i, i)] = alpha[i];
             if i > 0 {
                 t[(i, i - 1)] = beta[i - 1];
@@ -164,18 +248,36 @@ impl EigenvalueSolver {
             }
         }
 
-        // Use nalgebra's symmetric eigendecomposition
         let eigen = t.symmetric_eigen();
 
-        // Extract requested number of smallest eigenvalues and eigenvectors
-        let num_vals = num_eigenvalues.min(n);
-        let eigenvalues = DVector::from_fn(num_vals, |i, _| eigen.eigenvalues[i]);
-        let eigenvectors = DMatrix::from_fn(n, num_vals, |i, j| eigen.eigenvectors[(i, j)]);
+        // Lift the Ritz vectors out of the Krylov basis: x = V y.
+        let ritz_values = eigen.eigenvalues.clone();
+        let ritz_vectors = basis.columns(0, m) * &eigen.eigenvectors;
 
-        Ok((eigenvalues, eigenvectors))
+        Ok((ritz_values, ritz_vectors))
     }
 }
 
+/// Solve `L X = B` for `X`, with `L` lower triangular.
+fn solve_lower(l: &DMatrix<f64>, b: &DMatrix<f64>) -> FeaResult<DMatrix<f64>> {
+    l.clone().solve_lower_triangular(b).ok_or_else(|| {
+        SolverError::SolveError {
+            reason: "triangular solve failed: zero on the diagonal".to_string(),
+        }
+        .into()
+    })
+}
+
+/// Solve `U X = B` for `X`, with `U` upper triangular.
+fn solve_upper(u: &DMatrix<f64>, b: &DMatrix<f64>) -> FeaResult<DMatrix<f64>> {
+    u.clone().solve_upper_triangular(b).ok_or_else(|| {
+        SolverError::SolveError {
+            reason: "triangular solve failed: zero on the diagonal".to_string(),
+        }
+        .into()
+    })
+}
+
 /// Modal analysis results.
 #[derive(Debug, Clone)]
 pub struct ModalResults {
@@ -186,8 +288,25 @@ pub struct ModalResults {
 }
 
 impl ModalResults {
+    /// Derive frequencies and periods from the eigenvalues of `K φ = λ M φ`.
+    ///
+    /// `λ` is `ω²` in rad²/s², so `ω = √λ` and `f = √λ / (2π)`. Note this is
+    /// *not* `√(λ / 2π)` — the two agree only at `λ = 2π`, so a smoke test
+    /// that checks the frequency is merely positive will not tell them apart.
+    ///
+    /// A negative eigenvalue is not a frequency. It means the reduced
+    /// stiffness is not positive definite — buckling, or an unconverged
+    /// nonlinear step — and the corresponding entry is reported as NaN rather
+    /// than silently squared away into a plausible-looking number.
     pub fn from_eigenvalues(eigenvalues: DVector<f64>, eigenvectors: DMatrix<f64>) -> Self {
-        let frequencies = eigenvalues.map(|lambda| (lambda / (2.0 * std::f64::consts::PI)).sqrt());
+        let two_pi = 2.0 * std::f64::consts::PI;
+        let frequencies = eigenvalues.map(|lambda| {
+            if lambda < 0.0 {
+                f64::NAN
+            } else {
+                lambda.sqrt() / two_pi
+            }
+        });
         let periods = frequencies.map(|f| if f > 1e-12 { 1.0 / f } else { f64::INFINITY });
 
         Self {
@@ -199,78 +318,7 @@ impl ModalResults {
     }
 }
 
-#[cfg(disabled)]
-mod tests {
-    use super::*;
-    use approx::assert_relative_eq;
-
-    fn create_test_matrices() -> (SparseMatrix, SparseMatrix) {
-        // Create simple test matrices
-        // K = [4, -1, 0; -1, 4, -1; 0, -1, 4]  (tridiagonal stiffness)
-        // M = I (identity mass matrix)
-        let mut k = SparseMatrix::new(3, 3);
-        k.set(0, 0, 4.0).unwrap();
-        k.set(0, 1, -1.0).unwrap();
-        k.set(1, 0, -1.0).unwrap();
-        k.set(1, 1, 4.0).unwrap();
-        k.set(1, 2, -1.0).unwrap();
-        k.set(2, 1, -1.0).unwrap();
-        k.set(2, 2, 4.0).unwrap();
-
-        let mut m = SparseMatrix::new(3, 3);
-        m.set(0, 0, 1.0).unwrap();
-        m.set(1, 1, 1.0).unwrap();
-        m.set(2, 2, 1.0).unwrap();
-
-        (k, m)
-    }
-
-    #[test]
-    fn test_eigenvalue_solver() {
-        let (k, m) = create_test_matrices();
-        let solver = EigenvalueSolver::new(3);
-
-        let (eigenvalues, eigenvectors) = solver.solve(&k, &m).unwrap();
-
-        // Expected eigenvalues for this system (analytically known)
-        // λ₁ ≈ 2.0 - sqrt(2) ≈ 0.586
-        // λ₂ = 4.0
-        // λ₃ ≈ 2.0 + sqrt(2) ≈ 3.414
-
-        assert_eq!(eigenvalues.len(), 3);
-        assert_eq!(eigenvectors.ncols(), 3);
-
-        // Check that eigenvalues are sorted
-        for i in 1..eigenvalues.len() {
-            assert!(eigenvalues[i] >= eigenvalues[i - 1]);
-        }
-
-        // Check smallest eigenvalue is close to expected
-        assert_relative_eq!(eigenvalues[0], 4.0 - 2.0 * 2.0_f64.sqrt(), epsilon = 1e-6);
-    }
-
-    #[test]
-    fn test_modal_results() {
-        let eigenvalues = DVector::from_vec(vec![1.0, 4.0, 9.0]);
-        let eigenvectors = DMatrix::identity(3, 3);
-
-        let modal = ModalResults::from_eigenvalues(eigenvalues, eigenvectors);
-
-        // Check frequency calculation: f = sqrt(λ) / (2π)
-        assert_relative_eq!(
-            modal.frequencies[0],
-            1.0 / (2.0 * std::f64::consts::PI),
-            epsilon = 1e-10
-        );
-        assert_relative_eq!(
-            modal.frequencies[1],
-            2.0 / (2.0 * std::f64::consts::PI),
-            epsilon = 1e-10
-        );
-        assert_relative_eq!(
-            modal.frequencies[2],
-            3.0 / (2.0 * std::f64::consts::PI),
-            epsilon = 1e-10
-        );
-    }
-}
+// Tests live in tests/eigenvalue_closed_form.rs, where every expected value
+// is derived analytically. The `#[cfg(disabled)]` module that used to sit
+// here asserted a smallest eigenvalue of `4 - 2*sqrt(2)` for `tridiag(-1, 4,
+// -1)`, which is not an eigenvalue of that matrix at all.