ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 1// Dense symmetric eigensolver: Householder reduction to tridiagonal form followed by the
2// implicit QL algorithm with accumulation of the transformations (the tred2 / tql2 pair of
3// EISPACK, in the form given by Numerical Recipes), rearranged so that every inner loop runs
4// along a contiguous row. The eigenvectors are returned as the rows of vectorsT.
6const EPS = Number.EPSILON;
8function pythag(a: number, b: number): number {
9 return Math.hypot(a, b);
10}
12/** Householder tridiagonalisation of the symmetric row-major n x n matrix a (destroyed). On
13 * return d holds the diagonal, e the sub-diagonal (e[0] = 0) and a the orthogonal matrix Q. */
14function tred2(a: Float64Array, n: number, d: Float64Array, e: Float64Array, onProgress?: (frac: number) => void): void {
15 const g = new Float64Array(n);
16 for (let i = n - 1; i >= 1; i--) {
17 const l = i - 1;
18 const ri = i * n;
19 let h = 0, scale = 0;
20 if (l > 0) {
21 for (let k = 0; k <= l; k++) scale += Math.abs(a[ri + k]);
22 if (scale === 0) {
23 e[i] = a[ri + l];
24 } else {
25 for (let k = 0; k <= l; k++) {
26 a[ri + k] /= scale;
27 h += a[ri + k] * a[ri + k];
28 }
29 let f = a[ri + l];
30 const gg = f >= 0 ? -Math.sqrt(h) : Math.sqrt(h);
31 e[i] = scale * gg;
32 h -= f * gg;
33 a[ri + l] = f - gg;
34 f = 0;
35 // p = A u / h, using full rows of the (symmetric) leading block
36 for (let j = 0; j <= l; j++) {
37 a[j * n + i] = a[ri + j] / h;
38 let s = 0;
39 const rj = j * n;
40 for (let k = 0; k <= l; k++) s += a[rj + k] * a[ri + k];
41 e[j] = s / h;
42 f += e[j] * a[ri + j];
43 }
44 const hh = f / (h + h);
45 for (let j = 0; j <= l; j++) e[j] -= hh * a[ri + j];
46 // rank-two update of the full leading block: A -= u q^T + q u^T
47 for (let j = 0; j <= l; j++) {
48 const fj = a[ri + j], gj = e[j];
49 const rj = j * n;
50 for (let k = 0; k <= l; k++) a[rj + k] -= fj * e[k] + gj * a[ri + k];
51 }
52 }
53 } else {
54 e[i] = a[ri + l];
55 }
56 d[i] = h;
57 if (onProgress && (i & 63) === 0) onProgress((n - i) / n);
58 }
59 d[0] = 0;
60 e[0] = 0;
61 // accumulate the transformations into a (Q)
62 for (let i = 0; i < n; i++) {
63 const l = i - 1;
64 const ri = i * n;
65 if (d[i] !== 0) {
66 g.fill(0, 0, l + 1);
67 for (let k = 0; k <= l; k++) {
68 const aik = a[ri + k];
69 if (aik === 0) continue;
70 const rk = k * n;
71 for (let j = 0; j <= l; j++) g[j] += aik * a[rk + j];
72 }
73 for (let k = 0; k <= l; k++) {
74 const aki = a[k * n + i];
75 if (aki === 0) continue;
76 const rk = k * n;
77 for (let j = 0; j <= l; j++) a[rk + j] -= g[j] * aki;
78 }
79 }
80 d[i] = a[ri + i];
81 a[ri + i] = 1;
82 for (let j = 0; j <= l; j++) {
83 a[j * n + i] = 0;
84 a[ri + j] = 0;
85 }
86 }
87}
89/** Implicit QL on the tridiagonal (d, e), with the rotations applied to the rows of zT. */
90function tql2(d: Float64Array, e: Float64Array, n: number, zT: Float64Array, onProgress?: (frac: number) => void): void {
91 for (let i = 1; i < n; i++) e[i - 1] = e[i];
92 e[n - 1] = 0;
93 for (let l = 0; l < n; l++) {
94 let iter = 0;
95 let m: number;
96 do {
97 for (m = l; m < n - 1; m++) {
98 const dd = Math.abs(d[m]) + Math.abs(d[m + 1]);
99 if (Math.abs(e[m]) <= EPS * dd) break;
100 }
101 if (m !== l) {
102 if (iter++ === 100) throw new Error("tql2: too many iterations");
103 let g = (d[l + 1] - d[l]) / (2 * e[l]);
104 let r = pythag(g, 1);
105 g = d[m] - d[l] + e[l] / (g + (g >= 0 ? Math.abs(r) : -Math.abs(r)));
106 let s = 1, c = 1, p = 0;
107 let i: number;
108 for (i = m - 1; i >= l; i--) {
109 const f = s * e[i], b = c * e[i];
110 e[i + 1] = r = pythag(f, g);
111 if (r === 0) {
112 d[i + 1] -= p;
113 e[m] = 0;
114 break;
115 }
116 s = f / r;
117 c = g / r;
118 g = d[i + 1] - p;
119 r = (d[i] - g) * s + 2 * c * b;
120 p = s * r;
121 d[i + 1] = g + p;
122 g = c * r - b;
123 const ri = i * n, rj = (i + 1) * n;
124 for (let k = 0; k < n; k++) {
125 const f2 = zT[rj + k], z1 = zT[ri + k];
126 zT[rj + k] = s * z1 + c * f2;
127 zT[ri + k] = c * z1 - s * f2;
128 }
129 }
130 if (r === 0 && i >= l) continue;
131 d[l] -= p;
132 e[l] = g;
133 e[m] = 0;
134 }
135 } while (m !== l);
136 if (onProgress && (l & 15) === 0) onProgress(l / n);
137 }
138}
140export interface EigenResult {
141 values: Float64Array; // ascending
142 vectorsT: Float64Array; // (n, n): row k is the eigenvector of values[k]
143}
145/** All eigenpairs of the symmetric row-major matrix A (n x n, destroyed). */
146export function symmetricEigen(A: Float64Array, n: number, onProgress?: (phase: string, frac: number) => void): EigenResult {
147 const d = new Float64Array(n), e = new Float64Array(n);
148 tred2(A, n, d, e, onProgress ? (f) => onProgress("tridiagonalisation", f) : undefined);
149 // Q is in A (column j = j-th transformation); tql2 wants the transposed layout
150 const zT = new Float64Array(n * n);
151 for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) zT[j * n + i] = A[i * n + j];
152 tql2(d, e, n, zT, onProgress ? (f) => onProgress("QL iteration", f) : undefined);
153 // sort ascending
154 const order = Array.from({ length: n }, (_, i) => i).sort((p, q) => d[p] - d[q]);
155 const values = new Float64Array(n), vectorsT = new Float64Array(n * n);
156 order.forEach((k, i) => {
157 values[i] = d[k];
158 vectorsT.set(zT.subarray(k * n, (k + 1) * n), i * n);
159 });
160 return { values, vectorsT };
161}