1// Shift-invert Lanczos for the eigenvalues of the pencil (K, M) in one
2// slice [lo, hi] of the spectrum. With M diagonal the pencil is turned into
3// the standard symmetric problem A = M^{-1/2} K M^{-1/2}, and Lanczos is
4// run on (A - sigma I)^{-1} = M^{1/2} (K - sigma M)^{-1} M^{1/2} with the
5// banded factorization of banded.ts, full reorthogonalization, and a check
6// of the count of converged eigenvalues in the slice against the inertia
7// counts at its ends. The mode shapes returned are M^{-1/2} times the
8// eigenvectors of A, that is, mass-normalized.
10import { factorize, solve, type Factorization } from "./banded.ts";
12export interface SliceResult {
13 /** Eigenvalues omega^2 found in (lo, hi], ascending. */
14 lambda: number[];
15 /** Mass-normalized mode shapes, one Float64Array of length n per eigenvalue. */
16 shapes: Float64Array[];
17 /** Number of eigenvalues the inertia counts say lie in the slice. */
18 expected: number;
19 /** Lanczos steps used. */
20 steps: number;
21}
23export interface SliceOptions {
24 /** Relative tolerance on the shift-inverted Ritz values; 1e-9 by default. */
25 tol?: number;
26 /** Largest number of Lanczos steps per start; 3 count + 40 by default. */
27 maxSteps?: number;
28 /** Restarts with a fresh random vector when modes are missing (degenerate pairs); 3 by default. */
29 restarts?: number;
30 seed?: number;
31}
33/**
34 * Eigenvalues of (K, M) in (lo, hi] and their shapes. `countBelowLo` and
35 * `countBelowHi` are the inertia counts at the two ends, so that the
36 * expected number of eigenvalues is known.
37 */
38export function sliceModes(
39 band: Float64Array,
40 mass: Float64Array,
41 n: number,
42 bw: number,
43 lo: number,
44 hi: number,
45 countBelowLo: number,
46 countBelowHi: number,
47 opts: SliceOptions = {},
48): SliceResult {
49 const expected = countBelowHi - countBelowLo;
50 if (expected <= 0) return { lambda: [], shapes: [], expected, steps: 0 };
51 const tol = opts.tol ?? 1e-9;
52 const restarts = opts.restarts ?? 3;
53 const sigma = 0.5 * (lo + hi);
54 const fact = factorize(band, mass, n, bw, sigma);
55 const sqm = new Float64Array(n);
56 for (let i = 0; i < n; i++) sqm[i] = Math.sqrt(mass[i]);
57 let seed = (opts.seed ?? 4321) >>> 0;
58 const rnd = () => {
59 seed = (seed * 1664525 + 1013904223) >>> 0;
60 return seed / 4294967296 - 0.5;
61 };
63 const locked: { lambda: number; x: Float64Array }[] = [];
64 let steps = 0;
65 const maxSteps = Math.min(n, opts.maxSteps ?? 3 * expected + 40);
66 for (let attempt = 0; attempt <= restarts && locked.length < expected; attempt++) {
67 const found = lanczosRun(fact, sqm, n, sigma, lo, hi, expected - locked.length, locked.map((l) => l.x), maxSteps, tol, rnd);
68 steps += found.steps;
69 for (const f of found.pairs) locked.push(f);
70 }
71 locked.sort((a, b) => a.lambda - b.lambda);
72 const shapes = locked.map((l) => {
73 const phi = new Float64Array(n);
74 for (let i = 0; i < n; i++) phi[i] = l.x[i] / sqm[i];
75 return phi;
76 });
77 return { lambda: locked.map((l) => l.lambda), shapes, expected, steps };
78}
80function dot(a: Float64Array, b: Float64Array): number {
81 let s = 0;
82 for (let i = 0; i < a.length; i++) s += a[i] * b[i];
83 return s;
84}
86/** One Lanczos run from a random start, orthogonal to the locked vectors; returns the converged pairs in (lo, hi]. */
87function lanczosRun(
88 fact: Factorization,
89 sqm: Float64Array,
90 n: number,
91 sigma: number,
92 lo: number,
93 hi: number,
94 wanted: number,
95 locked: Float64Array[],
96 maxSteps: number,
97 tol: number,
98 rnd: () => number,
99): { pairs: { lambda: number; x: Float64Array }[]; steps: number } {
100 const Q: Float64Array[] = [];
101 const alpha: number[] = [], beta: number[] = [];
102 const tmp = new Float64Array(n);
103 // Classical Gram-Schmidt against the basis, repeated once when the first
104 // pass removed a large part of the vector ("twice is enough"); the second
105 // pass is rarely needed and the basis is the main memory traffic.
106 const orthogonalize = (r: Float64Array, against: Float64Array[]): void => {
107 if (against.length === 0) return;
108 const before = dot(r, r);
109 for (let pass = 0; pass < 2; pass++) {
110 for (const q of against) {
111 const c = dot(r, q);
112 if (c !== 0) for (let i = 0; i < n; i++) r[i] -= c * q[i];
113 }
114 if (dot(r, r) > 0.5 * before) break;
115 }
116 };
117 let q = new Float64Array(n);
118 for (let i = 0; i < n; i++) q[i] = rnd();
119 orthogonalize(q, locked);
120 let nrm = Math.sqrt(dot(q, q));
121 for (let i = 0; i < n; i++) q[i] /= nrm;
122 Q.push(q);
124 let pairs: { lambda: number; x: Float64Array }[] = [];
125 let steps = 0;
126 const checkEvery = Math.max(6, Math.round(wanted / 10));
127 const check = (): boolean => {
128 const m = Q.length;
129 const { values, vectors } = tridiagonalEigen(alpha, beta.slice(0, m - 1));
130 const betaM = beta[m - 1] ?? 0;
131 const conv: { lambda: number; s: Float64Array }[] = [];
132 for (let k = 0; k < m; k++) {
133 const theta = values[k];
134 if (theta === 0) continue;
135 const lambda = sigma + 1 / theta;
136 if (!(lambda > lo && lambda <= hi)) continue;
137 const err = Math.abs(betaM * vectors[k][m - 1]);
138 if (err <= tol * Math.abs(theta)) conv.push({ lambda, s: vectors[k] });
139 }
140 if (conv.length < wanted && m < maxSteps) return false;
141 pairs = conv.map(({ lambda, s }) => {
142 const x = new Float64Array(n);
143 for (let j = 0; j < m; j++) {
144 const sj = s[j], qj = Q[j];
145 if (sj === 0) continue;
146 for (let i = 0; i < n; i++) x[i] += sj * qj[i];
147 }
148 return { lambda, x };
149 });
150 return true;
151 };
153 for (let j = 0; j < maxSteps; j++) {
154 steps++;
155 // r = Op q_j with Op = M^{1/2} (K - sigma M)^{-1} M^{1/2}
156 for (let i = 0; i < n; i++) tmp[i] = sqm[i] * Q[j][i];
157 solve(fact, tmp, tmp);
158 const r = new Float64Array(n);
159 for (let i = 0; i < n; i++) r[i] = sqm[i] * tmp[i];
160 if (j > 0) for (let i = 0; i < n; i++) r[i] -= beta[j - 1] * Q[j - 1][i];
161 const a = dot(r, Q[j]);
162 alpha.push(a);
163 for (let i = 0; i < n; i++) r[i] -= a * Q[j][i];
164 orthogonalize(r, Q);
165 orthogonalize(r, locked);
166 const b = Math.sqrt(dot(r, r));
167 beta.push(b);
168 const last = j === maxSteps - 1 || Q.length >= n;
169 if (b < 1e-14 * Math.max(1, Math.abs(a))) {
170 // Invariant subspace: everything reachable has converged.
171 check();
172 break;
173 }
174 // Fewer than `wanted` steps cannot hold `wanted` converged Ritz values,
175 // and the tridiagonal eigenproblem costs m^3, so check only from there
176 // on and at a spacing that grows with the slice.
177 const m = j + 1;
178 if (last || (m >= wanted && (m - wanted) % checkEvery === 0)) {
179 if (check()) break;
180 }
181 if (Q.length >= n) break;
182 for (let i = 0; i < n; i++) r[i] /= b;
183 Q.push(r);
184 }
185 if (pairs.length === 0 && steps > 0) check();
186 return { pairs, steps };
187}
189/**
190 * Eigenvalues and eigenvectors of the symmetric tridiagonal matrix with
191 * diagonal `d` and off-diagonal `e` (implicit QL, after EISPACK tql2).
192 * vectors[k] is the eigenvector of values[k].
193 */
194export function tridiagonalEigen(dIn: number[], eIn: number[]): { values: number[]; vectors: Float64Array[] } {
195 const n = dIn.length;
196 const d = dIn.slice();
197 const e = eIn.slice();
198 e.push(0);
199 const z: Float64Array[] = [];
200 for (let i = 0; i < n; i++) {
201 const row = new Float64Array(n);
202 row[i] = 1;
203 z.push(row);
204 }
205 let f = 0, tst1 = 0;
206 const eps = 2 ** -52;
207 for (let l = 0; l < n; l++) {
208 tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l]));
209 let m = l;
210 while (m < n) {
211 if (Math.abs(e[m]) <= eps * tst1) break;
212 m++;
213 }
214 if (m > l) {
215 let iter = 0;
216 do {
217 iter++;
218 if (iter > 200) throw new Error("tridiagonal QL did not converge");
219 let g = d[l];
220 let p = (d[l + 1] - g) / (2 * e[l]);
221 let r = Math.hypot(p, 1);
222 if (p < 0) r = -r;
223 d[l] = e[l] / (p + r);
224 d[l + 1] = e[l] * (p + r);
225 const dl1 = d[l + 1];
226 let h = g - d[l];
227 for (let i = l + 2; i < n; i++) d[i] -= h;
228 f += h;
229 p = d[m];
230 let c = 1, c2 = c, c3 = c;
231 const el1 = e[l + 1];
232 let s = 0, s2 = 0;
233 for (let i = m - 1; i >= l; i--) {
234 c3 = c2; c2 = c; s2 = s;
235 g = c * e[i];
236 h = c * p;
237 r = Math.hypot(p, e[i]);
238 e[i + 1] = s * r;
239 s = e[i] / r;
240 c = p / r;
241 p = c * d[i] - s * g;
242 d[i + 1] = h + s * (c * g + s * d[i]);
243 for (let k = 0; k < n; k++) {
244 h = z[k][i + 1];
245 z[k][i + 1] = s * z[k][i] + c * h;
246 z[k][i] = c * z[k][i] - s * h;
247 }
248 }
249 p = (-s * s2 * c3 * el1 * e[l]) / dl1;
250 e[l] = s * p;
251 d[l] = c * p;
252 } while (Math.abs(e[l]) > eps * tst1);
253 }
254 d[l] += f;
255 e[l] = 0;
256 }
257 const order = Array.from({ length: n }, (_, i) => i).sort((a, b) => d[a] - d[b]);
258 const values = order.map((i) => d[i]);
259 const vectors = order.map((i) => {
260 const v = new Float64Array(n);
261 for (let k = 0; k < n; k++) v[k] = z[k][i];
262 return v;
263 });
264 return { values, vectors };
265}