Проблема с LU-разложением С++

1. Если взять алгоритм из Википедии (https://ru.wikipedia.org/wiki/LU-разложение), то реализовать его можно так, как показано ниже. Проверок там нет, - их стоит добавить.
2. Использование статических массивов, в данном случае, - не есть хорошо. Лучше использовать динамические, дабы программа могла оперировать с системами условно любой размерности.
3. Смешивать логику программы и интерфейс – тоже не есть хорошо. Логика – в одни функции / методы, ввод-вывод данных (включая сообщения) – в другие.
4. Если требуется граф. интерфейс (GUI), то прикрутить его к коду ниже (являющемуся кодом консольного приложения) не составит труда: в местах наполнения матриц и векторов случайными числами, а также в местах вывода результатов, нужно вставить соотв. функции / методы работы с элементами GUI.

Результаты работы программы:

A =
 42.000  68.000  35.000
  1.000  70.000  25.000
 79.000  59.000  63.000
B =
 65.000
  6.000
 46.000
L =
  1.000   0.000   0.000
  0.024   1.000   0.000
  1.881  -1.008   1.000
U =
 42.000  68.000  35.000
  0.000  68.381  24.167
  0.000   0.000  21.518
Y =
 65.000
  4.452
-71.775
X =
  2.313
  1.244
 -3.336

Проверка (https://ru.onlinemschool.com/math/assistance/equation/gaus/):

x1 = 142961/61801 = 2,313

x2 = 76876/61801 = 1,244

x3 = -206139/61801 = -3,336

Как видим, всё верно.

Собственно код:

#include "stdafx.h"
#include <iostream>
// печать матрицы
void printMatrix(float **matr, int n);
// печать вектор-столбца
void printVectorCol(float *v, int n);
//
float getSumU(int n, int i, int j, float **mL, float **mU);
//
float getSumL(int n, int i, int j, float **mL, float **mU);
// LU-разложение
void decomp(float **mA, int n, float **mL, float **mU);
//
float getSumY(int n, float **mL, float *vB, float *vY, int i);
//
float getSumX(int n, float **mU, float *vX, int i);
// решатель
void solver(float **mA, int n, float **mL, float **mU, float *vB, float *vY, float *vX);

int _tmain(int argc, _TCHAR* argv[])
{
    const int n = 3; // размерность системы
    // динамическое выделение памяти
    float **matrA = new float* [n];
    float **matrL = new float* [n];
    float **matrU = new float* [n];
    for (int i = 0; i < n; i++)
    {
        matrA[i] = new float [n];
        matrL[i] = new float [n];
        matrU[i] = new float [n];
    }
    float *vB = new float [n];
    float *vX = new float [n];
    float *vY = new float [n];
    // инициализация матриц и векторов нулями
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            matrA[i][j] = 0;
            matrL[i][j] = 0;
            matrU[i][j] = 0;
        }
        vB[i] = 0;
        vX[i] = 0;
        vY[i] = 0;
    }
    // заполнение матрицы A случайными числами из интервала [1; 100]
    for (int i = 0; i < n; i++) // по строкам
    {
        for (int j = 0; j < n; j++) // по столбцам
        {
            matrA[i][j] = rand() % 100 + 1;
        }
    }
    // задание случайных значений свободных коэффициентов
    for (int i = 0; i < n; i++)
    {
        vB[i] = rand() % 100 + 1;
    }
/***************************************************************************/
    // LU-разложение
    decomp(matrA, n, matrL, matrU);
    // решение системы
    solver(matrA, n, matrL, matrU, vB, vY, vX);
/***************************************************************************/
    // печать результатов
    printf("A = \n");
    printMatrix(matrA, n);
    printf("\n");
    printf("B = \n");
    printVectorCol(vB, n);
    printf("\n");
    printf("L = \n");
    printMatrix(matrL, n);
    printf("\n");
    printf("U = \n");
    printMatrix(matrU, n);
    printf("\n");
    printf("Y = \n");
    printVectorCol(vY, n);
    printf("\n");
    printf("X = \n");
    printVectorCol(vX, n);
    printf("\n");
/***************************************************************************/
    // освобождение памяти
    for (int i = 0; i < n; i++)
    {
        delete [] matrA[i];
        delete [] matrL[i];
        delete [] matrU[i];
    }
    delete [] matrA;
    delete [] matrL;
    delete [] matrU;
    delete [] vB;
    delete [] vX;
    delete [] vY;
    printf("\n");
    system("PAUSE");
    return 0;
}
// печать матрицы
void printMatrix(float **matr, int n)
{
    for (int i = 0; i < n; i++) // по строкам
    {
        for (int j = 0; j < n; j++) // по столбцам
        {
            printf("%7.3f ", matr[i][j]);
        }
        printf("\n"); // новая строка
    }
}
// печать вектор-столбца
void printVectorCol(float *v, int n)
{
    for (int i = 0; i < n; i++)
    {
        printf("%7.3f ", v[i]);
        printf("\n");
    }
}
//
float getSumU(int n, int i, int j, float **mL, float **mU)
{
    float res = 0;
    for (int k = 0; k < i; k++)
    {
        res += mL[i][k] * mU[k][j];
    }
    return res;
}
//
float getSumL(int n, int i, int j, float **mL, float **mU)
{
    float res = 0;
    for (int k = 0; k < i; k++)
    {
        res += mL[j][k] * mU[k][i];
    }
    return res;
}
// LU-разложение
void decomp(float **mA, int n, float **mL, float **mU)
{
    // i = 0
    for (int j = 0; j < n; j++)
    {
        mU[0][j] = mA[0][j];
        mL[j][0] = mA[j][0] / mU[0][0];
    }
    // i = 1 .. n-1
    for (int i = 1; i < n; i++)
    {
        for (int j = i; j < n; j++)
        {
            mU[i][j] = mA[i][j] - getSumU(n, i, j, mL, mU);
            mL[j][i] = (mA[j][i] - getSumL(n, i, j, mL, mU)) / mU[i][i];
        }
    }
}
//
float getSumY(int n, float **mL, float *vB, float *vY, int i)
{
    float res = 0;
    for (int k = 0; k < i; k++)
    {
        res += mL[i][k] * vY[k];
    }
    return res;
}
//
float getSumX(int n, float **mU, float *vX, int i)
{
    float res = 0;
    for (int k = n - 1; k > i; k--)
    {
        res += mU[i][k] * vX[k];
    }
    return res;
}
// решатель
void solver(float **mA, int n, float **mL, float **mU, float *vB, float *vY, float *vX)
{
    // Y. Решение прямой подстановкой
    for (int i = 0; i < n; i++)
    {
        vY[i] = (vB[i] - getSumY(n, mL, vB, vY, i)) / mL[i][i];
    }
    // X. Решение обратной подстановкой
    for (int i = n - 1; i >= 0; i--)
    {
        vX[i] = (vY[i] - getSumX(n, mU, vX, i)) / mU[i][i];
    }
}