/* simeq_lup.c                                        */
/* solve for X,  A X = Y, using L U = A decomposition */

#include <stdio.h>
#include <stdlib.h>
#undef  abs
#define abs(x) (((x)<0.0)?(-(x)):(x))

void simeq_lup(int n, double A[], double Y[], double X[])
{
  double *B;           /* [n][n+1]  WORKING MATRIX */
  int *ROW;            /* ROW INTERCHANGE INDICIES */
  double *by;          /* interchanged Y */
  int HOLD , I_PIVOT;  /* PIVOT INDICIES */
  double pivot;        /* PIVOT ELEMENT VALUE */
  int i, j, k, m, im, it;

  B = (double *)calloc(n*n, sizeof(double));
  ROW = (int *)calloc(n, sizeof(int));
  by = (double *)calloc(n, sizeof(double));

  for(m=0; m<n; m++) ROW[m] = m;
  for(m=0; m<n-1; m++)
  {
    pivot = abs(A[ROW[m]*n+m]);
    im = m;
    for(i=m+1; i<n; i++)
    {
      if(abs(A[ROW[i]*n+m])>pivot)
      {
        im = i;
        pivot = abs(A[ROW[i]*n+m]);
      }
    }
    it = ROW[im];
    ROW[im] = ROW[m];
    ROW[m] = it;
  }
  
  /* copy A into B with largest diagonal */
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      B[i*n+j] = A[ROW[i]*n+j];
    }
    by[i] = Y[ROW[i]];
  }

  for(m=0; m<n; m++)
  {
    /* compute column of L in B */
    for(i=m; i<n; i++)
    {
      B[i*n+m] = B[i*n+m];
      for(k=0; k<m; k++)
      {
        B[i*n+m] = B[i*n+m] - B[i*n+k]*B[k*n+m];
      }
    }
    /* compute row of U in B */
    for(j=m+1; j<n; j++)
    {
      B[m*n+j] = B[m*n+j];
      for(k=0; k<m; k++)
      {
        B[m*n+j] = B[m*n+j] - B[m*n+k]*B[k*n+j];
      }
      B[m*n+j] = B[m*n+j]/B[m*n+m];
    }
  }

  /* solve for X given by */
  for(i=0; i<n; i++)
  {
    for(k=0; k<i; k++)
    {
      by[i] = by[i] - B[i*n+k]*by[k];
    }
    by[i] = by[i]/B[i*n+i];
  }

  for(j=n-1; j>=0; j--)
  {
    X[j] = by[j];
    for(k=j+1; k<n; k++)
    {
      X[j] = X[j] - B[j*n+k]*X[k];
    }
  }
  free(B);
  free(ROW);
  free(by);
} /* end simeq_lup.c */