// determinant.c
#include <stdlib.h>

#define abs(x) ((x)<0.0?(-(x)):(x))
double determinant(int n, double A[n][n])
{
  double D = 1.0;  /* DETERMINANT */
  double *B;  /* WORKING MATRIX */
  int *ROW;  /* ROW INTERCHANGE INDICIES */
  int HOLD , I_PIVOT;  /* PIVOT INDICIES */
  double PIVOT;  /* PIVOT ELEMENT VALUE */
  double ABS_PIVOT;
  int i, j, k;

  B = (double *)calloc(n*n, sizeof(double));
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      B[i*n+j] = A[i][j];
    }
  }
  ROW = (int *)calloc(n, sizeof(int));
  for(k=0; k<n; k++)
  {
    ROW[k]= k;
  }
  for(k=0; k<n-1; k++)
  {
    PIVOT = B[ROW[k]*n+k];
    ABS_PIVOT = abs(PIVOT);
    I_PIVOT = k;
    for(i=k; i<n; i++)
    {
      if(abs(B[ROW[i]*n+k]) > ABS_PIVOT )
      {
        I_PIVOT = i;
        PIVOT = B[ROW[i]*n+k];
        ABS_PIVOT = abs(PIVOT);
      }
    }
    if(I_PIVOT != k)
    {
      HOLD = ROW[k];
      ROW[k] = ROW[I_PIVOT];
      ROW[I_PIVOT] = HOLD;
      D = - D;
    }
    if(ABS_PIVOT < 1.0E-20)
    {
      free(B);
      free(ROW);
      return ABS_PIVOT;
    }
    else
    {
      D = D * PIVOT;
      for(j=k+1; j<n; j++)
      {
        B[ROW[k]*n+j] = B[ROW[k]*n+j] / B[ROW[k]*n+k];
      }
      for(i=0; i<n; i++)
      {
        if(i != k)
        {
          for(j=k+1; j<n; j++)
          {
            B[ROW[i]*n+j] = B[ROW[i]*n+j] - B[ROW[i]*n+k]* B[ROW[k]*n+j];
          } /* end j */
        } /* end if */
      } /* end i */
    } /* end else */
  } /* end k */
  D = D * B[ROW[n-1]*n+n-1];
  free(B);
  free(ROW);
  return D;
} /* end determinant */