C++ 实现Cholesky分解

#include "iostream"
#include <vector>
using namespace std;

int chufa = 0;
int chengfa = 0;

//A是对称阵，U是上三角阵，D储存对角阵对角线上元素
void Cholesky(vector<vector<double> > A, vector<vector<double> >&U , vector<double>&D) { 
    int n = A.size();
    for (int i = 0; i < n; ++i) {
        D[i] = A[i][i];
        vector<double> a = A[i];
        a[i] = 1;
        for (int j = i + 1; j < n; j++) {
            a[j] /= D[i];
            chufa += 1;
        }
        U[i] = a;
        for (int j = 0; j < n; j++) {
            A[i][j] = 0;
            A[j][i] = 0;
        }
        for (int j = i + 1; j < n; j++) {
            for (int k = j; k < n; k++) {
                A[j][k] -= D[i] * a[j] * a[k];
                A[k][j] = A[j][k];
                chengfa += 2;
            }
        }
    }
}

//输入对角线为1的上三角阵，返回其逆阵
vector<vector<double> > invU(vector<vector<double> > U) {
    vector<vector<double> > V(U.size(), vector<double>(U.size(),0));
    for (int col = U.size() - 1; col >= 0; col--) {
        V[col][col] = 1;
        for (int row = col-1; row >= 0; row--) 
            for (int k = row + 1; k <= col; k++) {
                V[row][col] -= U[row][k] * V[k][col];
                chengfa += 1;
            }
    }
    return V;
}

vector<vector<double> > transposition(vector<vector<double> > A) {
    vector<vector<double> > B(A[0].size(), vector<double>(A.size()));
    for (int i = 0; i < A.size(); i++)
        for (int j = 0; j < A[0].size(); j++)
            B[j][i] = A[i][j];
    return B; //返回输入的转置
}

vector<double> operator*(const vector<vector<double> > A, const vector<double> x) {
    vector<double> y(x.size());
    for (int i = 0; i < A.size(); i++) {
        for (int j = 0; j < x.size(); j++) {
            y[i] += A[i][j] * x[j];
            chengfa += 1;
        }
            
    }
    return y;  //y是矩阵A和向量x的乘积
}

int main()
{
    vector<vector<double> > A = {
        {3.3, -0.6, 0.1, -0.5},
        {-0.6, 3.2, -3.3, 0.3},
        {0.1, -3.3, 10.1, -1.8},
        {-0.5, 0.3, -1.8, 1.4}
    };
    vector<double> x = { 1.3, 0.2, 0.8, 1.4 };
    vector<vector<double> > U(A.size(), vector<double>(A.size()));
    vector<double> D(A.size());

    Cholesky(A, U, D);

    vector<vector<double> >V = invU(U);

    vector<double> y = transposition(V) * x;
    for (int i = 0; i < y.size(); i++) {
        y[i] /= D[i];
        chufa += 1;
    }
    y = V * y;
    return 0;
}