/* invert.c  basic Gausian elimination version */
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#define abs(a) ((a)<0?(-(a)):(a))

void invert(int n, double AA[])
{
  int *ROW, *COL;
  double *TEMP;
  int HOLD , I_pivot , J_pivot;
  double pivot, abs_pivot;
  int i, j, k;

  ROW = (int *)calloc(n, sizeof(int));
  COL = (int *)calloc(n, sizeof(int));
  TEMP = (double *)calloc(n, sizeof(double));

  /*              set up row and column interchange vectors */
  for (k=0; k<n; k++)
  {
    ROW[k] = k ;
    COL[k] = k ;
  }

  /*               begin main reduction loop */
  for (k=0; k<n; k++)
  {
    /*             find largest element for pivot */
    pivot = AA[ROW[k]*n+COL[k]] ;
    I_pivot = k;
    J_pivot = k;
    for(i=k; i<n; i++)
    {
      for(j=k; j<n; j++)
      {
        abs_pivot = abs(pivot) ;
        if( abs(AA[ROW[i]*n+COL[j]]) > abs_pivot )
        {
          I_pivot = i ;
          J_pivot = j ;
          pivot = AA[ROW[i]*n+COL[j]] ;
        }
      }
    }
    if(abs(pivot) < 1.0E-65)
    {
      free(ROW);
      free(COL);
      free(TEMP);
      printf("MATRIX is SINGULAR !!! \n");
      return;
    }

    HOLD = ROW[k];
    ROW[k]= ROW[I_pivot];
    ROW[I_pivot] = HOLD ;
    HOLD = COL[k];
    COL[k]= COL[J_pivot];
    COL[J_pivot] = HOLD ;

    /*               reduce about pivot */
    AA[ROW[k]*n+COL[k]] = 1.0 / pivot ;
    for (j=0; j<n; j++)
    {
      if( j != k )
      {
        AA[ROW[k]*n+COL[j]] = AA[ROW[k]*n+COL[j]] * AA[ROW[k]*n+COL[k]];
      }
    }

    /*               inner reduction loop */
    for (i=0; i<n; i++)
    {
      if( k != i )
      {
        for (j=0; j<n; j++)
        {
          if( k != j )
          {
            AA[ROW[i]*n+COL[j]] = AA[ROW[i]*n+COL[j]] - AA[ROW[i]*n+COL[k]] *
                                  AA[ROW[k]*n+COL[j]] ;
          }
        }
        AA[ROW[i]*n+COL [k]] = - AA[ROW[i]*n+COL[k]] * AA[ROW[k]*n+COL[k]] ;
      }
    }
  }
  /*      end of main reduction loop */

  /*      unscramble rows */
  for (j=0; j<n; j++) {
    for (i=0; i<n; i++) {
      TEMP[COL[i]] = AA[ROW[i]*n+j];
    }
    for (i=0; i<n; i++) {
      AA[i*n+j] = TEMP[i] ;
    }
  }

  /*      unscramble columns */
  for (i=0; i<n; i++)
  {
    for (j=0; j<n; j++)
    {
      TEMP[ROW[j]] = AA[i*n+COL[j]] ;
    }
    for (j=0; j<n; j++)
    {
      AA[i*n+j] = TEMP[j] ;
    }
  }
  free(ROW);
  free(COL);
  free(TEMP);
} /* end invert */