Метод декомпозиции

#include <math.h>
#include <clocale> 
#include <iomanip>
#include <iostream>
using namespace std;
 
int main() {
    double A[4][4] = { {1, 2, 3, 4}, { 1,2,1,0 }, { -1,-2,2,0 }, { 0,1,1,0 } };
    double d[4] = { 7,1,1,2 };
    double B[4][4], C[4][4];
    const int n = 4;
    
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            B[i][j] = 0;
            C[i][j] = 0;
 
            if (i == j)
                C[i][j] = 1;
        }
    }
 
    for (int i = 0; i < n; i++)
    {
        B[i][0] = A[i][0];
        C[0][i] = A[0][i] / B[0][0];
    }
    
    for (int i = 1; i < n; i++)
    {
        double sum = 0;
        for (int j = i; j < n; j++) {
            for (int k = 0; k < j; k++)
            {
                sum += B[i][k] * C[k][j];
            }
            C[i][j] = A[i][j] - sum;
        }
 
        for (int j = i; j < n; j++) {
 
            for (int k = 0; k < i; k++)
            {
                sum += B[i][k] * C[k][j];
            }
            B[j][i] = (A[j][i] - sum) / C[i][i];
        }
    }
    
    printf("Матрица:\n");
    for (int i = 0; i < 4; i++)
    {
        for (int j = 0; j < 4; j++) {
            cout << setw(5u) << B[i][j];
        }
        printf("\n");
 
    }
 
    printf("Матрица:\n");
    for (int i = 0; i < 4; i++)
    {
        for (int j = 0; j < 4; j++) {
            cout << setw(5u) << C[i][j];
        }
        printf("\n");
 
    }
 
    return 0;
 
}