Решение СЛАУ

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (square root method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    From given matrix A we build two auxulary triangular matrixes S & D
//    Then solving equations Gy=b (where g[I][j]=s[j][I]*d[j]) and Sx=y
//    we obtain vector x.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N    4     // size of matrix
#define N1   N+1
float matrix[N][N1]=
    {{10.9,1.2,2.1,0.9,23.2},
  {1.2,11.2,1.5,2.5,38.1},
  {2.7,1.5,9.8,1.3,40.3},
  {0.9,2.5,1.3,12.1,58.2}
    };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  register char i,j,k;
  float s[N][N],x[N],y[N],d[N],sum;
  // Printing given matrix
  ShowMatrix();
  // Building matrixes S and D
  for (j=0;j<N;j++)
    for (i=0;i<=j;i++)
    {
      sum=matrix[i][j];
      for (k=0;k<i;k++)
        sum-=s[k][i]*d[k]*s[k][j];
      if (i==j)
      {
        d[i]=(sum>0?1:-1);
        s[i][j]=sqrt(fabs(sum));
      }
      else s[i][j]=sum/(d[i]*s[i][i]);
    }
  // Solving equation Gy=b (G: g[I][j]=s[j][I]*d[j])
  for (i=0;i<N;i++)
  {
    y[i]=matrix[i][N];
    for (j=0;j<i;j++) y[i]-=s[j][i]*d[j]*y[j];
    y[i]/=s[i][i]*d[i];
  }
  // Solving equation Sx=y
  for (i=N-1;i>=0;i--)
  {
    x[i]=y[i];
    for (j=i+1;j<N;j++) x[i]-=s[i][j]*x[j];
    x[i]/=s[i][i];
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
  printf("x%d=%f\n",i+1,x[i]);
}