cs-ind04@penelope.Oswego.EDU (Jason E. Forbes) (08/03/90)
Before I reinvent the wheel, does anyone have a C function that can find the inverse of a 4 by 4 matrix? Thanks, Jason Forbes cs-ind04@penelope.oswego.edu cs-ind04@oswego.oswego.edu Veterans Administration Hospital, Research Department Syracuse, NY -----------------------------------------------------
dodson@convex.COM (Dave Dodson) (08/03/90)
In article <1990Aug3.021200.21023@oswego.Oswego.EDU> cs-ind04@penelope.Oswego.EDU (Jason E. Forbes) writes: >Before I reinvent the wheel, does anyone have a C function that can >find the inverse of a 4 by 4 matrix? First, I would ask why you want to invert a matrix. While papers and books extensively use the notation A-inverse * b to denote the solution of the system of linear equations Ax = b, it is better to solve the system directly than to compute the inverse and multiply it times the right hand side. By _better_, I mean faster and more accurate. So if you are wanting to solve linear equations, do it directly. ---------------------------------------------------------------------- Dave Dodson dodson@convex.COM Convex Computer Corporation Richardson, Texas (214) 497-4234
ted@nmsu.edu (Ted Dunning) (08/04/90)
From: cs-ind04@penelope.Oswego.EDU (Jason E. Forbes)
Newsgroups: comp.sources.wanted,sci.math
Before I reinvent the wheel, does anyone have a C function that can
find the inverse of a 4 by 4 matrix?
...
the real answer to your question is to buy yourself a copy of
numerical recipes in c. just for fun, though, i revved maple up and
had it spit out the inversion in fortran and then used emacs to change
it into c. the total time including cursory debugging was about 30
minutes, and the resulting code obviously would provide a major
opportunity for a fancy optimizer, but would probably be a disaster on
a vectorizer. but, of course, 4x4 inversion isn't what you would be
doing on a vector machine anyway. no guarantees for numerical
stability, but anything at all should work pretty well for such small
matrices.
it would be an interesting question to to an operation count on the
resulting compiled code and comparing the result to what a looping
version would do. it would also be interesting to compare the number
of operations that a common subexpression optimizer could save.
so anyway, have fun with this.
/* invert a 4x4 matrix stored in a and store the result in b. both a
and b must be preallocated. this code was generated automatically
in fortran and then translated to c, so it should be thoroughly
tested before use */
invert_4_4(a,b)
double **a, **b;
{
double t1, t2, t4, t5, t6, t7, t8, t9, t11, t12, t14, t15, t17;
double t19, t22, t23, t26, t27, t29, t30, t32, t33, t48, t50, t53, t56;
double t58, t59, t61, t62, t66, t69, t71, t73, t76, t79, t84, t116, t118;
double t121, t124, t126, t128, t133, t135, t138, t141, t143, t145, t151;
double t153, t159, t161, t166, t168, t174, t176, t187, t189, t212, t214;
t1 = a[2][2]*a[3][3];
t2 = a[0][1]*t1;
t4 = a[2][3]*a[3][2];
t5 = a[0][1]*t4;
t6 = a[0][2]*a[3][3];
t7 = a[2][1]*t6;
t8 = a[0][3]*a[3][2];
t9 = a[2][1]*t8;
t11 = a[0][2]*a[2][3];
t12 = a[3][1]*t11;
t14 = a[0][3]*a[2][2];
t15 = a[3][1]*t14;
t17 = a[1][1]*t1;
t19 = a[1][1]*t4;
t22 = a[1][2]*a[3][3];
t23 = a[2][1]*t22;
t26 = a[1][3]*a[3][2];
t27 = a[2][1]*t26;
t29 = a[1][2]*a[2][3];
t30 = a[3][1]*t29;
t32 = a[1][3]*a[2][2];
t33 = a[3][1]*t32;
t48 = a[0][1]*t22;
t50 = a[0][1]*t26;
t53 = a[1][1]*t6;
t56 = a[1][1]*t8;
t58 = a[0][2]*a[1][3];
t59 = a[3][1]*t58;
t61 = a[0][3]*a[1][2];
t62 = a[3][1]*t61;
t66 = a[0][1]*t29;
t69 = a[0][1]*t32;
t71 = a[1][1]*t11;
t73 = a[1][1]*t14;
t76 = a[2][1]*t58;
t79 = a[2][1]*t61;
t84 = 1/(a[0][0]*t17-a[0][0]*t19-a[0][0]*t23+a[0][0]*t27+a[0][0]*t30-a[0][0]*t33-a[1][0]*t2+a[1][0]*t5+a[1][0]*t7-a[1][0]*t9-a[1][0]*t12+a[1][0]*t15+a[2][0]*t48-a[2][0]*t50-a[2][0]*t53+a[2][0]*t56+a[2][0]*t59-a[2][0]*t62-a[3][0]*t66+a[3][0]*t69+a[3][0]*t71-a[3][0]*t73-a[3][0]*t76+a[3][0]*t79);
t116 = a[2][1]*a[3][2];
t118 = a[2][2]*a[3][1];
t121 = a[0][1]*a[3][2];
t124 = a[0][2]*a[3][1];
t126 = a[0][1]*a[2][2];
t128 = a[0][2]*a[2][1];
t133 = a[2][1]*a[3][3];
t135 = a[2][3]*a[3][1];
t138 = a[0][1]*a[3][3];
t141 = a[0][3]*a[3][1];
t143 = a[0][1]*a[2][3];
t145 = a[0][3]*a[2][1];
t151 = a[1][1]*a[2][2];
t153 = a[1][2]*a[2][1];
t159 = a[0][1]*a[1][2];
t161 = a[0][2]*a[1][1];
t166 = a[1][1]*a[3][3];
t168 = a[1][3]*a[3][1];
t174 = a[0][1]*a[1][3];
t176 = a[0][3]*a[1][1];
t187 = a[1][1]*a[2][3];
t189 = a[1][3]*a[2][1];
t212 = a[1][1]*a[3][2];
t214 = a[1][2]*a[3][1];
b[0][1] = (-t2+t5+t7-t9-t12+t15)*t84;
b[0][2] = -(-t48+t50+t53-t56-t59+t62)*t84;
b[1][2] = -(a[0][0]*t22-a[0][0]*t26-a[1][0]*t6+a[1][0]*t8+a[3][0]*t58-a[3][0]*t61)*t84;
b[1][3] = -(-a[0][0]*t29+a[0][0]*t32+a[1][0]*t11-a[1][0]*t14-a[2][0]*t58+a[2][0]*t61)*t84;
b[3][1] = (a[0][0]*t116-a[0][0]*t118-a[2][0]*t121+a[2][0]*t124+a[3][0]*t126-a[3][0]*t128)*t84;
b[2][1] = -(a[0][0]*t133-a[0][0]*t135-a[2][0]*t138+a[2][0]*t141+a[3][0]*t143-a[3][0]*t145)*t84;
b[3][3] = (a[0][0]*t151-a[0][0]*t153-a[1][0]*t126+a[1][0]*t128+a[2][0]*t159-a[2][0]*t161)*t84;
b[2][2] = (a[0][0]*t166-a[0][0]*t168-a[1][0]*t138+a[1][0]*t141+a[3][0]*t174-a[3][0]*t176)*t84;
b[2][0] = (a[1][0]*t133-a[1][0]*t135-a[2][0]*t166+a[2][0]*t168+a[3][0]*t187-a[3][0]*t189)*t84;
b[1][0] = -(a[1][0]*t1-a[1][0]*t4-a[2][0]*t22+a[2][0]*t26+a[3][0]*t29-a[3][0]*t32)*t84;
b[0][3] = -(t66-t69-t71+t73+t76-t79)*t84;
b[3][2] = -(a[0][0]*t212-a[0][0]*t214-a[1][0]*t121+a[1][0]*t124+a[3][0]*t159-a[3][0]*t161)*t84;
b[0][0] = (t17-t19-t23+t27+t30-t33)*t84;
b[2][3] = -(a[0][0]*t187-a[0][0]*t189-a[1][0]*t143+a[1][0]*t145+a[2][0]*t174-a[2][0]*t176)*t84;
b[3][0] = -(a[1][0]*t116-a[1][0]*t118-a[2][0]*t212+a[2][0]*t214+a[3][0]*t151-a[3][0]*t153)*t84;
b[1][1] = (a[0][0]*t1-a[0][0]*t4-a[2][0]*t6+a[2][0]*t8+a[3][0]*t11-a[3][0]*t14)*t84;
}
/* multiply a times b to get c. c must be preallocated. code
generated in fortran by maple and then converted semi-automatically
to c. it may well be buggy. */
mul_4_4(a,b,c)
double **a, **b, **c;
{
c[0][1] = a[0][0]*b[0][1]+a[0][1]*b[1][1]+a[0][2]*b[2][1]+a[0][3]*b[3][1];
c[0][2] = a[0][0]*b[0][2]+a[0][1]*b[1][2]+a[0][2]*b[2][2]+a[0][3]*b[3][2];
c[1][2] = a[1][0]*b[0][2]+a[1][1]*b[1][2]+a[1][2]*b[2][2]+a[1][3]*b[3][2];
c[1][3] = a[1][0]*b[0][3]+a[1][1]*b[1][3]+a[1][2]*b[2][3]+a[1][3]*b[3][3];
c[3][1] = a[3][0]*b[0][1]+a[3][1]*b[1][1]+a[3][2]*b[2][1]+a[3][3]*b[3][1];
c[2][1] = a[2][0]*b[0][1]+a[2][1]*b[1][1]+a[2][2]*b[2][1]+a[2][3]*b[3][1];
c[3][3] = a[3][0]*b[0][3]+a[3][1]*b[1][3]+a[3][2]*b[2][3]+a[3][3]*b[3][3];
c[2][2] = a[2][0]*b[0][2]+a[2][1]*b[1][2]+a[2][2]*b[2][2]+a[2][3]*b[3][2];
c[2][0] = a[2][0]*b[0][0]+a[2][1]*b[1][0]+a[2][2]*b[2][0]+a[2][3]*b[3][0];
c[1][0] = a[1][0]*b[0][0]+a[1][1]*b[1][0]+a[1][2]*b[2][0]+a[1][3]*b[3][0];
c[0][3] = a[0][0]*b[0][3]+a[0][1]*b[1][3]+a[0][2]*b[2][3]+a[0][3]*b[3][3];
c[3][2] = a[3][0]*b[0][2]+a[3][1]*b[1][2]+a[3][2]*b[2][2]+a[3][3]*b[3][2];
c[0][0] = a[0][0]*b[0][0]+a[0][1]*b[1][0]+a[0][2]*b[2][0]+a[0][3]*b[3][0];
c[2][3] = a[2][0]*b[0][3]+a[2][1]*b[1][3]+a[2][2]*b[2][3]+a[2][3]*b[3][3];
c[3][0] = a[3][0]*b[0][0]+a[3][1]*b[1][0]+a[3][2]*b[2][0]+a[3][3]*b[3][0];
c[1][1] = a[1][0]*b[0][1]+a[1][1]*b[1][1]+a[1][2]*b[2][1]+a[1][3]*b[3][1];
}
/* simple cursory test of the inversion and multiplication routines */
main()
{
double **a,**b,**c;
void *calloc();
int i,j;
a = calloc(4,sizeof(a[0]));
b = calloc(4,sizeof(b[0]));
c = calloc(4,sizeof(c[0]));
for (i=0;i<4;i++) {
a[i] = calloc(4,sizeof(a[0][0]));
b[i] = calloc(4,sizeof(b[0][0]));
c[i] = calloc(4,sizeof(c[0][0]));
for (j=0;j<4;j++) {
a[i][j] = random()%1000;
}
}
invert_4_4(a,b);
mul_4_4(a,b,c);
for (i=0;i<4;i++) {
for (j=0;j<4;j++) {
printf("%19.17lf ", c[i][j]);
}
printf("\n");
}
}
--
Offer void except where prohibited by law.