代码¶

// ── 核心：One-sided Jacobi SVD (Hankel 矩阵用) ──
static int jacobi_svd(float *A, int m, int n, float *S, float *Vt, int max_iter)
{
    // Vt = I
    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++) {
                // 计算列 p, 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;

                // 计算 Givens 旋转角
                float zeta = (beta - alpha) / (2.0f * gamma);
                float t = (zeta >= 0 ? 1 : -1)
                          / (fabsf(zeta) + sqrtf(1 + zeta*zeta));
                float c = 1.0f / sqrtf(1 + t*t), s = c * t;

                // 左右旋转 A 和 Vt
                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;
    }
    // S[k] = 列范数, 降序排列
}

// ── 平衡实现 ──
static void era_build_A(const float *H0_svd, int rows, int cols_h0,
                         const float *H1, const float *S, const float *Vt,
                         int order, float *A_out)
{
    // A = Σ⁻¹/² · Uᵀ · H₁ · V · Σ⁻¹/²   (n×n)
    // 特征值分解 → 离散极点 → 连续模态参数
}