I am trying to invert some 10X10 matrices using the subroutine "gaussj" (a Gauss-Jordan elimination method), but I am not able to recover the identity matrix with any accuracy. I am able to compile and build the NR example program "xgaussj" and reproduce the results of the NR sample data. It appears to be more accurate with matrices smaller than 5X5, but still not at the level I am satisfied with. There is not much information on the limitations and applicability of this subroutine. Does any one know if there is a size limitation in using this routine? I would be surprised if there was. I have been able to recover the identity matrix using other numerical software. I've attached my original matrix, inverse using guassj, and actual, confirmed inverse.

Thanks for the help, Charlie

Hi Charlie,

The matrix you are working with is fairly ill-conditioned (close to singular). You can see this from the fact that the typical entries of the inverse are so much bigger than the entries of the original matrix. Run the example in double precision and it will work fine.

Saul Teukolsky

Saul,

Thanks for the info. I try DP and let you know how it turns out.

Charlie Weil