Реализация алгоритма решения СЛАУ методом Гаусса - C (СИ)

long double *gauss(int size, long double **main_mat, long double *main_b)
{
    int i, j, k;
    long double **mat = NULL;
    long double *b = NULL;
    long double *x = NULL;
    int *order = NULL;
    //делаем копии переданных данных
    if(!(mat = malloc(sizeof(long double*) * size))
        || !(b = malloc(sizeof(long double) * size))
        || !(x = malloc(sizeof(long double) * size))
        || !(order = malloc(sizeof(int) * size)))
    {
        free(mat);
        free(b);
        free(x);
        free(order);
        puts("!");
        return NULL;
    }
    for(i = 0; i < size; i++)
    {
        if(!(mat[i] = malloc(sizeof(long double) * size)))
        {
            for(; i >=0; i--)
                free(mat[i]);
            puts("!");
            return NULL;
        }
        memcpy(mat[i], main_mat[i], sizeof(long double) * size);
        order[i] = i;
    }
    memcpy(b, main_b, sizeof(long double) * size);
    
    puts("\nМетод Гаусса");
    for(k = 0; k < size; k++)
    {
        putchar('\n');
        long double max = fabsl(mat[k][k]);
        int max_i = k, max_j = k;
        if(k != size - 1)
        {
            //находим максимальный каэффициент
            for(i = k; i < size; i++)
            {
                for(j = k; j < size; j++)
                {
                    if(fabsl(mat[i][order[j]]) > max)
                    {
                        max = fabsl(mat[i][order[j]]);
                        max_i = i;
                        max_j = j;
                    }
                }
            }
            //он нулевой?..
            if(max == 0)
            {
                puts("Система ворожденная...");
                for(i = 0; i < size; i++)
                {
                    free(mat[i]);
                }
                free(mat);
                free(b);
                free(x);
                free(order);
                return NULL;
            }
            printf("Оптимальное основание [%d][%d] = %Lf\n",
                   max_i, max_j, mat[max_i][order[max_j]]);
        }
        else
            puts("Итоговая матрица");
        //по не обходимости меняем местами строки
        if(max_i != k)
        {
            long double *ptr = mat[k];
            mat[k] = mat[max_i];
            mat[max_i] = ptr;
            long double tmp = b[k];
            b[k] = b[max_i];
            b[max_i] = tmp;
        }
        //и/или переставляем столбцы
        if(max_j != k)
        {
            int tmp = order[k];
            order[k] = order[max_j];
            order[max_j] = tmp;
        }
        //делим на основание
        for(j = k + 1; j < size; j++)
            mat[k][order[j]] /= mat[k][order[k]];
        b[k] /= mat[k][order[k]];
        mat[k][order[k]] = 1;
        //вычитаем из оставшихся строк получившуюя с соответствующим коэффициентом 
        for(i = k + 1; i < size; i++)
        {
            long double factor = mat[i][order[k]];
            for(j = k; j < size; j++)
                mat[i][order[j]] -= mat[k][order[j]] * factor;
            b[i] -= b[k] * factor;
        }
        for(i = 0; i < size; i++)
            printf(" x%-7d ", order[i]);
        puts(" b");
        for(i = 0; i < size; i++)
        {
            for(j = 0; j < size; j++)
                printf("% 8Lf ", mat[i][order[j]]);
            printf("% 8Lf\n", b[i]);
        }
    }
    //совстевенно, считаем искомый вектор
    for(i = size - 1; i >= 0; i--)
    {
        for(j = i + 1; j < size; j++)
            b[i] -= x[order[j]] * mat[i][order[j]];
        x[order[i]] = b[i];
    }
    for(i = 0; i < size; i++)
    {
        free(mat[i]);
    }
    free(mat);
    free(b);
    free(order);
    return x;
}