/* inverseg.c  Gausian elimination version that keeps going when singular */
/*             this may be called the general inverse, that can return    */
/*             rows and columns of all zeros                              */
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#define abs(a) ((a)<0?(-(a)):(a))

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

  degree = 0; /* degree of singularity */
  ROW = (int *)calloc(n, sizeof(int));
  COL = (int *)calloc(n, sizeof(int));
  TEMP = (double *)calloc(n, sizeof(double));
  memcpy(AA,A,n*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]] ;
        }
      }
    }
    HOLD = ROW[k]; /* could be null interchange, yet may be needed */
    ROW[k]= ROW[I_pivot];
    ROW[I_pivot] = HOLD ;
    HOLD = COL[k];
    COL[k]= COL[J_pivot];
    COL[J_pivot] = HOLD ;
    if(abs(pivot) < 1.0E-24) /* may be set smaller */
    {
      degree++;
      for(ip=0; ip<n; ip++) /* general inverse of singular matirx */
      {
        AA[ROW[ip]*n+COL[k]] = 0.0;
        AA[ROW[k]*n+COL[ip]] = 0.0;
      }
    }
    else /* not singular here */
    {
      /*               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 singular detection */
  }
  /*      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] ;
    }
  }
  if(degree>0)
    printf("warning inverseg detected singularity degree %d of %d\n",
           degree, n);
  free(ROW);
  free(COL);
  free(TEMP);
} /* end inverse */