Num. Recipes matrix inversion chokes on this singular matrix
Ajay Shah
ajayshah at alhena.usc.edu
Wed May 1 09:34:05 AEST 1991
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 at usc.edu
The more things change, the more they stay insane.
_______________________________________________________________________________
More information about the Comp.lang.c
mailing list