Методом Гаусса решить систему n линейных алгебраических уравнений

#include <windows.h>//malloc, system("pause")
#include <stdio.h>  //i/o
#include <conio.h>  //getch
#include <math.h>
 
void ShowVector(int n, double * vec);
void PryamoiHod(int n, double **a, double *b);
void ObratniHod(int n, double **a, double *b, double *x);
 
int main()
{
    int i,j,n;
    double **a, *b, *x;
    do
    {
        printf("Enter NUM of equations: ");
        scanf("%d",&n);
        //Выделяем память под матрицу А и векторы В и Х
        a = (double **)malloc(n*sizeof(double));
        b = (double *)malloc(n*sizeof(double));
        x = (double *)malloc(n*sizeof(double));
        for(i = 0; i < n; i++)
        {
            a[i] = (double *)malloc(n*sizeof(double));
            //Ввод a
            for(j = 0; j < n; j++)
            {
                printf("a[%d][%d] = ",i + 1,j + 1);
                scanf("%lf",&a[i][j]);
            }
        }
        //Ввод b
        for(i = 0; i < n; i++)
        {
            printf("b[%d] = ",i + 1);
            scanf("%lf",&b[i]);
        }
        
        printf("\tSee input\r\n");
        printf("Matrix A:\r\n");
        for(i = 0; i < n; i++)
            ShowVector(n, a[i]);
        printf("Vector B:\r\n");
        ShowVector(n, b);
        
        printf("\tSolving on Gauss method\r\n");
        PryamoiHod(n, a, b);
        printf("Forvard Gauss course\r\n");//Прямой ход
        printf("Matrix A:\r\n");
        for(i = 0; i < n; i++)
            ShowVector(n, a[i]);
        printf("Vector B:\r\n");
        ShowVector(n, b);
 
        ObratniHod(n, a, b, x);
        printf("Back Gauss course\r\n");//Обратный ход
        printf("Matrix A:\r\n");
        for(i = 0; i < n; i++)
            ShowVector(n, a[i]);
        printf("Vector B:\r\n");
        ShowVector(n, b);
 
        printf("Results :\r\n");
        ShowVector(n, x);
 
        printf("Press Y for new input\r\n");
        //Чистим память
        free((void *)a);
        free((void *)b);
        free((void *)x);
    }
    while(toupper(getch()) == 'Y');
    return 0;
}
 
void ShowVector(int n, double * vec)
{
    for(int i = 0; i < n; i++)
        printf("%.3f ",vec[i]);
    printf("\r\n");
}
 
void PryamoiHod(int n, double **a, double *b)
{
        double v;
        for(int k = 0,i,j,im; k < n - 1; k++)
        {
                im = k;
                for(i = k + 1; i < n; i++)
                {
                        if(fabs(a[im][k]) < fabs(a[i][k]))
                        {
                                im = i;
                        }
                }
                if(im != k)
                {
                        for(j = 0; j < n; j++)
                        {
                                v                = a[im][j];
                                a[im][j] = a[k][j];
                                a[k][j]  = v;
                        }
                        v     = b[im];
                        b[im] = b[k];
                        b[k]  = v;
                }
                for(i = k + 1; i < n; i++)
                {
                        v               = 1.0*a[i][k]/a[k][k];
                        a[i][k] = 0;
                        b[i]    = b[i] - v*b[k];
                        if(v != 0)
                        for(j = k + 1; j < n; j++)
                        {
                                a[i][j] = a[i][j] - v*a[k][j];
                        }
                }
        }
}
 
void ObratniHod(int n, double **a, double *b, double *x)
{
        double s = 0;
        x[n - 1] = 1.0*b[n - 1]/a[n - 1][n - 1];
        for(int i = n - 2, j; 0 <= i; i--)
        {
                s = 0;
                for(j = i + 1; j < n; j++)
                {
                        s = s+a[i][j]*x[j];
                }
                x[i] = 1.0*(b[i] - s)/a[i][i];
        }
}