Xref: utzoo comp.lang.c:17765 misc.wanted:4636
Path: utzoo!utgpu!jarvis.csri.toronto.edu!mailrus!ncar!ames!haven!adm!smoke!gwyn
From: gwyn@smoke.BRL.MIL (Doug Gwyn)
Newsgroups: comp.lang.c,misc.wanted
Subject: Re: Need matrix inversion C routine.
Keywords: C matrix invert
Message-ID: <10087@smoke.BRL.MIL>
Date: 21 Apr 89 16:46:46 GMT
References: <5785@cbnews.ATT.COM>
Reply-To: gwyn@brl.arpa (Doug Gwyn)
Organization: Ballistic Research Lab (BRL), APG, MD.
Lines: 48

In article <5785@cbnews.ATT.COM> wkb@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.