Re: Free42 matrixmatrix division bug Message #3 Posted by Thomas Okken on 10 Nov 2004, 10:16 a.m., in response to message #2 by hugh steers
According to the HP Journal article (May 1983, page 3334),
if a calculated diagonal element of U, which we call a pivot, is found to be zero during the LU decomposition, rather than aborting the matrix calculation and reporting the input matrix to be singular, the HP15C replaces the zero pivot by a small positive number and continues with the calculation. This number is usually small compared to the rounding errors involved in the calculations. Specifically, it will be about 10^10 times the largest absolute value of any element in that column of the original matrix.
The article also points out why testing for singularity by looking for zero pivots is unreliable, and why continuing despite finding zero pivots is useful in certain applications.
When I copied the LU decomposition code from Numerical Recipes in Fortran, I noticed something similar; I replaced it with a fatal error because I didn't like what the NR code did: it uses an arbitrary small value, which is not in any way guaranteed to have the right scale, and there's no explanation why this is done, either.
I think I'll use the HP15C approach instead, now. I'll have to study the issue a bit more but when I make this change to lu_decomp_r() and lu_decomp_c(), that would probably also be a good moment to add proper overflow checking and reporting to those routines.
UPDATE  I just checked and the 15C and 42S have another interesting property: they won't even report that a matrix is singular even if it contains nothing but zeroes. This is a different case than what is mentioned in the NR book and the HPJ article. Apparently, a zero column in the original matrix simply results in a zero column with a large number (1e99 on the 15C, 9.999...e499 on the 42S) on the diagonal position. Just wondering: does anybody know what that is for?
Edited: 10 Nov 2004, 11:19 a.m.
