ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 1// Banded symmetric positive-definite linear algebra for the static condensation: a reverse
2// Cuthill-McKee ordering of the node graph to make the bandwidth small, a banded Cholesky
3// factorisation, and forward/back substitution for many right-hand sides at once.
5/** Reverse Cuthill-McKee ordering of an undirected graph given as CSR adjacency (no self loops
6 * required). Returns perm with perm[newIndex] = oldIndex. */
7export function reverseCuthillMcKee(n: number, ptr: Uint32Array, idx: Uint32Array): Uint32Array {
8 const deg = new Uint32Array(n);
9 for (let i = 0; i < n; i++) deg[i] = ptr[i + 1] - ptr[i];
10 const visited = new Uint8Array(n);
11 const order: number[] = [];
12 const byDegree = Array.from({ length: n }, (_, i) => i).sort((a, b) => deg[a] - deg[b]);
13 for (const start of byDegree) {
14 if (visited[start]) continue;
15 // pseudo-peripheral start: walk to the last node of a BFS level structure twice
16 let s = start;
17 for (let rep = 0; rep < 2; rep++) {
18 const seen = new Uint8Array(n);
19 let frontier = [s];
20 seen[s] = 1;
21 let last = s;
22 while (frontier.length) {
23 const next: number[] = [];
24 for (const u of frontier) {
25 for (let k = ptr[u]; k < ptr[u + 1]; k++) {
26 const v = idx[k];
27 if (!seen[v] && !visited[v]) { seen[v] = 1; next.push(v); }
28 }
29 }
30 if (next.length) { next.sort((a, b) => deg[a] - deg[b]); last = next[0]; }
31 frontier = next;
32 }
33 s = last;
34 }
35 const queue: number[] = [s];
36 visited[s] = 1;
37 let head = 0;
38 while (head < queue.length) {
39 const u = queue[head++];
40 order.push(u);
41 const nb: number[] = [];
42 for (let k = ptr[u]; k < ptr[u + 1]; k++) {
43 const v = idx[k];
44 if (!visited[v]) { visited[v] = 1; nb.push(v); }
45 }
46 nb.sort((a, b) => deg[a] - deg[b]);
47 for (const v of nb) queue.push(v);
48 }
49 }
50 order.reverse();
51 return Uint32Array.from(order);
52}
54/** Symmetric band matrix, lower part stored row-wise: entry (i, j) with i - b <= j <= i lives at
55 * data[i * (b + 1) + (j - i + b)]. */
56export interface BandMatrix {
57 n: number;
58 b: number;
59 data: Float64Array;
60}
62export function bandMatrix(n: number, b: number): BandMatrix {
63 return { n, b, data: new Float64Array(n * (b + 1)) };
64}
66export function bandAdd(A: BandMatrix, i: number, j: number, v: number): void {
67 if (j > i) { const t = i; i = j; j = t; }
68 const off = j - i + A.b;
69 if (off < 0) throw new Error(`entry (${i}, ${j}) outside the band ${A.b}`);
70 A.data[i * (A.b + 1) + off] += v;
71}
73/** In-place Cholesky A = L L^T of a banded SPD matrix; returns false if a pivot is not positive. */
74export function bandedCholesky(A: BandMatrix, onProgress?: (frac: number) => void): boolean {
75 const { n, b, data } = A;
76 const w = b + 1;
77 for (let i = 0; i < n; i++) {
78 const ri = i * w;
79 const j0 = Math.max(0, i - b);
80 for (let j = j0; j <= i; j++) {
81 const rj = j * w;
82 let s = data[ri + (j - i + b)];
83 const k0 = Math.max(j0, j - b);
84 // s -= sum_k L[i,k] L[j,k] for k in [k0, j)
85 let pi = ri + (k0 - i + b), pj = rj + (k0 - j + b);
86 for (let k = k0; k < j; k++) s -= data[pi++] * data[pj++];
87 if (j === i) {
88 if (!(s > 0)) return false;
89 data[ri + b] = Math.sqrt(s);
90 } else {
91 data[ri + (j - i + b)] = s / data[rj + b];
92 }
93 }
94 if (onProgress && (i & 255) === 0) onProgress(i / n);
95 }
96 return true;
97}
99/** Solve L L^T X = R for X, with R given row-major (n x m) and overwritten by X. */
100export function bandedSolve(L: BandMatrix, R: Float64Array, m: number, onProgress?: (frac: number) => void): void {
101 const { n, b, data } = L;
102 const w = b + 1;
103 // forward: L y = r
104 for (let i = 0; i < n; i++) {
105 const ri = i * w, rowI = i * m;
106 const k0 = Math.max(0, i - b);
107 for (let k = k0; k < i; k++) {
108 const l = data[ri + (k - i + b)];
109 if (l === 0) continue;
110 const rowK = k * m;
111 for (let c = 0; c < m; c++) R[rowI + c] -= l * R[rowK + c];
112 }
113 const d = 1 / data[ri + b];
114 for (let c = 0; c < m; c++) R[rowI + c] *= d;
115 if (onProgress && (i & 255) === 0) onProgress(0.5 * i / n);
116 }
117 // backward: L^T x = y (column i of L^T is row i of L)
118 for (let i = n - 1; i >= 0; i--) {
119 const ri = i * w, rowI = i * m;
120 const d = 1 / data[ri + b];
121 for (let c = 0; c < m; c++) R[rowI + c] *= d;
122 const k0 = Math.max(0, i - b);
123 for (let k = k0; k < i; k++) {
124 const l = data[ri + (k - i + b)];
125 if (l === 0) continue;
126 const rowK = k * m;
127 for (let c = 0; c < m; c++) R[rowK + c] -= l * R[rowI + c];
128 }
129 if (onProgress && (i & 255) === 0) onProgress(0.5 + 0.5 * (n - i) / n);
130 }
131}