1// Symmetric banded LDL^T factorization of K - sigma M, its solve, and its
2// inertia (the number of negative pivots, which by Sylvester's law equals
3// the number of eigenvalues of the pencil (K, M) below sigma). The band
4// layout is that of assemble.ts: row i holds K[i, i - d] at offset d.
6export interface Factorization {
7 n: number;
8 bw: number;
9 /** Unit lower band factor L, same layout as the input (the diagonal offset is unused). */
10 L: Float64Array;
11 /** Pivots D. */
12 D: Float64Array;
13 /** Number of negative pivots: eigenvalues of (K, M) below the shift. */
14 negative: number;
15 shift: number;
16}
18/** LDL^T of K - shift * diag(mass). Throws on an exactly zero pivot (a shift equal to an eigenvalue). */
19export function factorize(band: Float64Array, mass: Float64Array, n: number, bw: number, shift: number): Factorization {
20 const s = bw + 1;
21 const L = new Float64Array(band);
22 for (let i = 0; i < n; i++) L[i * s] -= shift * mass[i];
23 const D = new Float64Array(n);
24 // W[j][k] = D[k] L[j][k], kept alongside L so that the inner loop of the
25 // factorization multiplies two loads rather than three.
26 const W = new Float64Array(n * s);
27 let negative = 0;
28 for (let i = 0; i < n; i++) {
29 const rowI = i * s;
30 const kmin = Math.max(0, i - bw);
31 for (let j = kmin; j < i; j++) {
32 const rowJ = j * s;
33 let acc = L[rowI + (i - j)];
34 const k0 = Math.max(kmin, j - bw);
35 // acc -= sum_k L[i][k] W[j][k]; L[i][k] sits at rowI + i - k, W[j][k] at rowJ + j - k.
36 const offI = rowI + i, offJ = rowJ + j;
37 for (let k = k0; k < j; k++) acc -= L[offI - k] * W[offJ - k];
38 L[rowI + (i - j)] = acc / D[j];
39 }
40 let d = L[rowI];
41 for (let k = kmin; k < i; k++) {
42 const l = L[rowI + (i - k)];
43 const w = l * D[k];
44 W[rowI + (i - k)] = w;
45 d -= l * w;
46 }
47 if (d === 0 || !isFinite(d)) throw new Error(`zero pivot at row ${i} for shift ${shift}`);
48 D[i] = d;
49 if (d < 0) negative++;
50 }
51 return { n, bw, L, D, negative, shift };
52}
54/** Solves (K - shift M) x = b in place (x may be b). */
55export function solve(f: Factorization, b: Float64Array, x: Float64Array): void {
56 const { n, bw, L, D } = f;
57 const s = bw + 1;
58 if (x !== b) x.set(b);
59 // Forward: L y = b
60 for (let i = 0; i < n; i++) {
61 const rowI = i * s;
62 const kmin = Math.max(0, i - bw);
63 let acc = x[i];
64 for (let k = kmin; k < i; k++) acc -= L[rowI + (i - k)] * x[k];
65 x[i] = acc;
66 }
67 for (let i = 0; i < n; i++) x[i] /= D[i];
68 // Backward: L^T x = y
69 for (let i = n - 1; i >= 0; i--) {
70 const xi = x[i];
71 const rowI = i * s;
72 const kmin = Math.max(0, i - bw);
73 for (let k = kmin; k < i; k++) x[k] -= L[rowI + (i - k)] * xi;
74 }
75}
77/** Number of eigenvalues of (K, M) below the shift. */
78export function inertia(band: Float64Array, mass: Float64Array, n: number, bw: number, shift: number): number {
79 return factorize(band, mass, n, bw, shift).negative;
80}