20/02/2004 @21:02:25 ^22:13:42

void *black_hole(void);

is the worst computational astrophysics joke ever

matrix_invert

It starts with the following innocuous remarks

16/02 14:08 - 16/02 14:10
 <Timberwolf> Thankfully, with rotation matrices, to invert a rotation all
              you need to do is transpose the matrix, not invert it.
        <RjY> write a generalised matrix inverter you lazy faggot
    <EvilHRF> fazy laggot?
 <Timberwolf> I tell you what.
 <Timberwolf> You write it.
 <Timberwolf> Then I'll put it in my project.

So about a day later there you go. It was doable because I understood the algorithm from doing it in Linear Algebra back in the day, and I was able to visualise the structures and functions necessary for the calculation instead of sitting there for an hour going "uh what do I do now guys help"

So anyway what you do is this

• Matrix *M = matrix_new(N); to create a matrix of size NxN initialised to zero
• Fill it up with elements using matrix_set_elt(Matrix *M, int i, int j, Real v) where 0<=i,j<N. There is a matrix_get_elt as well. Apparently it's better to have functions to do this rather than fiddling about with the data structure in main memory, although I am very lazy when it comes to writing error checking so whatever.
• Then you can call matrix_det(M) to get its determinant (recursion, whee) and finally Matrix *N = matrix_invert(M); will create a new matrix containing the inverse. Don't forget to matrix_free(M) when you've finished using it.
• There are a few other things you can do; for example, the pleasantly trivial function that makes a copy of a matrix. Also you can do addition, scalar multiplication, matrix multiplication, dump a matrix to stdout for debugging, and so forth. Obviously all of these are in matrix.h.

There is an example main.c supplied which reads a matrix off the command line. Unpack it somewhere, type "make", then you supply the size, then the rows one after the other, as follows.

rjy@baron:~/src/matrix_invert$ ./matrix_invert 2 1 2 3 4
Let M be the matrix
(   1.000   2.000)
(   3.000   4.000)
with determinant -2.000000

Therefore M has inverse N
(  -2.000   1.000)
(   1.500  -0.500)

Check that MN is the identity: MN = 
(   1.000   0.000)
(   0.000   1.000)

It's not very fast or efficient because it uses the determinant-adjoint method. A faster implementation would use the Gaussian method, but that's harder to code. I might still do it but I don't know. Please download it and see what you think. I've never written object-like code in C before; the interface design, which I pretty much made up off the top of my head probably having been roughly inspired by things from the GLib documentation, might be a whole load of crap.

The code is placed into the public domain. Do what you like with it but don't come crying to me when your shit gets fucked up later.

In conclusion, I learned a valuable lesson about the difference between a 2d array double x[10][10]; and a pointer to a pointer double **x; Furthermore "Timberwolf" has to put my shitty code in his final year project, as well as

16/02 14:10 - 16/02 14:10
 <Timberwolf> And we can have a credits screen that flashes in all-caps, 
              blue and yellow, "MATRIX INVERTER BY RjY!!!!!"