Задача методом Гаусса

//-----------------------------------------------------
//                        gauss
//  Решение системы линейных уравнений методом Гаусса
//
// Обращение: int gauss (int N, double far *A,
//                       double far *B, double far *X)
// Входные параметры:
// int N     - количество уравнений;
// double far *A - указатель на матрицу коэффициентов
//             размерностью NxN;
// double far *B - указатель на столбец свободных членов
//             размерностью N.
//
// Выходные параметры:
// double far *X - указатель на столбец неизвестных
//             размерностью N.
// Процедура возвращает: 0 - нормальное завершение;
//                      -1 - система несовместна.
//
//----------------------------------------------------------------
 
int gauss(int N, double far *A, double far *B, double far *X)
    {
     int i, j, k;
     double c;
 
 for(k = 0; k < N-1; k++)
    {j = 0; while(A[k*N+k] == 0.0)
        {if(j++ > N-k-1) return(-1);
         for(i = 0; i <   N; i++) X[i] = A[k*N+i]; c = B[k];
         for(i = k; i < N-1; i++)
            {for(j = 0; j < N; j++) A[i*N+j] = A[(i+1)*N+j];
                                      B[i] = B[i+1];}
         for(i = 0; i < N; i++) A[(N-1)*N+i] = X[i]; B[N-1] = c;
        }
 
 for(i = k+1; i < N; i++)
    {c = A[i*N+k]/A[k*N+k]; A[i*N+k] = 0.0;
     for(j = k+1; j < N; j++) A[i*N+j] -= c*A[k*N+j];
     B[i] -= c*B[k];
     }
 
 if(A[(N-1)*N+N-1] == 0.0) return(-1);
 X[N-1] = B[N-1]/A[(N-1)*N+N-1];
 for(i = N-2; i >= 0; i--)
    {c = 0.0;
     for(j = i+1; j < N; j++) c += A[i*N+j]*X[j];
     X[i] = (B[i]-c)/A[i*N+i];
    } }
  return(0);
}