Skip to content

CODE

ERA Implementation — NExT → Hankel SVD → Balanced Realization → Modal Parameters

static int jacobi_svd(float *A, int m, int n, float *S, float *Vt, int max_iter)
{
    for (int iter = 0; iter < max_iter; iter++) {
        int changed = 0;
        for (int p = 0; p < n - 1; p++)
            for (int q = p + 1; q < n; q++) {
                float alpha=0, beta=0, gamma=0;
                for (int r = 0; r < m; r++) {
                    float ap=A[r*n+p], aq=A[r*n+q];
                    alpha += ap*ap; beta += aq*aq; gamma += ap*aq;
                }
                if (fabsf(gamma) < 1e-6f * sqrtf(alpha * beta)) continue;
                float zeta = (beta-alpha)/(2*gamma);
                float t = (zeta>=0?1:-1)/(fabsf(zeta)+sqrtf(1+zeta*zeta));
                float c = 1/sqrtf(1+t*t), s = c*t;
                for (int r = 0; r < m; r++) {
                    float ap=A[r*n+p], aq=A[r*n+q];
                    A[r*n+p]=c*ap-s*aq; A[r*n+q]=s*ap+c*aq;
                }
                for (int r = 0; r < n; r++) {
                    float vp=Vt[p*n+r], vq=Vt[q*n+r];
                    Vt[p*n+r]=c*vp-s*vq; Vt[q*n+r]=s*vp+c*vq;
                }
                changed = 1;
            }
        if (!changed) break;
    }
}