Skip to content

Eigen — Implementation

File Structure

eigen/
├── tiny_eigen.h        (15 lines — eigenvalue API)
├── tiny_eigen.c        (210 lines — eigenvalue implementation)
├── tiny_modal.h        (24 lines — modal analysis API)
└── tiny_modal.c        (157 lines — modal analysis implementation)

Dependencies: tiny_cfloat.h, tiny_decomp.h, <math.h>.

Eigenvalue Pipeline

2×2 Closed-Form

Simplest path — no iteration, no workspace:

float tr  = a00 + a11;
float det = (a00 - a11)*(a00 - a11) + 4*a01*a10;

if (det >= 0) {              // real eigenvalues
    lam = (tr ± sqrt(det)) / 2;
} else {                     // complex conjugate pair
    re = tr / 2;
    im = sqrt(-det) / 2;
}

Used as the base case for deflation in the Francis QR iteration.

Schur Extraction

After Francis QR iteration, \(A\) is in real Schur form (quasi-triangular):

Real Schur (n=4):
┌───┬───┬───┬───┐
│ x │ x │ x │ x │
│ 0 │ x │ x │ x │  ← sub-diagonal elements are 0 (or nearly 0)
│ 0 │ 0 │ x │ x │
│ 0 │ 0 │ x │ x │  ← 2×2 block → complex conjugate pair
└───┴───┴───┴───┘

The function walks the diagonal:

  • If T[k+1][k] ≈ 0 → real eigenvalue at T[k][k]
  • If T[k+1][k] ≠ 0 → \(2 \times 2\) block → extract complex conjugate pair via tiny_eig_2x2_f32

General Eigendecomposition

tiny_eig_general_f32(A, n, eigs, ws, ws_size)

Algorithm flow:

1. Allocate workspace (n×n matrix + n-vector + n² scratch)
2. Copy A → H (row-major, contiguous)
3. Hessenberg reduction on H  (Householder reflections)
4. Francis double-shift QR on H  (implicit shift, Wilkinson shift)
5. Extract eigenvalues from H (quasi-triangular)
6. Return eigenvalues in eigs[]

Workspace layout:

Section Size Purpose
H (work matrix) \(n \times n\) floats Copy of A, transformed in-place
Hessenberg scratch \(n\) floats Householder vector store
QR scratch \(n^2\) floats Rotation accumulators

Pole Conversion

// Continuous-time pole from discrete eigenvalue λ:
p = ln(λ) * fs

// Natural frequency and damping:
ω_n = |p|              // rad/s
f  = ω_n / (2π)        // Hz
ζ  = -Re(p) / |p|      // damping ratio (0 = undamped, 1 = critically damped)

MAC Computation

MAC(φ_a, φ_b) = |φ_a^T φ_b|² / ((φ_a^T φ_a)(φ_b^T φ_b))
  • Pair: direct dot product of two vectors, \(O(n_{\text{dof}})\)
  • Matrix: \(n_1 \times n_2\) pairwise combinations, \(O(n_1 \cdot n_2 \cdot n_{\text{dof}})\)

Stable Pole Filtering

A pole is considered "stable" when:

\[\zeta < \zeta_{\max} \quad \land \quad f_{\min} < f < f_{\max}\]

The function counts how many poles satisfy these criteria — caller then extracts the corresponding subset.

Numerical Considerations

  • Francis QR uses Wilkinson shifts for accelerated convergence — shifts are computed from the bottom-right \(2 \times 2\) sub-matrix
  • Deflation threshold: \(10^{-14} \cdot \max(1, |T_{kk}| + |T_{k+1,k+1}|)\)
  • Real matrices only: The algorithm preserves real arithmetic — no complex overhead until eigenvalue extraction
  • Float32 precision limits eigenvalue accuracy to \(\sim 10^{-4}\) for typical SSI matrices (condition \(\sim 10^6\))