[comp.lang.c] Num. Recipes matrix inversion chokes on this singular matrix

ajayshah@alhena.usc.edu (Ajay Shah) (05/01/91)

I used this matrix for testing:

	2 3 4
	3 4 5
	4 5 6

which is singular.  The program chokes on it.  It works for
many other test matrices I rigged up.  That's bad because I'd
like it to detect a singular matrix and not blow up on it.

Here is my inversion code:

int d_Inverse(double **A, double **Y, int n)
{
	int *indexes;
	int i, j;
	double d;
	double *onecolumn;

	indexes = ivector(1, n);
	if (d_ludcmp(A, n, indexes, &d) == 1) return 1;
	onecolumn = dvector(1, n);

	for (i=1; i<=n; i++) {
		for (j=1; j<=n; j++) onecolumn[j] = 0.0; onecolumn[i] = 1.0;
		d_lubksb(A, n, indexes, onecolumn);
		for (j=1; j<=n; j++) Y[j][i] = onecolumn[j];
	}
	free_ivector(indexes, 1, n); free_dvector(onecolumn, 1, n);
	return 0;
}

d_ludcmp is the double precision version of the original ludcmp.
Ditto for d_lubksb.  I'm using TINY = 1.0e-100 for double
precision.  Is that a good idea?  Or should I just stick to
TINY=1.0e-20 and get more robust code?

Thanks!

	-ans.

-- 
_______________________________________________________________________________
Ajay Shah, (213)734-3930, ajayshah@usc.edu
                             The more things change, the more they stay insane.
_______________________________________________________________________________