Programs for matrix inversion

The fact that C requires relatively fixed dimensions for arrays when working with multiple indices leads to developing special matrix routines for each dimension. GEOM needs different dimensions for different applications -- three dimensions for routine vector analysis in space, but four and five dimensions for the projective geometry of spheres. The following collection inverts $4\times 4$ matrices.

The method employed is Gaussian elimination with pivoting on the largest elements.

/* - - - - - - - - - M A T R I X   I N V E R S E - - - - - - - -  */

/* the subroutines required to invert a matrix of fixed dimension */

/* geomi4 - 4x4 matrix inverse */

void geomi4(z,u) double z[4][4], u[4][4]; {
int    i, j, n, ii[4], geoge4();
double f, w[4][4];
n=4;
geocm4(w,u);
geoum4(z);
for (i=0; i<n; i++) {
  ii[i]=geoge4(w,i);
  geoxr4(w,i,ii[i]);
  for (j=0; j<n; j++) {
    if (i==j) continue;
    f=-w[i][j]/w[i][i];
    geoac4(w,j,j,f,i);
    geoac4(z,j,j,f,i);
    }
  }
for (i=0; i<n; i++) geosm4(z,1.0/w[i][i],i);
for (i=0; i<n; i++) {j=n-i-1; geoxc4(z,j,ii[j]);}
}

/* geoge4 - greatest element in the nth column of 4x4 matrix */

int geoge4(p,n) double p[4][4]; int n; {double g, h; int i, m;
m=n;
g=p[n][n];
g=(g<0.0?-g:g);
for (i=n; i<4; i++) {
  h=p[i][n];
  h=(h<0.0?-h:h);
  if (h<g) continue;
  g=h; m=i;
  }
return m;
}

/* geocm4 - copy 4x4 matrix */

void geocm4(z,x) double z[4][4], x[4][4]; {int i, j;
for (i=0; i<4; i++) for (j=0; j<4; j++) z[i][j]=x[i][j];}

/* geoum4 - 4x4 unit matrix */

void geoum4(z) double z[4][4]; {int i, j;
for (i=0; i<4; i++) {
  for (j=0; j<4; j++) z[i][j]=0.0;
  z[i][i]=1.0;
  }
}

/* geoxc4 - exchange ith and jth columns of a 4x4 matrix */

void geoxc4(z,i,j) double z[4][4]; int i, j; {int k; double t;
if (i==j) return;
for (k=0; k<4; k++) {t=z[k][i]; z[k][i]=z[k][j]; z[k][j]=t;}
}

/* geoxr4 - exchange ith and jth rows of a 4x4 matrix */

void geoxr4(z,i,j) double z[4][4]; int i, j; {int k; double t;
if (i==j) return;
for (k=0; k<4; k++) {t=z[i][k]; z[i][k]=z[j][k]; z[j][k]=t;}
}

/* geoxtc4 - extract nth column from 4x4 matrix */

void geoxtc4(p,z,n) double p[4], z[4][4]; int n; {int i;
for (i=0; i<4; i++) p[i]=z[i][n];}


/* geoac4 - augment column of a 4x4 matrix */

void geoac4(z,i,j,f,k) double z[4][4], f; int i, j, k; {int l;
for (l=0; l<4; l++) z[l][i]=z[l][j]+f*z[l][k];}

/* geosm4 - multiply ith column of 4x4 matrix by a factor */

void geosm4(z,f,i) double z[4][4], f; int i; {int j;
for (j=0; j<4; j++) z[j][i]*=f;}

/* geoome4 - "one minus epsilon" for iterative */
/*	     refinement of 4x4 inverse	       */

void geoome4(z,x) double z[4][4], x[4][4]; {int i, j;
for (i=0; i<4; i++)
for (j=0; j<4; j++)
  if (i==j) z[i][j]=2.0-x[i][j]; else z[i][j]=-x[i][j];
}

/* geomp4 - product of two 4x4 matrices */

void geomp4(z,x,y) double z[4][4], x[4][4], y[4][4]; {
int i, j, k;
double u[4][4], v[4][4];
geocm4(u,x);
geocm4(v,y);
for (i=0; i<4; i++) for (j=0; j<4; j++) {
  z[i][j]=0.0;
  for (k=0; k<4; k++) z[i][j]+=u[i][k]*v[k][j];
  }
}