// With gcc compile thus:
// gcc invert.c -std=gnu99 -Wall
#include <math.h>
/* Inversion uses twice the necessary storage and does nearly twice the
   necessary arithmetic. Many of the results are known to be zero
   or one. These are computed and stored as well. */

// slv below solves a set of n simultaneous linear equations
// for m different sets of n constants.
// The calculations overwrite all of the inputs
// and the answers are found where the constants were provided.
int slv(int n, int m, double A[n][n+m]); // prototype
int slv(int n, int m, double a[n][n+m]){
  for(int i=0; i<n; i++) {double big=0; int im=n;
    for(int j=i; j<n; j++) if(big<fabs(a[j][i])) {im=j; big=fabs(a[j][i]);}
    if(im==n) return 1;
    if(im!=i) for (int j=i; j<n+m; j++)
      {double q = a[i][j]; a[i][j]=a[im][j]; a[im][j]=q;}
    for(int j=i+1; j<n; j++) {double g = a[j][i]/a[i][i];
      for(int k=0; k<n+m; k++) a[j][k] -= g*a[i][k];}}
  for(int i=0; i<n; i++) {double g = 1.0/a[i][i];
    for(int j=i; j<n+m; j++) a[i][j]*=g; 
    for(int j=0; j<i; j++) {double g=a[j][i];
      for(int k=0; k<n+m; k++) a[j][k] -=g*a[i][k];}}
return 0;}

// invert below deposits the inverse of square matrix A into matrix B.
// A and B may overlap. A is unmodified otherwise.
int invert(int n, double A[n][n], double B[n][n]); // prototype
int invert(int n, double A[n][n], double B[n][n]){
  double a[n][2*n];
  for(int i=0; i<n; i++) for(int j=0; j<n; j++)
    {a[j][i] = A[j][i]; a[j][n+i]= i==j;}
  int r = slv(n, n, a);
  for(int i=0; i<n; i++) for(int j=0; j<n; j++) B[i][j]=a[i][j+n];
  return r;}

// solve below solves a set of linear equations:
// for k [0<=k<n] (sum[0<=j<n] A[k][j]*x[j] = B[k])
// The answer, x[k], overwrites B[k].
int solve(int n, double A[n][n], double B[n]);
int solve(int n, double A[n][n], double B[n]){
  double a[n][n+1];
  for(int i=0; i<n; i++) {
    for(int j=0; j<n; j++) a[i][j] = A[i][j]; a[i][n]= B[i];}
  int r = slv(n, 1, a);
  for(int i=0; i<n; ++i) B[i] = a[i][n];
  return r;}

#if 1
// main below is merely a simple demonstration and test.
#define n 10
#include <stdio.h>
#include <stdlib.h>
static void pr(char * m, double x[n][n])
   {for(int i=0; i<n; ++i) {for (int j=0; j<n; ++j) printf("%16.12f ", x[i][j]);
      printf(" <- %s\n", m);} printf("\n");}

int main(){if(1){double a[n][n], x[n][n], p[n][n];
  for(int i=0; i<n; ++i) for (int j=0; j<n; ++j) a[i][j] = random()/1.e9;
  invert(n, a, x);
  for(int i=0; i<n; ++i) for (int j=0; j<n; ++j)
    {p[i][j] = 0; for(int k=0; k<n; ++k) p[i][j] += a[i][k]*x[k][j];}
  pr("a", a); pr("x", x); pr("p", p);}
else {double a[][3] = {{3, 5, 4}, {2, 6, 6}, {3, -2, 5}};
  double v[] = {4, 2.6, -3.1};
  solve(3, a, v);
  for(int j=0; j<3; ++j){double t = 0;
    for(int k=0; k<3; ++k) t += v[k]*a[j][k];
    printf("v=%16.12f  ", t);} printf("\n");}
return 0;}
#endif