/* * This code was (mostly) written by Ken Turkowski, who said: * * Oh, that. I wrote it in college the first time. It's open source - I think I * posted it after seeing so many people solve equations by inverting matrices * by computing minors naïvely. * -Ken * * The views represented here are mine and are not necessarily shared by * my employer. Ken Turkowski turk@apple.com Immersive Media Technologist http://www.worldserver.com/turk/ Apple Computer, Inc. 1 Infinite Loop, MS 302-3VR Cupertino, CA 95014 */ /* Matinv() inverts the matrix A using LU decomposition. Arguments: * A - the (n x n) matrix to be inverted * Ainv - the (n x n) inverted matrix * n - the order of the matrices A and Ainv */ #include <stdlib.h> #include <render.h> extern int lu_decompose(double **a, int n); extern void lu_solve(double *x, double *b, int n); int matinv(double **A, double **Ainv, int n) { register int i, j; double *b, temp; /* Decompose matrix into L and U triangular matrices */ if (lu_decompose(A, n) == 0) return(0); /* Singular */ /* Invert matrix by solving n simultaneous equations n times */ b = N_NEW(n,double); for (i = 0; i < n; i++) { for(j = 0; j < n; j++) b[j] = 0.0; b[i] = 1.0; lu_solve(Ainv[i], b, n); /* Into a row of Ainv: fix later */ } free(b); /* Transpose matrix */ for (i = 0; i < n; i++) { for (j = 0; j < i; j++) { temp = Ainv[i][j]; Ainv[i][j] = Ainv[j][i]; Ainv[j][i] = temp; } } return(1); }