Коэффициенты системы линейных уравнений

#include <iostream>
#include <cmath>
 
void quit()
{
    std::cout << "Матрица системы - вырожденная" << std::endl;
    exit(EXIT_SUCCESS);
}
 
int main()
{
    setlocale(LC_ALL, "Russian");
 
    const size_t size=4;
    double A1[size][size];
    double B1[size];
    //Матрица системы
    double A[size][size]={ 
        {6.1, 6.2, -6.3, 6.4},
        {1.1, -1.5, 2.2, -3.8},
        {5.1, -5, 4.9, -4.8},
        {1.8, 1.9, 2, -2.1}};
 
    //Вектор, определяющий правую часть с-мы
    double B[size]={6.5, 1.2, 4.7, 2.2};
    //Число перестановок строк
    size_t rCnt=0;
    //Ведущий элемент
    double ldn;
    //Номер ведущего элемента
    size_t index;
    //Неизвестные корни
    double x[size];
    //Вектор невязки
    double r[size]={0, 0, 0, 0};
    //Вспомогательный вектор
    double asv[size]={0, 0, 0, 0};
    //Определитель
    double det;
 
    //Копирование исходных массивов
    for(size_t i=0; i<size; ++i)
    {
        for(size_t j=0; j<size; ++j)
            A1[i][j]=A[i][j];
        B1[i]=B[i];
    }
 
    //Прямой ход
    for(size_t j=0; j<size-1; ++j)
    {
        //Выбор ведущего элемента столбца
        ldn=fabs(A[j][j]);
        index=j;
        for(size_t i=j; i<size; ++i)
            if(fabs(A[i][j])>ldn)
            {
                ldn=fabs(A[i][j]);
                index=i;
            }
        //Если ведущий элемент равен нулю
        if(!ldn)
            quit();//Выход
 
        //Если ведущий элемент - не на гл. диагонали
        //меняем строки местами
        if(index!=j)
        {
            ++rCnt;
            double temp=B[j];
            B[j]=B[index];
            B[index]=temp;
            for(size_t i=0; i<size; ++i)
            {
                temp=A[j][i];
                A[j][i]=A[index][i];
                A[index][i]=temp;
            }
        }
        //Обнуление элементов ниже гл. диагонали
        for(size_t i=j+1; i<size; ++i)
        {
            double coeff=-(A[i][j]/A[j][j]);
            for(size_t k=j+1; k<size; ++k)
                A[i][k]=A[i][k]+A[j][k]*coeff;
            B[i]=B[i]+B[j]*coeff;
        }
    }
 
    det=1;
    //Подсчет определителя матрицы
    for(size_t i=0; i<size; i++)
        det*=A[i][i];
    det*=(rCnt%2)? (-1) : 1;
 
    //Обратная подстановка
    for(int j=size-1; j>=0; --j)
    {
        x[j]=B[j]/A[j][j];
        if(j!=3)
            for(int i=size-1; i>j; --i)
                x[j]=x[j]-A[j][i]*x[i]/A[j][j];
    }
 
    //Вычисление вектора невязки
    for(size_t i=0; i<size; ++i)
    {
        //Вычисление вспомогательного вектора
        for(size_t j=0; j<size; ++j)
            asv[i]+=x[j]*A1[i][j];
        r[i]=B1[i]-asv[i];
    }
 
    //Вывод значений
    std::cout << "\tРешение системы линейных уравнений методом Гаусса"
                 "\n\tс частичным выбором ведущего элемента"
                 "\n\nРешение системы:" << std::endl;
    for(size_t i=0; i<size; ++i, std::cout << std::endl)
        std::cout << "x[" << i << "]=" << x[i];
    std::cout << "\nОпределитель равен " << det
              << "\n\nВектор невязки:" << std::endl;   
    for(size_t i=0; i<size; ++i, std::cout << std::endl)
        std::cout << "r[" << i << "]=" << r[i];
 
    system("pause");
    return EXIT_SUCCESS;
}