#include <math.h>
#include <malloc.h>


/*
 * Invert the matrix xa.
 *
 * A should be dimensioned (n x n).
 *
 * Also, A _must_ be symmetric. The algorithm relies
 * upon this being true.
 */

fast_matrix_invert(xa,n)
	double **xa;
	int n;
{
  register int i,j,k;
  register double **a,**pj;
  register double *v,*pi,*pk;
  register double scale;
  

  
  /* Okay, I have seen so-called optimizing compilers screw      */
  /* optimizations like this up... that's why this line is here. */
  /* In general, it's a good idea not to leave everything to the */
  /* compiler... because unless you wrote the compiler, you      */
  /* shouldn't assume that it does all those nice things you     */
  /* might think it does  (even with the -O option on...)        */
  
  a = xa;
  
  /* V can be thought of as an extra column for the elimination. */

  v = (double *) malloc((long)n * sizeof(double));
  if (!v) return(-1);
  

  /* First work down and eliminate. */

  for (j=0;j<n;j++) {		
    
    /* Prepare appended column. */
/*
    for (k=j;k<n;k++) v[k] = 0;
*/
    for (pk=(v+(j<<3));pk<(v+(n<<3));pk++) *pk = 0;
    
    /* Second, normalize jth row w.r.t. jth column */
    scale = 1 / a[j][j];
    for (i=0;i<n;i++) a[j][i] *= scale;
    v[j] =  scale;

    /* Third, eliminate */
    for (i=j+1; i<n; i++) {
      scale = a[i][j];
      for (k=0;k<n;k++) {	
	a[i][k] -= scale * a[j][k];
      }
      v[i] -= scale * v[j];		/* Do the appended column, too. */
    }
    /* Now move the appended column onto the eliminated column */
    for (i=j;i<n;i++) {
      a[i][j] = v[i];
    }
  }

  /* Now work up backwards killing off columns in the original matrix. */
  /* But we only have to calculate values in the inverted matrix. */

  for (j=n-1;j>=0;j--) {

    for (i=j-1;i>=0;i--) {
      scale = a[i][j];
      for (k=i;k>=0;k--) {
	a[i][k] -= scale * a[j][k];
      }
      a[i][j] = a[j][i];
    }
  }

  free(v);
}
