Cholesky decomposition/Разложение Холецкого на C++

cholesky_decomposition.cpp

#include <iostream>
#include <cmath>
using namespace std;

void FillMatrix(double **matrix, const int size){
    for (int i=0; i<size; i++){
        for (int j=0; j<size; j++){
            cout << "a[" << i+1 << "][" << j+1 << "] = ";
            cin >> matrix[i][j];
        }
    }
}

void PrintMatrix(double **matrix, const int size){
    for (int i=0; i<size; i++){
        for (int j=0; j<size; j++){
            cout << matrix[i][j] << "\t";
        }
        cout << endl;
    }
}

void clearMemory(double** matrix, int size){
    for (int i = 0; i < size; i++) {
        delete[] matrix[i];
    }
    delete [] matrix;
}

void Transpose(double **matrix, double **transposed, const int size){
    for (int i=0; i<size; i++){
        for (int j=0; j<size; j++){
            transposed[i][j] = matrix[j][i];
        }
    }
}

int main(){
    int size = 3;

    double **A = new double *[size];
    for (int i = 0; i < size; i++) {
        A[i] = new double [size];
    }
    double **Transposed = new double *[size];
    for (int i = 0; i < size; i++) {
        Transposed[i] = new double [size];
    }
    double **L = new double *[size];
    for (int i = 0; i < size; i++) {
        L[i] = new double [size];
    }

    //FillMatrix(A, size);
    A[0][0] =4;    A[0][1] =12;    A[0][2] =-16;
    A[1][0] =12;    A[1][1] =37;    A[1][2] =-43;
    A[2][0] =-16;    A[2][1] =-43;    A[2][2] =98;

    Transpose(A, Transposed, size);
    cout << "A =" << endl;
    PrintMatrix(A, size);

    for (int i=0; i<size; i++){
        for (int j=0; j<size; j++){
            if (Transposed[i][j] != A[i][j]){
                cout << "Матрица не является симметричной." << endl;
                exit(0);
            }
        }
    }

    for (int i = 0; i < size; i++) {
        for (int j = 0; j <= i; j++) {
            double sum = 0;
            for (int k = 0; k < j; k++) {
                sum += L[i][k] * L[j][k];
            }
            if (i == j) {
                L[i][j] = sqrt(A[i][i] - sum);
                if ((L[i][i] <= 0) || isnan(L[i][i])) {
                    cout << "Матрица не положительно определенная." << endl;
                    exit(0);
                }
            }
            else {
                L[i][j] = (1.0 / L[j][j] * (A[i][j] - sum));
            }
        }
    }

    cout << endl << "L =" << endl;
    PrintMatrix(L, size);

    clearMemory(A,size);
    clearMemory(Transposed, size);
    clearMemory(L, size);
    return 0;
}