// crout1.c  solve simultaneous equation  A x = y  by Crout method
//              converted from crout1.java, A[] single subscript


void crout1(int n, double A[], 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*n+k]*A[k*n+m];
      } // k
      A[i*n+m] = A[i*n+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*n+k]*A[k*n+j];
      } // k
      A[m*n+j] = (A[m*n+j]-sum)/A[m*n+m];
    } // j
    sum = 0.0;
    for(k=0; k<m; k++)
    {
      sum = sum + A[m*n+k]*y[k];
    } // k
    y[m] = (y[m]-sum)/A[m*n+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*n+k]*x[k];
    } // k
    x[i] = y[i] - sum;
  } // i
} // end crout1.c