solve with LU decomposition

#include <stdio.h>
#define N (4)
int main(void) {
  double a[N][N] = {{4., -6., -10., -2.},
                    {-6., 8., 21., -1.},
                    {-10., 21., -20., 23.},
                    {-2., -1., 23., -18.}};
  double b[N] = {-46., 69., 64., -7.};
  double L[N][N] = {
      {0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}};
  double U[N][N] = {
      {1., 0., 0., 0.}, {0., 1., 0., 0.}, {0., 0., 1., 0.}, {0., 0., 0., 1.}};
  for (int k = 0; k < N; ++k) {
    for (int i = k; i < N; ++i) {
      double s = 0.;
      for (int j = 0; j < k; ++j) {
        s += L[i][j] * U[j][k];
      }
      L[i][k] = a[i][k] - s;
    }
    for (int j = k + 1; j < N; ++j) {
      double s = 0.;
      for (int i = 0; i < k; ++i) {
        s += L[k][i] * U[i][j];
      }
      U[k][j] = (a[k][j] - s) / L[k][k];
    }
  }
  double x[N], y[N];
  for (int i = 0; i < N; ++i) {
    double s = 0.;
    for (int j = 0; j < i; ++j) {
      s += L[i][j] * y[j];
    }
    y[i] = (b[i] - s) / L[i][i];
  }
  for (int i = N - 1; 0 <= i; --i) {
    double s = 0.;
    for (int j = i + 1; j < N; ++j) {
      s += U[i][j] * x[j];
    }
    x[i] = y[i] - s;
  }
  printf("x = {");
  for (int i = 0; i < N; ++i) {
    printf("%s%f", (i == 0 ? "" : ", "), x[i]);
  }
  puts("}");
}