Need matrix inversion C routine.

[Doug Gwyn]

In article <5785 at cbnews.ATT.COM> wkb at cbnews.ATT.COM (Wm. Keith Brummett) writes:
>    I have need of a small, fast routine in C language that will invert
>    matrices of order <= 4.

/* The following do not attempt to handle singularities or ill-conditioning. */
void
Inv1( double a[1][1], double b[1][1] )
	{
	b[0][0] = 1.0 / a[0][0];
	}
void
Inv2( double a[2][2], double b[2][2] )
	{
	register double	d = a[0][0] * a[1][1] - a[0][1] * a[1][0];
	b[0][0] =  a[1][1] / d;
	b[0][1] = -a[1][0] / d;
	b[1][0] = -a[0][1] / d;
	b[1][1] =  a[0][0] / d;
	}
void
Inv3( double a[3][3], double b[3][3] )
	{
	double		m00 = a[1][1] * a[2][2] - a[2][1] * a[1][2];
	double		m01 = a[1][2] * a[2][0] - a[2][2] * a[1][0];
	double		m02 = a[1][0] * a[2][1] - a[2][0] * a[1][1];
	register double	d = a[0][0] * m00 + a[0][1] * m01 + a[0][2] * m02;
	b[0][0] = m00 / d;
	b[0][1] = m01 / d;
	b[0][2] = m01 / d;
	b[1][0] = (a[2][1] * a[0][2] - a[0][1] * a[2][2]) / d;
	b[1][1] = (a[2][2] * a[0][0] - a[0][2] * a[2][0]) / d;
	b[1][2] = (a[2][0] * a[0][1] - a[0][0] * a[2][1]) / d;
	b[2][0] = (a[0][1] * a[1][2] - a[1][1] * a[0][2]) / d;
	b[2][1] = (a[0][2] * a[1][0] - a[1][2] * a[0][0]) / d;
	b[2][2] = (a[0][0] * a[1][1] - a[1][0] * a[0][1]) / d;
	}
void
Inv4( double a[4][4], double b[4][4] )
	{
	/* XXX -- you provide this yourself, I'm getting tired */
	}

>    BTW, can anyone tell me why it is that every language except C seems
>    to have standard subroutines to do this?

Very few languages have standard matrix inversion functions.
APL is about the only one I know of.
Certainly not Fortran or Pascal.