1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 1// Small dense linear algebra: symmetric eigenvalues by cyclic Jacobi, and a least-squares plane fit.
3/** Eigenvalues (ascending) of a symmetric n x n matrix given row-major; the input is copied. */
4export function symmetricEigenvalues(A: Float64Array | number[], n: number): Float64Array {
5 const a = Float64Array.from(A);
6 for (let sweep = 0; sweep < 60; sweep++) {
7 let off = 0;
8 for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) off += a[p * n + q] * a[p * n + q];
9 let diag = 0;
10 for (let p = 0; p < n; p++) diag += a[p * n + p] * a[p * n + p];
11 if (off <= 1e-30 * (diag + 1e-300)) break;
12 for (let p = 0; p < n; p++) {
13 for (let q = p + 1; q < n; q++) {
14 const apq = a[p * n + q];
15 if (apq === 0) continue;
16 const app = a[p * n + p], aqq = a[q * n + q];
17 const theta = (aqq - app) / (2 * apq);
18 const t = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(theta * theta + 1));
19 const c = 1 / Math.sqrt(t * t + 1), s = t * c;
20 for (let k = 0; k < n; k++) {
21 const akp = a[k * n + p], akq = a[k * n + q];
22 a[k * n + p] = c * akp - s * akq; a[k * n + q] = s * akp + c * akq;
23 }
24 for (let k = 0; k < n; k++) {
25 const apk = a[p * n + k], aqk = a[q * n + k];
26 a[p * n + k] = c * apk - s * aqk; a[q * n + k] = s * apk + c * aqk;
27 }
28 }
29 }
30 }
31 const ev = new Float64Array(n);
32 for (let i = 0; i < n; i++) ev[i] = a[i * n + i];
33 return ev.sort();
34}
76fcbe0GPU: rigid modes evaluated from node positions; convert older states on importJeremy Magland 36/** Inverse of an n x n matrix (row-major) by Gauss-Jordan elimination with partial pivoting. */
37export function invertMatrix(A: ArrayLike<number>, n: number): Float64Array {
38 const a = Array.from({ length: n }, (_, i) => [...Array.from({ length: n }, (_, j) => A[n * i + j]), ...Array.from({ length: n }, (_, j) => (i === j ? 1 : 0))]);
39 for (let c = 0; c < n; c++) {
40 let piv = c;
41 for (let r = c + 1; r < n; r++) if (Math.abs(a[r][c]) > Math.abs(a[piv][c])) piv = r;
42 [a[c], a[piv]] = [a[piv], a[c]];
43 const d = a[c][c];
44 for (let j = 0; j < 2 * n; j++) a[c][j] /= d;
45 for (let r = 0; r < n; r++) if (r !== c) { const f = a[r][c]; for (let j = 0; j < 2 * n; j++) a[r][j] -= f * a[c][j]; }
46 }
47 return Float64Array.from(a.flatMap((row) => row.slice(n)));
48}
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 50/** Solve the 3 x 3 system A x = b (row-major A) by Gaussian elimination with pivoting. */
51export function solve3(A: number[], b: number[]): number[] {
52 const M = [A.slice(0, 3).concat(b[0]), A.slice(3, 6).concat(b[1]), A.slice(6, 9).concat(b[2])];
53 for (let c = 0; c < 3; c++) {
54 let piv = c;
55 for (let r = c + 1; r < 3; r++) if (Math.abs(M[r][c]) > Math.abs(M[piv][c])) piv = r;
56 [M[c], M[piv]] = [M[piv], M[c]];
57 if (Math.abs(M[c][c]) < 1e-300) throw new Error("singular 3x3 system");
58 for (let r = 0; r < 3; r++) {
59 if (r === c) continue;
60 const f = M[r][c] / M[c][c];
61 for (let k = c; k < 4; k++) M[r][k] -= f * M[c][k];
62 }
63 }
64 return [M[0][3] / M[0][0], M[1][3] / M[1][1], M[2][3] / M[2][2]];
65}
6d7f34aRank-safe plane fit for contact patchesJeremy Magland 67/** Least-squares fit z = a x + b y + c through points (x_i, y_i, z_i); returns the residuals z - fit.
68 *
69 * The fit is done by Gram-Schmidt on the columns (1, x, y), dropping a column that is (nearly)
70 * dependent on the previous ones, which is what numpy's minimum-norm least squares does in the
71 * Python package: the collinear points of a stick's line get a line fit, a single point a zero
72 * residual. The normal equations would be singular in those cases. */
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 73export function planeResiduals(x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number>): Float64Array {
74 const n = x.length;
6d7f34aRank-safe plane fit for contact patchesJeremy Magland 75 const res = Float64Array.from(z as number[]);
76 const basis: Float64Array[] = [];
77 for (const col of [Array.from({ length: n }, () => 1), Array.from(x), Array.from(y)]) {
78 const v = Float64Array.from(col);
79 let norm0 = 0;
80 for (let i = 0; i < n; i++) norm0 += v[i] * v[i];
81 for (const b of basis) {
82 let d = 0;
83 for (let i = 0; i < n; i++) d += b[i] * v[i];
84 for (let i = 0; i < n; i++) v[i] -= d * b[i];
85 }
86 let norm = 0;
87 for (let i = 0; i < n; i++) norm += v[i] * v[i];
88 if (norm === 0 || norm <= 1e-20 * norm0) continue;
89 const inv = 1 / Math.sqrt(norm);
90 for (let i = 0; i < n; i++) v[i] *= inv;
91 basis.push(v);
92 }
93 for (const b of basis) {
94 let d = 0;
95 for (let i = 0; i < n; i++) d += b[i] * res[i];
96 for (let i = 0; i < n; i++) res[i] -= d * b[i];
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 97 }
98 return res;
99}