diff --git a/examples/sp/file.cir b/examples/sp/file.cir index d34491bf0..7a047bea8 100644 --- a/examples/sp/file.cir +++ b/examples/sp/file.cir @@ -30,22 +30,33 @@ ri net1 sin 50 X1 sin sp_model_out 0 s2p_generic touchstone="137MHz_BPF.s2p" ro sp_model_out 0 50 +* A version derived from extracted s-parameters using the nport + +VP4 net2 0 dc 0 ac 0.316228 ; Same as VP1 but that must not be loaded. +ri2 net2 sin2 50 +N1 sin2 np_model_out 0 mymodel +.model mymodel nport(file="new_137MHz_BPF.nport") +ro2 np_model_out 0 50 + * Control section, run simulations. .control -* Run an SP analysis, derive group delay and write a Touchstone file - +* Run an SP analysis, derive group delay, write a Touchstone file +* Use this Touchstone file for vector fitting by pre_snp, thus get a file +* new_137MHz_BPF.nport SP LIN 71 107MEG 167MEG let Group_Delay = -1*deriv(cph(S_2_1))/(2*pi) wrs2p new_137MHz_BPF.s2p +pre_snp -native new_137MHz_BPF.s2p * Now run AC analysis alter vp2 ac = 0 AC LIN 71 107MEG 167MEG set xbrushwidth=2 -plot group_delay(filter_out) sp1.Group_Delay group_delay(sp_model_out) -plot db(filter_out) cph(filter_out) db(sp_model_out) cph(sp_model_out) +display +plot group_delay(filter_out) sp1.Group_Delay group_delay(sp_model_out) group_delay(np_model_out) +plot db(filter_out) cph(filter_out) db(sp_model_out) cph(sp_model_out) db(np_model_out) cph(np_model_out) .endc * Subcircuit for SP behavioral device diff --git a/examples/sp/file_nport.cir b/examples/sp/file_nport.cir new file mode 100644 index 000000000..3a7514c97 --- /dev/null +++ b/examples/sp/file_nport.cir @@ -0,0 +1,52 @@ +Example circuit to read a Touchstone-format parameter file. + +* This file contains an an analog filter circuit; SP and AC analyses are run +* and the S-parameters are saved in a file. +* It also uses vetor fitting for a *.nport file serving the nport instance +* in a circuit that reproduces the original from its saved S-parameters. +* Results from the two versions are plotted together for comparison. + +* Derived from: Qucs 1.0.2 C:/Users/Tom Hajjar/.qucs - S/S-parameter_files_ngspice_prj/137MHz_BPF.sch + +* Original circuit for SP analysis + +VP2 filter_out 0 dc 0 ac 0.316228 SIN(0 0.316228 1G) portnum 2 z0 50 +VP1 _net1 0 dc 0 ac 0.316228 SIN(0 0.316228 1G) portnum 1 z0 50 +C1 _net2 _net1 4.7P +C5 _net3 filter_out 4.7P +C3 _net3 _net2 1P +C4 _net3 0 15P +C2 0 _net2 15P +L2 0 _net3 66N +L1 0 _net2 66N + +* A version fitted from extracted s-parameters using the nport + +VP4 net2 0 dc 0 ac 0.316228 ; Same as VP1 but that must not be loaded. +ri2 net2 sin2 50 +N1 sin2 np_model_out 0 mymodel +.model mymodel nport(file="new_137MHz_BPF.nport") +ro2 np_model_out 0 50 + +* Control section, run simulations. + +.control +* Run an SP analysis, derive group delay, write a Touchstone file +* Use this Touchstone file for vector fitting by pre_snp, thus get a file +* new_137MHz_BPF.nport +SP LIN 71 107MEG 167MEG +let Group_Delay = -1*deriv(cph(S_2_1))/(2*pi) +wrs2p new_137MHz_BPF.s2p +pre_snp -native new_137MHz_BPF.s2p + +* Now run AC analysis + +alter vp2 ac = 0 +AC LIN 71 107MEG 167MEG +set xbrushwidth=2 +display +plot group_delay(filter_out) sp1.Group_Delay group_delay(np_model_out) +plot db(filter_out) cph(filter_out) db(np_model_out) cph(np_model_out) +.endc + +.END diff --git a/src/frontend/snp2va.cpp b/src/frontend/snp2va.cpp deleted file mode 100644 index 41c516ebd..000000000 --- a/src/frontend/snp2va.cpp +++ /dev/null @@ -1,988 +0,0 @@ -/* Enhancement-200: Touchstone (.sNp) -> Verilog-A n-port converter, in C. - * - * A C port of the pure-Python snp2va.py (Enhancement-199): parse a Touchstone - * S-parameter file, convert S -> Y, fit every Y_ij(f) with a common-pole rational - * (Gustavsen vector fitting), and emit a Verilog-A n-port realized with laplace_nd, - * so through OpenVAF/OSDI it works in AC and transient. Used by the `pre_snp` - * front-end command, which then invokes openvaf-r to compile the emitted .va. - * - * The numerical core is self-contained (only stdio/stdlib/string/math/complex), - * so it does not pull in ngspice's own complex.h. Public entry point: - * int snp2va_convert(const char *snp, const char *va, const char *module, - * char *msg, int msglen); - * returns 0 on success, non-zero on failure (msg gets a one-line status). - */ -#ifdef _MSC_VER -#include -#include - -extern "C" { -#include "snp2va.h" -#include -}; -#else -#include -#include -#include -#include -#include -#include -#include "snp2va.h" -#endif - -#ifdef _MSC_VER -typedef std::complex cplx; -#define cabs abs -#define cpow pow -#define cimag imag -#define creal real -#define cexp exp -#define I cplx(0.0, 1.0) -#define strcasecmp _stricmp -#else -typedef double _Complex cplx; -#endif - -#ifndef M_PI -#define M_PI 3.14159265358979323846 -#endif - -/* ============================ linear algebra ============================ */ - -/* Least-squares min||A x - b|| for a REAL system via Householder QR. A is - * row-major m*n, b length m, x length n. Returns 0 on ok. Handles the - * underdetermined case m=m) at zero rather than reading past the m-row A and b - * (a heap-buffer-overflow reachable from `pre_snp` on a Touchstone file with very - * few frequency points, where the vector fit's stacked system has fewer rows than - * poles). For the normal overdetermined path (m>=n) the guard never fires. */ -static int lstsq_real(double *A, double *b, int m, int n, double *x) -{ - int i, j, k; - for (k = 0; k < n; k++) { - double norm = 0.0; - for (i = k; i < m; i++) norm += A[i*n+k]*A[i*n+k]; - norm = sqrt(norm); - if (norm == 0.0) continue; - double alpha = (A[k*n+k] >= 0.0) ? -norm : norm; - double *v = (double*) calloc((size_t) m, sizeof(double)); - v[k] = A[k*n+k] - alpha; - for (i = k+1; i < m; i++) v[i] = A[i*n+k]; - double vn2 = 0.0; - for (i = k; i < m; i++) vn2 += v[i]*v[i]; - if (vn2 == 0.0) { free(v); continue; } - for (j = k; j < n; j++) { - double s = 0.0; - for (i = k; i < m; i++) s += v[i]*A[i*n+j]; - s = s*2.0/vn2; - for (i = k; i < m; i++) A[i*n+j] -= s*v[i]; - } - double s = 0.0; - for (i = k; i < m; i++) s += v[i]*b[i]; - s = s*2.0/vn2; - for (i = k; i < m; i++) b[i] -= s*v[i]; - free(v); - } - for (i = n-1; i >= 0; i--) { - if (i >= m) { x[i] = 0.0; continue; } /* unknown i has no constraining row */ - double acc = b[i]; - for (j = i+1; j < n; j++) acc -= A[i*n+j]*x[j]; - x[i] = (A[i*n+i] != 0.0) ? acc/A[i*n+i] : 0.0; - } - return 0; -} - -/* In-place inverse of an n x n complex matrix (row-major), Gauss-Jordan w/ pivot. - * Returns 0 on ok, 1 if singular. */ -static int mat_inv_c(cplx *M, int n) -{ - int i, j, c, p; - cplx *A = (cplx*) malloc((size_t) n*2*n*sizeof(cplx)); - for (i = 0; i < n; i++) { - for (j = 0; j < n; j++) A[i*2*n+j] = M[i*n+j]; - for (j = 0; j < n; j++) A[i*2*n+n+j] = (i==j) ? 1.0 : 0.0; - } - for (c = 0; c < n; c++) { - p = c; double best = cabs(A[c*2*n+c]); - for (i = c+1; i < n; i++) { double v = cabs(A[i*2*n+c]); if (v > best) { best = v; p = i; } } - if (best == 0.0) { free(A); return 1; } - if (p != c) for (j = 0; j < 2*n; j++) { cplx t = A[c*2*n+j]; A[c*2*n+j]=A[p*2*n+j]; A[p*2*n+j]=t; } - cplx d = A[c*2*n+c]; - for (j = 0; j < 2*n; j++) A[c*2*n+j] /= d; - for (i = 0; i < n; i++) if (i != c) { - cplx f = A[i*2*n+c]; - if (f != 0.0) for (j = 0; j < 2*n; j++) A[i*2*n+j] -= f*A[c*2*n+j]; - } - } - for (i = 0; i < n; i++) for (j = 0; j < n; j++) M[i*n+j] = A[i*2*n+n+j]; - free(A); - return 0; -} - -/* Monic polynomial (DESCENDING, len nr+1) from roots. coef must hold nr+1. */ -static void poly_from_roots(const cplx *roots, int nr, cplx *coef) -{ - int i, k; - coef[0] = 1.0; - for (i = 1; i <= nr; i++) coef[i] = 0.0; - for (k = 0; k < nr; k++) { - for (i = k+1; i >= 1; i--) coef[i] = coef[i] - roots[k]*coef[i-1]; - } -} - -static cplx poly_eval(const cplx *c, int deg, cplx x) -{ - cplx r = 0.0; int i; - for (i = 0; i <= deg; i++) r = r*x + c[i]; - return r; -} - -/* All roots of a DESCENDING-coeff polynomial (deg = len-1) via Durand-Kerner. - * roots[] must hold deg. Returns 0 on ok. */ -static int poly_roots(const cplx *cin, int deg, cplx *roots) -{ - int i, j, it; - if (deg <= 0) return 0; - cplx *c = (cplx*) malloc((size_t)(deg+1)*sizeof(cplx)); - for (i = 0; i <= deg; i++) c[i] = cin[i] / cin[0]; /* monic */ - cplx seed = 0.4 + 0.9*I; - for (i = 0; i < deg; i++) roots[i] = cpow(seed, (double) i); - for (it = 0; it < 500; it++) { - double maxstep = 0.0; - for (i = 0; i < deg; i++) { - cplx num = poly_eval(c, deg, roots[i]); - cplx den = 1.0; - for (j = 0; j < deg; j++) if (j != i) den *= (roots[i]-roots[j]); - cplx step = (cabs(den) > 1e-300) ? num/den : 0.0; - roots[i] -= step; - if (cabs(step) > maxstep) maxstep = cabs(step); - } - if (maxstep < 1e-14) break; - } - free(c); - return 0; -} - -/* ============================ Touchstone I/O ============================ */ - -typedef struct { double *freqs; cplx *S; int nf; int N; double z0; char ptype; } TS; - -static void ts_free(TS *t) { free(t->freqs); free(t->S); } - -/* returns 0 on ok */ -static int parse_touchstone(const char *fn, TS *out, char *msg, int msglen) -{ - FILE *f = fopen(fn, "r"); - if (!f) { snprintf(msg, (size_t) msglen, "cannot open '%s'", fn); return 1; } - double fmul = 1e9, z0 = 50.0; char ptype = 'S'; char fmt[3] = "MA"; - /* collect all numeric tokens after the '#' options line(s) */ - double *nums = NULL; long ncap = 0, nn = 0; - char line[4096]; - while (fgets(line, sizeof line, f)) { - char *h = strchr(line, '!'); if (h) *h = '\0'; - char *p = line; - while (*p && isspace((unsigned char)*p)) p++; - if (*p == '\0') continue; - if (*p == '#') { - char *tok = strtok(p+1, " \t\r\n"); - while (tok) { - if (!strcasecmp(tok,"HZ")) fmul=1.0; - else if (!strcasecmp(tok,"KHZ")) fmul=1e3; - else if (!strcasecmp(tok,"MHZ")) fmul=1e6; - else if (!strcasecmp(tok,"GHZ")) fmul=1e9; - else if (!strcasecmp(tok,"S")||!strcasecmp(tok,"Y")||!strcasecmp(tok,"Z")) ptype=(char)toupper((unsigned char)tok[0]); - else if (!strcasecmp(tok,"MA")||!strcasecmp(tok,"DB")||!strcasecmp(tok,"RI")) { fmt[0]=(char)toupper((unsigned char)tok[0]); fmt[1]=(char)toupper((unsigned char)tok[1]); } - else if (!strcasecmp(tok,"R")) { char *z=strtok(NULL," \t\r\n"); if (z) z0=atof(z); } - tok = strtok(NULL, " \t\r\n"); - } - continue; - } - /* numeric data line */ - char *tok = strtok(p, " \t\r\n"); - while (tok) { - char *end; double v = strtod(tok, &end); - if (end != tok) { - if (nn >= ncap) { ncap = ncap ? ncap*2 : 1024; nums = (double*) realloc(nums, (size_t) ncap*sizeof(double)); } - nums[nn++] = v; - } - tok = strtok(NULL, " \t\r\n"); - } - } - fclose(f); - /* infer port count N: try the file extension .sNp, else brute force */ - int N = 0; - const char *dot = strrchr(fn, '.'); - if (dot && (dot[1]=='s'||dot[1]=='S') && (fn[strlen(fn)-1]=='p'||fn[strlen(fn)-1]=='P')) { - N = atoi(dot+2); - /* Enhancement-227: reject an implausible port count from the filename - * (e.g. `.s2147483647p`). N is stored in out->N and used to size the - * downstream N x N vector fit; a huge N over-allocates / overflows and - * corrupts the heap. Real Touchstone files have few ports -- above the - * brute-force limit, drop back to inferring N from the data. */ - if (N > 512) - N = 0; - } - if (N <= 0) { - int c; - for (c = 1; c <= 512; c++) if (nn % (1 + 2*c*c) == 0) { N = c; break; } - } - if (N <= 0) { free(nums); snprintf(msg,(size_t)msglen,"cannot determine port count"); return 1; } - int rec = 1 + 2*N*N; - int nf = (int)(nn / rec); - if (nf < 2) { free(nums); snprintf(msg,(size_t)msglen,"too few frequency points (%d)", nf); return 1; } - out->freqs = (double*) malloc((size_t) nf*sizeof(double)); - out->S = (cplx*) malloc((size_t) nf*N*N*sizeof(cplx)); - out->nf = nf; out->N = N; out->z0 = z0; out->ptype = ptype; - cplx *pv = (cplx*) malloc((size_t) N*N*sizeof(cplx)); /* heap: N may be large */ - int r, kk; - for (r = 0; r < nf; r++) { - double *chunk = nums + (long) r*rec; - out->freqs[r] = chunk[0]*fmul; - double *vals = chunk+1; - for (kk = 0; kk < N*N; kk++) { - double a = vals[2*kk], b = vals[2*kk+1]; - if (!strcmp(fmt,"MA")) pv[kk] = a*cexp(I*b*M_PI/180.0); - else if (!strcmp(fmt,"DB")) pv[kk] = pow(10.0,a/20.0)*cexp(I*b*M_PI/180.0); - else pv[kk] = a + I*b; - } - /* Touchstone: N=2 order is S11 S21 S12 S22; general is row-major */ - cplx *M = out->S + (long) r*N*N; - if (N == 2) { M[0]=pv[0]; M[2]=pv[1]; M[1]=pv[2]; M[3]=pv[3]; } - else for (kk = 0; kk < N*N; kk++) M[kk] = pv[kk]; - } - free(pv); free(nums); - return 0; -} - -/* S/Y/Z -> Y (row-major per frequency), Yout must hold nf*N*N. */ -static int to_Y(const TS *t, cplx *Yout) -{ - int N = t->N, r, i, j; - cplx *tmp = (cplx*) malloc((size_t) N*N*sizeof(cplx)); - for (r = 0; r < t->nf; r++) { - const cplx *M = t->S + (long) r*N*N; - cplx *Y = Yout + (long) r*N*N; - if (t->ptype == 'Y') { for (i=0;iptype == 'Z') { for (i=0;i Y = (1/z0)(I-S)(I+S)^-1 */ - cplx *IpS = tmp; - for (i=0;iz0; - } - } - } - free(tmp); - return 0; -} - -/* ============================ vector fitting ============================ */ - -/* complex partial-fraction basis (real-valued cc combos), Ns x Np, row-major */ -static void build_basis(const cplx *s, int Ns, const cplx *poles, int Np, cplx *A) -{ - int r, i; - for (r = 0; r < Ns; r++) { - i = 0; - while (i < Np) { - if (fabs(cimag(poles[i])) < 1e-9*fabs(creal(poles[i]))+1e-30) { - A[r*Np+i] = 1.0/(s[r]-poles[i]); i += 1; - } else { - cplx p = poles[i]; - A[r*Np+i] = 1.0/(s[r]-p) + 1.0/(s[r]-conj(p)); - A[r*Np+i+1] = I/(s[r]-p) - I/(s[r]-conj(p)); - i += 2; - } - } - } -} - -/* complex residues from real ctil coeffs, per real/cc layout */ -static void ctil_to_cres(const double *ctil, const cplx *poles, int Np, cplx *cres) -{ - int i = 0; - while (i < Np) { - if (fabs(cimag(poles[i])) < 1e-9*fabs(creal(poles[i]))+1e-30) { cres[i] = ctil[i]; i += 1; } - else { cres[i] = ctil[i] + I*ctil[i+1]; cres[i+1] = conj(cres[i]); i += 2; } - } -} - -/* Canonicalize a pole set into the layout every consumer here assumes: real poles - * (exact Im=0) first, then complex poles as ADJACENT exact conjugate pairs. The - * Durand-Kerner roots of the (real) sigma numerator are conjugate-symmetric in - * theory, but numerical noise can leave a "pair" split or a near-real pole with a - * tiny nonzero imag -- which makes build_basis / ctil_to_cres walk one slot past - * the array (heap overflow). Rebuilding pairs from the Im>0 representatives (using - * exact conj) guarantees the structure and cleans the asymmetry. */ -static void canon_poles(cplx *p, int Np) -{ - int i, no = 0; - cplx *out = (cplx*) malloc((size_t) Np*sizeof(cplx)); - for (i = 0; i < Np; i++) - if (fabs(cimag(p[i])) < 1e-6*cabs(p[i])) p[i] = creal(p[i]); /* snap near-real */ - for (i = 0; i < Np; i++) - if (cimag(p[i]) == 0.0 && no < Np) out[no++] = p[i]; /* reals first */ - for (i = 0; i < Np; i++) - if (cimag(p[i]) > 0.0) { /* one per pair */ - if (no+1 < Np) { out[no++] = p[i]; out[no++] = conj(p[i]); } - else if (no < Np) out[no++] = creal(p[i]); /* no room -> real */ - } - while (no < Np) { out[no] = -fabs(creal(p[no])) ; no++; } /* pad (safety) */ - for (i = 0; i < Np; i++) p[i] = out[i]; - free(out); -} - -/* Householder-reduce B (m x nB, row-major) in place, applying every reflector to - * the tail T (m x nT) and rhs g (m). After this the rows [nB..m) of T and g are a - * reduced least-squares system in the remaining (shared) unknowns only -- the per- - * element unknowns spanned by B have been projected out. (Fast Vector Fitting: - * Deschrijver et al. 2008 -- avoids ever forming the full block-arrow matrix.) */ -static void hh_reduce(double *B, double *T, double *g, int m, int nB, int nT) -{ - int i, j, k; - double *v = (double*) malloc((size_t) m * sizeof(double)); - for (k = 0; k < nB; k++) { - double norm = 0.0; - for (i = k; i < m; i++) norm += B[i*nB+k]*B[i*nB+k]; - norm = sqrt(norm); - if (norm == 0.0) continue; - double alpha = (B[k*nB+k] >= 0.0) ? -norm : norm; - for (i = 0; i < k; i++) v[i] = 0.0; - v[k] = B[k*nB+k] - alpha; - for (i = k+1; i < m; i++) v[i] = B[i*nB+k]; - double vn2 = 0.0; - for (i = k; i < m; i++) vn2 += v[i]*v[i]; - if (vn2 == 0.0) continue; - for (j = k; j < nB; j++) { - double sdot = 0.0; for (i = k; i < m; i++) sdot += v[i]*B[i*nB+j]; - sdot = sdot*2.0/vn2; for (i = k; i < m; i++) B[i*nB+j] -= sdot*v[i]; - } - for (j = 0; j < nT; j++) { - double sdot = 0.0; for (i = k; i < m; i++) sdot += v[i]*T[i*nT+j]; - sdot = sdot*2.0/vn2; for (i = k; i < m; i++) T[i*nT+j] -= sdot*v[i]; - } - { double sdot = 0.0; for (i = k; i < m; i++) sdot += v[i]*g[i]; - sdot = sdot*2.0/vn2; for (i = k; i < m; i++) g[i] -= sdot*v[i]; } - } - free(v); -} - -/* One vector-fit run (fixed pole count) over the element indices in elems[0..Ne). - * s,F normalized; F is [N*N][Ns], res/d/e are written only for the listed elements. - * Common poles are identified with the FAST (block-reduced) pole solve, so cost is - * O(Ne*Ns*Np^2) with O(Ne*Ns*Np) memory instead of the O(Ne^2) dense stack. Stops - * early once the poles stop moving. */ -static void vector_fit(const cplx *s, int Ns, const cplx *F, - const int *elems, int Ne, int Np, - cplx *poles, cplx *res, double *d, double *e, int maxiter) -{ - int iter, i, j, k, r, ke; - cplx *A = (cplx*) malloc((size_t) Ns*Np*sizeof(cplx)); - int m = 2*Ns, nB = Np+2, nT = Np, redrows = 2*Ns - (Np+2); - if (redrows < 0) redrows = 0; - long stackrows = (long) redrows * Ne; - double *SC = (double*) malloc((size_t) stackrows * Np * sizeof(double)); - double *Sb = (double*) malloc((size_t) stackrows * sizeof(double)); - double *B = (double*) malloc((size_t) m * nB * sizeof(double)); - double *T = (double*) malloc((size_t) m * nT * sizeof(double)); - double *g = (double*) malloc((size_t) m * sizeof(double)); - - for (iter = 0; iter < maxiter; iter++) { - build_basis(s, Ns, poles, Np, A); - long sr = 0; - for (ke = 0; ke < Ne; ke++) { - k = elems[ke]; - for (r = 0; r < Ns; r++) { - cplx Fkr = F[(long)k*Ns + r]; - for (j = 0; j < Np; j++) { - cplx a = A[r*Np+j], neg = -Fkr*a; - B[r*nB + j] = creal(a); B[(r+Ns)*nB + j] = cimag(a); - T[r*nT + j] = creal(neg); T[(r+Ns)*nT + j] = cimag(neg); - } - B[r*nB + Np] = 1.0; B[(r+Ns)*nB + Np] = 0.0; /* d */ - B[r*nB + Np+1] = creal(s[r]); B[(r+Ns)*nB + Np+1] = cimag(s[r]); /* e*s */ - g[r] = creal(Fkr); g[r+Ns] = cimag(Fkr); - } - hh_reduce(B, T, g, m, nB, nT); - for (i = nB; i < m; i++) { - for (j = 0; j < Np; j++) SC[sr*Np + j] = T[i*nT + j]; - Sb[sr] = g[i]; - sr++; - } - } - double *ctil = (double*) calloc((size_t) Np, sizeof(double)); - lstsq_real(SC, Sb, (int) sr, Np, ctil); - cplx *cres = (cplx*) malloc((size_t) Np*sizeof(cplx)); - ctil_to_cres(ctil, poles, Np, cres); - /* relocate: roots of D(s) + sum cres_i * D(s)/(s-a_i) */ - cplx *D = (cplx*) malloc((size_t)(Np+1)*sizeof(cplx)); - poly_from_roots(poles, Np, D); - cplx *numsig = (cplx*) malloc((size_t)(Np+1)*sizeof(cplx)); - for (i = 0; i <= Np; i++) numsig[i] = D[i]; - cplx *Di = (cplx*) malloc((size_t) Np*sizeof(cplx)); - cplx *sub = (cplx*) malloc((size_t) Np*sizeof(cplx)); /* poles minus i */ - for (i = 0; i < Np; i++) { - int t = 0; for (j = 0; j < Np; j++) if (j != i) sub[t++] = poles[j]; - poly_from_roots(sub, Np-1, Di); /* len Np, degree Np-1 */ - for (j = 0; j < Np; j++) numsig[j+1] += cres[i]*Di[j]; - } - cplx *newp = (cplx*) malloc((size_t) Np*sizeof(cplx)); - poly_roots(numsig, Np, newp); - for (i = 0; i < Np; i++) if (creal(newp[i]) > 0) newp[i] = -creal(newp[i]) + I*cimag(newp[i]); - /* sort: real poles first, then by real, imag (keeps cc pairs adjacent-ish) */ - for (i = 0; i < Np; i++) for (j = i+1; j < Np; j++) { - int swap = 0; - double ki = (fabs(cimag(newp[i]))>1e-6)?1:0, kj = (fabs(cimag(newp[j]))>1e-6)?1:0; - if (kj < ki) swap = 1; - else if (kj == ki) { if (creal(newp[j]) < creal(newp[i]) - 1e-30) swap = 1; - else if (fabs(creal(newp[j])-creal(newp[i]))<1e-30 && cimag(newp[j]) convergence */ - for (i = 0; i < Np; i++) { - double dm = cabs(newp[i]-poles[i])/(cabs(poles[i])+1e-300); - if (dm > mv) mv = dm; - } - for (i = 0; i < Np; i++) poles[i] = newp[i]; - free(ctil); free(cres); free(D); free(numsig); free(Di); free(sub); free(newp); - if (mv < 1e-4) break; /* poles settled -- stop early (Fix #3) */ - } - /* final residues (fixed poles), per fitted element */ - build_basis(s, Ns, poles, Np, A); - for (ke = 0; ke < Ne; ke++) { - k = elems[ke]; - int ncol = Np+2, nrow = Ns; - double *M = (double*) calloc((size_t)(2*nrow)*ncol, sizeof(double)); - double *bb = (double*) calloc((size_t)(2*nrow), sizeof(double)); - for (r = 0; r < Ns; r++) { - cplx Fkr = F[(long)k*Ns+r]; - for (j = 0; j < Np; j++) { M[r*ncol+j]=creal(A[r*Np+j]); M[(r+nrow)*ncol+j]=cimag(A[r*Np+j]); } - M[r*ncol+Np]=1.0; - M[r*ncol+Np+1]=creal(s[r]); M[(r+nrow)*ncol+Np+1]=cimag(s[r]); - bb[r]=creal(Fkr); bb[r+nrow]=cimag(Fkr); - } - double *x = (double*) calloc((size_t) ncol, sizeof(double)); - lstsq_real(M, bb, 2*nrow, ncol, x); - cplx *cr = (cplx*) malloc((size_t) Np*sizeof(cplx)); - ctil_to_cres(x, poles, Np, cr); - for (j = 0; j < Np; j++) res[(long)k*Np+j] = cr[j]; - d[k] = x[Np]; e[k] = x[Np+1]; - free(M); free(bb); free(x); free(cr); - } - free(A); free(SC); free(Sb); free(B); free(T); free(g); -} - -/* Reciprocal network? (S/Y symmetric -> fit only the upper triangle.) F is the - * Y data [N*N][Ns]; compare Y_ij vs Y_ji across a frequency subset. */ -static int is_reciprocal(const cplx *F, int N, int Ns) -{ - double maxd = 0.0, maxv = 0.0; - int i, j, r; - int step = Ns/16 > 0 ? Ns/16 : 1; - for (i = 0; i < N; i++) for (j = i+1; j < N; j++) - for (r = 0; r < Ns; r += step) { - cplx a = F[(long)(i*N+j)*Ns+r], b = F[(long)(j*N+i)*Ns+r]; - double dd = cabs(a-b), va = cabs(a); - if (dd > maxd) maxd = dd; - if (va > maxv) maxv = va; - } - return maxd <= 1e-6*(maxv+1e-300); -} - -static void seed_poles(double fmin, double fmax, int npair, cplx *p) -{ - int k; - for (k = 0; k < npair; k++) { - double beta = (npair==1) ? 2*M_PI*sqrt(fmin*fmax) - : 2*M_PI*fmin*pow(fmax/fmin, (double)k/(npair-1)); - double a = beta/100.0; - p[2*k] = -a + I*beta; - p[2*k+1] = -a - I*beta; - } -} - -/* ============================ emit VA ============================ */ -static void emit_arr(FILE *f, const double *v, int n) -{ - int i; fprintf(f, "'{"); - for (i = 0; i < n; i++) { - char buf[64]; - snprintf(buf, sizeof buf, "%.12g", v[i]); - /* Force a REAL literal: `%.12g` prints 1.0 as "1", and OpenVAF crashes on an - * integer literal inside a laplace_nd coefficient array assigned to a - * variable (index-out-of-bounds in its lowering). Append ".0" if needed. */ - fprintf(f, "%s%s%s", i?", ":"", buf, strpbrk(buf, ".eE") ? "" : ".0"); - } - fprintf(f, "}"); -} - -/* term separator inside one `I(p) <+ ...;` contribution */ -static void term_sep(FILE *f, int *first) -{ - if (*first) *first = 0; - else fprintf(f, "\n + "); -} - -/* Symmetric eigen-decomposition via cyclic Jacobi. A (n x n, row-major) is - * overwritten; eigenvalues -> w, orthonormal eigenvectors (columns) -> V. */ -static void jacobi_sym(double *A, int n, double *w, double *V) -{ - int i, j, p, q, sweep; - for (i = 0; i < n; i++) { for (j = 0; j < n; j++) V[i*n+j] = (i==j)?1.0:0.0; } - for (sweep = 0; sweep < 100; sweep++) { - double off = 0.0; - for (p = 0; p < n; p++) for (q = p+1; q < n; q++) off += A[p*n+q]*A[p*n+q]; - if (off < 1e-300) break; - for (p = 0; p < n; p++) for (q = p+1; q < n; q++) { - double apq = A[p*n+q]; - if (fabs(apq) < 1e-300) continue; - double app = A[p*n+p], aqq = A[q*n+q]; - double phi = 0.5*(aqq-app)/apq; - double t = (phi>=0?1.0:-1.0)/(fabs(phi)+sqrt(phi*phi+1.0)); - double c = 1.0/sqrt(t*t+1.0), sn = t*c; - for (i = 0; i < n; i++) { - double aip = A[i*n+p], aiq = A[i*n+q]; - A[i*n+p] = c*aip - sn*aiq; A[i*n+q] = sn*aip + c*aiq; - } - for (i = 0; i < n; i++) { - double api = A[p*n+i], aqi = A[q*n+i]; - A[p*n+i] = c*api - sn*aqi; A[q*n+i] = sn*api + c*aqi; - } - for (i = 0; i < n; i++) { - double vip = V[i*n+p], viq = V[i*n+q]; - V[i*n+p] = c*vip - sn*viq; V[i*n+q] = sn*vip + c*viq; - } - } - } - for (i = 0; i < n; i++) w[i] = A[i*n+i]; -} - -/* Project the improper (e*s) capacitance matrix onto the symmetric PSD cone. - * Each Y_ij is fit independently, so E=[e_ij] carries no passivity constraint; - * a negative eigenvalue is a "negative capacitance" that makes the transient DAE - * unstable (it diverges) even though every pole is in the LHP and AC is exact. - * Symmetrizing and clamping negative eigenvalues to 0 yields a passive C-matrix - * -> stable transient, while genuinely-improper (shunt-C) networks keep their - * positive eigenvalues. e[] is indexed i*N+j. */ -static void psd_project_E(double *e, int N) -{ - double *A = (double*) malloc((size_t) N*N*sizeof(double)); - double *V = (double*) malloc((size_t) N*N*sizeof(double)); - double *w = (double*) malloc((size_t) N*sizeof(double)); - int i, j, k; - for (i = 0; i < N; i++) for (j = 0; j < N; j++) A[i*N+j] = 0.5*(e[i*N+j]+e[j*N+i]); - jacobi_sym(A, N, w, V); - for (i = 0; i < N; i++) for (j = 0; j < N; j++) { - double acc = 0.0; - for (k = 0; k < N; k++) { double wk = w[k] > 0.0 ? w[k] : 0.0; acc += V[i*N+k]*wk*V[j*N+k]; } - e[i*N+j] = acc; - } - free(A); free(V); free(w); -} - -/* One-sided Jacobi SVD of a square real n x n matrix: A = U * diag(S) * V^T, - * S descending. U, S, V preallocated (n*n, n, n*n). Robust, no external deps; - * used to detect and factor low-rank coupling in the emit (Enhancement-205). */ -static void jacobi_svd(const double *Ain, int n, double *U, double *S, double *V) -{ - int i, p, q, sweep; - for (i = 0; i < n*n; i++) U[i] = Ain[i]; - for (i = 0; i < n; i++) { int t; for (t = 0; t < n; t++) V[i*n+t] = (i==t)?1.0:0.0; } - for (sweep = 0; sweep < 60; sweep++) { - int rotated = 0; - for (p = 0; p < n-1; p++) for (q = p+1; q < n; q++) { - double alpha=0, beta=0, gamma=0; - for (i = 0; i < n; i++) { double ap=U[i*n+p], aq=U[i*n+q]; - alpha += ap*ap; beta += aq*aq; gamma += ap*aq; } - if (fabs(gamma) <= 1e-15*sqrt(alpha*beta)) continue; - rotated = 1; - double zeta = (beta - alpha) / (2.0*gamma); - double t = (zeta>=0?1.0:-1.0) / (fabs(zeta) + sqrt(1.0+zeta*zeta)); - double c = 1.0/sqrt(1.0+t*t), sn = c*t; - for (i = 0; i < n; i++) { - double ap=U[i*n+p], aq=U[i*n+q]; - U[i*n+p] = c*ap - sn*aq; U[i*n+q] = sn*ap + c*aq; - double vp=V[i*n+p], vq=V[i*n+q]; - V[i*n+p] = c*vp - sn*vq; V[i*n+q] = sn*vp + c*vq; - } - } - if (!rotated) break; - } - for (q = 0; q < n; q++) { - double nrm = 0.0; for (i = 0; i < n; i++) nrm += U[i*n+q]*U[i*n+q]; - nrm = sqrt(nrm); S[q] = nrm; - if (nrm > 1e-300) for (i = 0; i < n; i++) U[i*n+q] /= nrm; - } - for (p = 0; p < n-1; p++) { /* sort columns by S descending (n small) */ - int mx = p; for (q = p+1; q < n; q++) if (S[q] > S[mx]) mx = q; - if (mx != p) { double t=S[p]; S[p]=S[mx]; S[mx]=t; - for (i = 0; i < n; i++) { double a=U[i*n+p]; U[i*n+p]=U[i*n+mx]; U[i*n+mx]=a; - double b=V[i*n+p]; V[i*n+p]=V[i*n+mx]; V[i*n+mx]=b; } } - } -} - -/* Emit the "num, den" coefficient arrays of one pole section's laplace_nd: - * kind 0 = real pole (1/(s-p)); 1 = conj-pair "cos" (s-sigma)/D; 2 = "sin" om/D. */ -static void emit_filter(FILE *fo, cplx pole, int kind) -{ - if (kind == 0) { - double num[1] = { 1.0 }, de[2] = { -creal(pole), 1.0 }; - emit_arr(fo, num, 1); fprintf(fo, ", "); emit_arr(fo, de, 2); - } else { - double sig = creal(pole), om = cimag(pole); - double de[3] = { sig*sig+om*om, -2.0*sig, 1.0 }; - if (kind == 1) { double num[2] = { -sig, 1.0 }; emit_arr(fo, num, 2); } - else { double num[1] = { om }; emit_arr(fo, num, 1); } - fprintf(fo, ", "); emit_arr(fo, de, 3); - } -} - -/* ============================ public API ============================ */ -/* Shared front half: parse Touchstone, S->Y, common-pole vector fit with order - * selection, reciprocal mirror, PSD-project E. On success returns 0 and hands the - * caller freshly-owned fit arrays P[Np], res[N*N*Np], d[N*N], e[N*N] (caller frees - * them); the parse scaffolding is freed here. Layout is exactly what both emitters - * and the native `.nport` device expect: poles canonicalized (real first, then - * adjacent conjugate pairs); res indexed (i*N+j)*Np+k. Leak-tolerant for discarded - * candidate fits, as before -- a one-shot conversion tool. */ -static int snp_fit(const char *snpfile, int *pN, int *pNp, - cplx **pP, cplx **pRes, double **pD, double **pE, - double *pErr, char *msg, int msglen) -{ - TS ts; int i, j, k, r; - if (parse_touchstone(snpfile, &ts, msg, msglen)) return 1; - int N = ts.N, nf = ts.nf; - cplx *Y = (cplx*) malloc((size_t) nf*N*N*sizeof(cplx)); - if (to_Y(&ts, Y)) { snprintf(msg,(size_t)msglen,"S->Y conversion failed (singular)"); ts_free(&ts); free(Y); return 1; } - /* Y entries in row-major i*N+j order, each a length-nf vector */ - int Nf = N*N; - cplx *s = (cplx*) malloc((size_t) nf*sizeof(cplx)); - for (r = 0; r < nf; r++) s[r] = 2*M_PI*ts.freqs[r]*I; - double wn = sqrt(cabs(s[0])*cabs(s[nf-1])); - cplx *sn = (cplx*) malloc((size_t) nf*sizeof(cplx)); - for (r = 0; r < nf; r++) sn[r] = s[r]/wn; - cplx *F = (cplx*) malloc((size_t) Nf*nf*sizeof(cplx)); - for (i = 0; i < N; i++) for (j = 0; j < N; j++) for (r = 0; r < nf; r++) - F[(long)(i*N+j)*nf+r] = Y[(long)r*N*N + i*N+j]; - - /* reciprocity: a passive network gives symmetric Y, so fit only the upper - * triangle (N(N+1)/2 elements) and mirror -- ~2x fewer LS solves (Fix #2). */ - int reciprocal = is_reciprocal(F, N, nf); - int *elems = (int*) malloc((size_t) Nf*sizeof(int)); - int Ne = 0; - if (reciprocal) { for (i = 0; i < N; i++) for (j = i; j < N; j++) elems[Ne++] = i*N+j; } - else { for (i = 0; i < Nf; i++) elems[Ne++] = i; } - - /* ---- order selection: climb, keep best STABLE fit, knee near a floor ---- */ - double fmin = ts.freqs[0], fmax = ts.freqs[nf-1], tol = 1e-3; - cplx *bestP=NULL,*bestRes=NULL; double *bestD=NULL,*bestE=NULL; int bestNp=0; double bestErr=1e300; - cplx *prevP=NULL,*prevRes=NULL; double *prevD=NULL,*prevE=NULL; int prevNp=0; double prevErr=-1, firstErr=-1; - int chosenP=0; cplx *chP=NULL,*chRes=NULL; double *chD=NULL,*chE=NULL; - int npair; - for (npair = 1; npair <= 12; npair++) { - int Np = 2*npair; - cplx *P = (cplx*) malloc((size_t) Np*sizeof(cplx)); - cplx *res = (cplx*) malloc((size_t) Nf*Np*sizeof(cplx)); - double *dd = (double*) malloc((size_t) Nf*sizeof(double)); - double *ee = (double*) malloc((size_t) Nf*sizeof(double)); - cplx *p0 = (cplx*) malloc((size_t) Np*sizeof(cplx)); - seed_poles(fmin, fmax, npair, p0); - for (i = 0; i < Np; i++) P[i] = p0[i]/wn; - vector_fit(sn, nf, F, elems, Ne, Np, P, res, dd, ee, 10); - /* un-normalize (fitted elements only) */ - for (i = 0; i < Np; i++) P[i] *= wn; - for (int ei = 0; ei < Ne; ei++) { k = elems[ei]; for (i = 0; i < Np; i++) res[(long)k*Np+i] *= wn; ee[k] /= wn; } - /* rms rel error + stability (over fitted elements; Y symmetric so representative) */ - double err = 0.0; int stable = 1; - for (i = 0; i < Np; i++) if (creal(P[i]) > 1e-6) stable = 0; - for (int ei = 0; ei < Ne; ei++) { - k = elems[ei]; - double numr=0, denr=0; - for (r = 0; r < nf; r++) { - cplx fit = dd[k] + s[r]*ee[k]; - for (i = 0; i < Np; i++) fit += res[(long)k*Np+i]/(s[r]-P[i]); - cplx dif = fit - F[(long)k*nf+r]; - numr += creal(dif)*creal(dif)+cimag(dif)*cimag(dif); - cplx fv = F[(long)k*nf+r]; - denr += creal(fv)*creal(fv)+cimag(fv)*cimag(fv); - } - double e2 = sqrt(numr)/(sqrt(denr)+1e-300); - if (e2 > err) err = e2; - } - if (!(err==err)) stable = 0; /* NaN */ - free(p0); - if (firstErr < 0) firstErr = err; - int keep_best = stable && err < bestErr; - if (keep_best) { - /* Don't free the old best buffers if `prev` still aliases them (that - * happens after any prior keep_best iteration, where best==prev==P): - * the prev-shift below frees them exactly once. Freeing here would - * leave prev dangling and double-free at line ~490. */ - if (bestP != prevP) { free(bestP);free(bestRes);free(bestD);free(bestE); } - bestP=P;bestRes=res;bestD=dd;bestE=ee;bestNp=Np;bestErr=err; - } - if (!stable) { if(!keep_best){free(P);free(res);free(dd);free(ee);} break; } - if (err < tol) { chosenP=Np; chP=P;chRes=res;chD=dd;chE=ee; if(keep_best){/*owned by best too*/} break; } - int near_floor = (err < 0.1*firstErr) || (err < 0.05); - if (prevErr >= 0 && err > 0.7*prevErr && near_floor) { /* knee at floor -> use prev */ - chosenP=prevNp; chP=prevP;chRes=prevRes;chD=prevD;chE=prevE; - if (!keep_best) { free(P);free(res);free(dd);free(ee); } - break; - } - /* shift prev <- current (free old prev unless it is the best) */ - if (prevP && prevP!=bestP) { free(prevP);free(prevRes);free(prevD);free(prevE); } - prevP=P;prevRes=res;prevD=dd;prevE=ee;prevNp=Np;prevErr=err; - } - /* pick chosen, else best, else prev */ - cplx *P; cplx *res; double *dd,*ee; int Np; - if (chP) { P=chP;res=chRes;dd=chD;ee=chE;Np=chosenP; } - else if (bestP) { P=bestP;res=bestRes;dd=bestD;ee=bestE;Np=bestNp; } - else { P=prevP;res=prevRes;dd=prevD;ee=prevE;Np=prevNp; } - - /* mirror the fitted upper triangle into the lower one (reciprocal case, Fix #2) */ - if (reciprocal) { - for (i = 0; i < N; i++) for (j = i+1; j < N; j++) { - long u = i*N+j, l = j*N+i; - for (k = 0; k < Np; k++) res[l*Np+k] = res[u*Np+k]; - dd[l] = dd[u]; ee[l] = ee[u]; - } - } - - /* force the improper (e*s) capacitance matrix passive so transient is stable */ - psd_project_E(ee, N); - - *pN = N; *pNp = Np; *pP = P; *pRes = res; *pD = dd; *pE = ee; *pErr = bestErr; - free(Y); free(s); free(sn); free(F); free(elems); ts_free(&ts); - return 0; -} - -/* ---------- Touchstone -> Verilog-A (structured laplace_nd realization) ---------- */ -int snp2va_convert(const char *snpfile, const char *vafile, const char *module, - char *msg, int msglen) -{ - int N, Np, i, j, k; - cplx *P, *res; double *dd, *ee, bestErr; - if (snp_fit(snpfile, &N, &Np, &P, &res, &dd, &ee, &bestErr, msg, msglen)) return 1; - - /* ---- emit VA (shared-pole realization; Fix #4) ---- - * All N^2 elements share the SAME poles, so realize the pole-filters ONCE per - * input port and form each output current as a cheap weighted sum, instead of - * one independent laplace_nd bank per element. That is N*Np filter sections and - * O(N*Np) OSDI state, not O(N^2*Np) -- the difference between a model that - * compiles/simulates at large N and one that does not. Each section is still a - * well-conditioned 1st/2nd-order laplace_nd (a real pole -> res/(s-p); a conj - * pair -> a real "cos" basis (s-sigma)/D and "sin" basis omega/D, with the - * residue entering as a real weight), so coefficients stay <= O(|p|^2) and the - * transient is stable. d is a plain conductance, e*s a ddt (PSD-projected). */ - int nsec = 0; - int *sc_pole = (int*) malloc((size_t) Np*sizeof(int)); - int *sc_kind = (int*) malloc((size_t) Np*sizeof(int)); /* 0 real, 1 cos, 2 sin */ - for (k = 0; k < Np; ) { - int is_pair = (k+1 < Np) && (fabs(cimag(P[k])) > 1e-6*cabs(P[k])); - if (!is_pair) { sc_pole[nsec]=k; sc_kind[nsec]=0; nsec++; k += 1; } - else { sc_pole[nsec]=k; sc_kind[nsec]=1; nsec++; - sc_pole[nsec]=k; sc_kind[nsec]=2; nsec++; k += 2; } - } - /* Enhancement-205: structured (diagonal + low-rank) emit. Build each - * "channel" weight matrix W (N x N, real) -- the conductance d, the - * capacitance e, and every pole section's residue weights -- and pick, per - * channel, the cheaper of a DENSE emit (a term per significant entry; one - * laplace_nd per input port) or a LOW-RANK emit W = U*V^T (rank r). Because - * laplace is linear, the low-rank form filters the r COMBINED inputs - * u_m = sum_j V[j][m]*V(p_j) ONCE (r filters, not N) and distributes them to - * the outputs via U -- O(N*r) terms and O(r*Np) filters. A full-rank device - * keeps the dense form (identical model, modulo dropping <1e-8 fit noise); - * a device whose ports couple through a few shared modes (multi-port - * filters/cavities/packages) collapses to a compact, fast-compiling model. */ - int force_dense = (getenv("PRE_SNP_DENSE") != NULL); /* escape hatch / A-B test */ - const double tol_drop = 1e-9; /* dense: drop only sub-1e-9 fit noise */ - const double tol_rank = 1e-7; /* low-rank: keep modes above 1e-7 (well below any - * verify tolerance, so a genuinely low-rank device - * with a clean singular-value gap compresses fully - * and is reproduced essentially exactly, while a - * slowly-decaying channel just stays dense). */ - int nchan = 2 + nsec, c, m; - double **chW = (double**) malloc((size_t) nchan*sizeof(double*)); - double **chU = (double**) calloc((size_t) nchan, sizeof(double*)); - double **chV = (double**) calloc((size_t) nchan, sizeof(double*)); - double *chMx = (double*) malloc((size_t) nchan*sizeof(double)); - int *ch_lr = (int*) calloc((size_t) nchan, sizeof(int)); - int *ch_r = (int*) calloc((size_t) nchan, sizeof(int)); - int *ch_kind = (int*) malloc((size_t) nchan*sizeof(int)); /* 0=V, 1=ddt, 2=filtered */ - int *ch_sec = (int*) malloc((size_t) nchan*sizeof(int)); - double *svU=(double*)malloc((size_t)N*N*sizeof(double)); - double *svS=(double*)malloc((size_t)N*sizeof(double)); - double *svV=(double*)malloc((size_t)N*N*sizeof(double)); - int nfilt = 0; - for (c = 0; c < nchan; c++) { - double *W = (double*) malloc((size_t) N*N*sizeof(double)); - if (c == 0) { ch_kind[c]=0; ch_sec[c]=-1; for (i=0;imx) mx=fabs(W[i]); - chMx[c]=mx; - if (mx==0.0) continue; /* zero channel */ - int nnz=0; for (i=0;itol_drop*mx) nnz++; - jacobi_svd(W, N, svU, svS, svV); - int rnk=0; for (i=0;i tol_rank*svS[0]) rnk++; if (rnk<1) rnk=1; - int filt = (ch_kind[c]==2); - long lr_cost = 2L*N*rnk + (filt? rnk : 0); - long de_cost = (long) nnz + (filt? N : 0); - if (!force_dense && lr_cost < de_cost) { - ch_lr[c]=1; ch_r[c]=rnk; nfilt += filt? rnk : 0; - double *Uc=(double*)malloc((size_t)N*rnk*sizeof(double)); - double *Vc=(double*)malloc((size_t)N*rnk*sizeof(double)); - for (i=0;i1e-30){ term_sep(fo,&first); - if (ch_kind[c]==2) fprintf(fo,"(%.12g)*h%d_%d",u,c,m); else fprintf(fo,"(%.12g)*g%d_%d",u,c,m); } } - } else { - double thr = tol_drop*chMx[c]; - for (j=0;jthr){ term_sep(fo,&first); - if (ch_kind[c]==0) fprintf(fo,"(%.12g)*V(p%d,ref)",w,j+1); - else if (ch_kind[c]==1) fprintf(fo,"(%.12g)*ddt(V(p%d,ref))",w,j+1); - else fprintf(fo,"(%.12g)*w%d_%d",w,c,j); } } - } - } - if (first) fprintf(fo, "0.0"); fprintf(fo, ";\n"); - } - fprintf(fo, " end\nendmodule\n"); - fclose(fo); - for (c=0;c native `.nport` fit file (Enhancement-242) ---------- - * Same shared parse+fit, written directly as the compact pole/residue description - * the built-in n-port device reads -- no Verilog-A, no openvaf-r compile. The - * device evaluates the identical model Y_ij(s)=d+s*e+sum_k res/(s-p_k), so a - * `-native` run reproduces a `-osdi` run to fit accuracy. */ -int snp2nport_convert(const char *snpfile, const char *nportfile, - char *msg, int msglen) -{ - int N, Np, i, j, k; - cplx *P, *res; double *dd, *ee, bestErr; - if (snp_fit(snpfile, &N, &Np, &P, &res, &dd, &ee, &bestErr, msg, msglen)) return 1; - - FILE *fo = fopen(nportfile, "w"); - if (!fo) { snprintf(msg,(size_t)msglen,"cannot write '%s'", nportfile); return 1; } - fprintf(fo, "NPORT 1\n"); - fprintf(fo, "# generated by pre_snp -native from %s\n", snpfile); - fprintf(fo, "nports %d\nnpoles %d\n", N, Np); - fprintf(fo, "poles\n"); - for (k = 0; k < Np; k++) - fprintf(fo, " %.15e %.15e\n", creal(P[k]), cimag(P[k])); - fprintf(fo, "d\n"); - for (i = 0; i < N; i++) { - for (j = 0; j < N; j++) fprintf(fo, " %.15e", dd[i*N+j]); - fprintf(fo, "\n"); - } - fprintf(fo, "e\n"); - for (i = 0; i < N; i++) { - for (j = 0; j < N; j++) fprintf(fo, " %.15e", ee[i*N+j]); - fprintf(fo, "\n"); - } - fprintf(fo, "res\n"); - for (i = 0; i < N; i++) for (j = 0; j < N; j++) - for (k = 0; k < Np; k++) { - long idx = (long)(i*N+j)*Np + k; - fprintf(fo, " %.15e %.15e\n", creal(res[idx]), cimag(res[idx])); - } - fclose(fo); - snprintf(msg,(size_t)msglen,"%d-port, %d poles, rms rel err %.2e", N, Np, bestErr<1e300?bestErr:0.0); - return 0; -} - -#ifdef SNP2VA_TEST -int main(int argc, char **argv) -{ - setbuf(stderr, NULL); - if (argc < 3) { fprintf(stderr, "usage: %s file.sNp out.va [module]\n", argv[0]); return 2; } - char msg[256]; - int rc = snp2va_convert(argv[1], argv[2], argc>3?argv[3]:"nport", msg, sizeof msg); - fprintf(stderr, "snp2va: %s\n", msg); - return rc; -} -#endif