// crout.c  solve simultaneous equation  A x = y  by Crout method
//              converted from crout.java, A[N][N] double subscript
//              N is dimensioned size of matrix, n is number of equations 


void crout(int n, int N, double A[N][N], double y[], double x[])
{
  double sum;
  int i, j, k, m;

  for(m=0; m<n; m++)
  {
    // below diagonal computations
    for(i=m; i<n; i++)
    {
      sum = 0.0;
      for(k=0; k<m; k++)
      {
        sum = sum + A[i][k]*A[k][m];
      } // k
      A[i][m] = A[i][m] - sum;
    } // i

    // above diagonal computations
    for(j=m+1; j<n; j++)
    {
      sum = 0.0;
      for(k=0; k<m; k++)
      {
        sum = sum + A[m][k]*A[k][j];
      } // k
      A[m][j] = (A[m][j]-sum)/A[m][m];
    } // j
    sum = 0.0;
    for(k=0; k<m; k++)
    {
      sum = sum + A[m][k]*y[k];
    } // k
    y[m] = (y[m]-sum)/A[m][m];
  } // m

  // back substitute
  x[0] = 0.0;
  for(i=n-1; i>=0; i--)
  {
    sum = 0.0;
    for(k=i+1; k<n; k++)
    {
      sum = sum + A[i][k]*x[k];
    } // k
    x[i] = y[i] - sum;
  } // i
} // end crout.c