c++ double** jacobi SVD

以下是C++中使用Jacobi SVD算法实现矩阵奇异值分解的示例代码：

#include <iostream>
#include <cmath>

void jacobiSVD(double** A, int m, int n, double** U, double* D, double** Vt) {
    // 初始化临时变量D、U、Vt
    for (int i = 0; i < n; i++) {
        D[i] = 0;
        for (int j = 0; j < m; j++) {
            U[j][i] = 0;
        }
        for (int j = 0; j < n; j++) {
            Vt[i][j] = 0;
        }
    }

    // 迭代计算
    double* B = new double[n];
    double* Z = new double[n];
    double* V = new double[n];
    double f, h, c, s, x, y, z;

    for (int i = 0; i < n; i++) {
        B[i] = D[i] = A[i][i];
        U[i][i] = 1;
        V[i] = 0;
    }

    for (int k = 0; k < 50; k++) {
        double sm = 0;
        for (int i = 0; i < n - 1; i++) {
            for (int j = i + 1; j < n; j++) {
                sm += std::abs(A[i][j]);
            }
        }

        if (sm == 0) {
            break;
        }

        double tresh = (k < 3) ? 0.2 * sm / (n * n) : 0;

        for (int i = 0; i < n - 1; i++) {
            for (int j = i + 1; j < n; j++) {
                double g = 100 * std::abs(A[i][j]);
                if (k > 3 && std::abs(D[i]) + g == std::abs(D[i]) && std::abs(D[j]) + g == std::abs(D[j])) {
                    A[i][j] = 0;
                }
                else if (std::abs(A[i][j]) > tresh) {
                    h = D[j] - D[i];
                    if (std::abs(h) + g == std::abs(h)) {
                        t = A[i][j] / h;
                    }
                    else {
                        theta = 0.5 * h / A[i][j];
                        t = 1 / (std::abs(theta) + std::sqrt(1 + theta * theta));
                        if (theta < 0) {
                            t = -t;
                        }
                    }
                    c = 1 / std::sqrt(1 + t * t);
                    s = t * c;
                    tau = s / (1 + c);
                    h = t * A[i][j];
                    Z[i] -= h;
                    Z[j] += h;
                    D[i] -= h;
                    D[j] += h;
                    A[i][j] = 0;

                    for (int k = 0; k < i; k++) {
                        g = A[k][i];
                        h = A[k][j];
                        A[k][i] = g - s * (h + g * tau);
                        A[k][j] = h + s * (g - h * tau);
                    }

                    for (int k = i + 1; k < j; k++) {
                        g = A[i][k];
                        h = A[k][j];
                        A[i][k] = g - s * (h + g * tau);
                        A[k][j] = h + s * (g - h * tau);
                    }

                    for (int k = j + 1; k < n; k++) {
                        g = A[i][k];
                        h = A[j][k];
                        A[i][k] = g - s * (h + g * tau);
                        A[j][k] = h + s * (g - h * tau);
                    }

                    for (int k = 0; k < n; k++) {
                        g = U[k][i];
                        h = U[k][j];
                        U[k][i] = g - s * (h + g * tau);
                        U[k][j] = h + s * (g - h * tau);
                    }
                }
            }
        }

        D[i] += Z[i];
        B[i] = D[i];
        Z[i] = 0;
    }

    delete[] B;
    delete[] Z;
    delete[] V;
}

int main() {
    // 示例用法
    int m = 3;
    int n = 3;
    double** A = new double*[n];
    double** U = new double*[m];
    double* D = new double[n];
    double** Vt = new double*[n];
    for (int i = 0; i < n; i++) {
        A[i] = new double[n];
        U[i] = new double[m];
        Vt[i] = new double[n];
    }

    // 初始化A矩阵
    A[0][0] = 1;
    A[0][1] = 2;
    A[0][2] = 3;
    A[1][0] = 4;
    A[1][1] = 5;
    A[1][2] = 6;
    A[2][0] = 7;
    A[2][1] = 8;
    A[2][2] = 9;

    jacobiSVD(A, m, n, U, D, Vt);

    // 输出结果
    std::cout << "U:" << std::endl;
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
            std::cout << U[i][j] << " ";
        }
        std::cout << std::endl;
    }

    std::cout << "D:" << std::endl;
    for (int i = 0; i < n; i++) {
        std::cout << D[i] << " ";
    }
    std::cout << std::endl;

    std::cout << "Vt:" << std::endl;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            std::cout << Vt[i][j] << " ";
        }
        std::cout << std::endl;
    }

    // 释放内存
    for (int i = 0; i < n; i++) {
        delete[] A[i];
        delete[] U[i];
        delete[] Vt[i];
    }
    delete[] A;
    delete[] U;
    delete[] D;
    delete[] Vt;

    return 0;
}