From 1211a039d4b9e1f901193ac6c31b37d7fcba22d5 Mon Sep 17 00:00:00 2001 From: Meisam Date: Thu, 30 Jul 2026 11:00:30 +0200 Subject: [PATCH] pa-202 --- src/maths/dense/dense.c | 84 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 74 insertions(+), 10 deletions(-) diff --git a/src/maths/dense/dense.c b/src/maths/dense/dense.c index 874812569..bde3cbdb5 100644 --- a/src/maths/dense/dense.c +++ b/src/maths/dense/dense.c @@ -801,22 +801,86 @@ Mat* inverse(Mat* A) { return C; } -CMat* cinverse(CMat* A) { - CMat* B = cadjoint(A); - cplx de = cinv(cdet(A)); +/* Complex matrix inverse by Gauss-Jordan elimination with partial pivoting, in + * O(n^3). This replaces the former adjugate/determinant (Cramer's-rule) method, + * whose cdet() is a recursive cofactor expansion -- O(n!) -- and cadjoint() does + * n^2 of those, so cinverse was O(n*n!). Because the .sp S-parameter analysis + * inverts n x n port matrices at every frequency (CKTspCalcSMatrix, for S, Y and + * Z), that made S-parameter extraction explode by a factor of ~n per added port + * (an 8-port took ~13 s, a 10-port ~18 min); with this it is O(n^3) per inverse. + * Writes cinv_gj() into `out` (n x n, preallocated); returns 0 on ok, 1 if the + * matrix is singular. */ +static int cinv_gj(CMat* A, CMat* out) { + int n = A->row, i, j, c, p; + if (A->col != n) return 1; + /* augmented work: a = copy(A) reduced to I, b = I reduced to A^-1 */ + double* ar = (double*)malloc((size_t)n * n * sizeof(double)); + double* ai = (double*)malloc((size_t)n * n * sizeof(double)); + double* br = (double*)malloc((size_t)n * n * sizeof(double)); + double* bi = (double*)malloc((size_t)n * n * sizeof(double)); + if (!ar || !ai || !br || !bi) { free(ar); free(ai); free(br); free(bi); return 1; } + for (i = 0; i < n; i++) for (j = 0; j < n; j++) { + ar[i*n+j] = A->d[i][j].re; ai[i*n+j] = A->d[i][j].im; + br[i*n+j] = (i == j) ? 1.0 : 0.0; bi[i*n+j] = 0.0; + } + for (c = 0; c < n; c++) { + /* partial pivot on largest |a[.][c]| in rows >= c */ + p = c; double best = ar[c*n+c]*ar[c*n+c] + ai[c*n+c]*ai[c*n+c]; + for (i = c+1; i < n; i++) { + double v = ar[i*n+c]*ar[i*n+c] + ai[i*n+c]*ai[i*n+c]; + if (v > best) { best = v; p = i; } + } + if (best == 0.0) { free(ar); free(ai); free(br); free(bi); return 1; } + if (p != c) for (j = 0; j < n; j++) { + double t; + t = ar[c*n+j]; ar[c*n+j] = ar[p*n+j]; ar[p*n+j] = t; + t = ai[c*n+j]; ai[c*n+j] = ai[p*n+j]; ai[p*n+j] = t; + t = br[c*n+j]; br[c*n+j] = br[p*n+j]; br[p*n+j] = t; + t = bi[c*n+j]; bi[c*n+j] = bi[p*n+j]; bi[p*n+j] = t; + } + /* scale pivot row by 1/pivot */ + double pr = ar[c*n+c], pi = ai[c*n+c], den = pr*pr + pi*pi; + double sr = pr/den, sinv = -pi/den; /* 1/pivot = (pr - i pi)/|p|^2 */ + for (j = 0; j < n; j++) { + double xr, xi; + xr = ar[c*n+j]*sr - ai[c*n+j]*sinv; xi = ar[c*n+j]*sinv + ai[c*n+j]*sr; + ar[c*n+j] = xr; ai[c*n+j] = xi; + xr = br[c*n+j]*sr - bi[c*n+j]*sinv; xi = br[c*n+j]*sinv + bi[c*n+j]*sr; + br[c*n+j] = xr; bi[c*n+j] = xi; + } + /* eliminate column c from all other rows */ + for (i = 0; i < n; i++) { + if (i == c) continue; + double fr = ar[i*n+c], fi = ai[i*n+c]; + if (fr == 0.0 && fi == 0.0) continue; + for (j = 0; j < n; j++) { + ar[i*n+j] -= fr*ar[c*n+j] - fi*ai[c*n+j]; + ai[i*n+j] -= fr*ai[c*n+j] + fi*ar[c*n+j]; + br[i*n+j] -= fr*br[c*n+j] - fi*bi[c*n+j]; + bi[i*n+j] -= fr*bi[c*n+j] + fi*br[c*n+j]; + } + } + } + for (i = 0; i < n; i++) for (j = 0; j < n; j++) { + out->d[i][j].re = br[i*n+j]; out->d[i][j].im = bi[i*n+j]; + } + free(ar); free(ai); free(br); free(bi); + return 0; +} - CMat* C = complexmultiply(B, de); - freecmat(B); +CMat* cinverse(CMat* A) { + CMat* C = newcmatnoinit(A->row, A->col); + if (C == NULL) return NULL; /* true allocation failure (E_NOMEM) */ + if (cinv_gj(A, C)) { /* singular: zero-fill (old contract never returned NULL) */ + for (int i = 0; i < C->row; i++) + for (int j = 0; j < C->col; j++) { C->d[i][j].re = 0.0; C->d[i][j].im = 0.0; } + } return C; } void cinversedest(CMat* A, CMat* dest) { - CMat* B = cadjoint(A); - cplx de = cinv(cdet(A)); - - complexmultiplydest(B, de, dest); - freecmat(B); + cinv_gj(A, dest); return; }