https://bugzilla.redhat.com/show_bug.cgi?id=2253356
--- Comment #7 from Benson Muite benson_muite@emailplus.org --- The relevant section of code is: https://github.com/nwh/lusol/blob/master/src/lusol8b.f#L739-L786 A similar section of code is: https://github.com/nwh/lusol/blob/master/src/lusol.f90#L6479-L6528
In lusol.f90 in routine lu8rpc jrep corresponds to the column being replaced in A and krep is the modified column in U, A=LU. The main steps are i) declare singular if all entries in the column krep are zero ii) if not all entries in column krep are zero, declare singular if the index of the column number in U being modified is not the same as jrep, the column being replaced in A iii) if not declared singular, declare singular if the diagonal entry of U corresponding to the column index is less than some threshold, otherwise it is not singular.
In this case singular is approximately singular, so entry is close enough to zero for numerical problems.
In lusol8b.f the routine lu8mod updates the LU factorization A = L*U when the m by n matrix A is subjected to a rank-one modification to become A + beta * v * w(transpose) kfirst to klast are the modified columns in U. i) declare singular if all entries in the column klast are zero ii) jrep is used here but not declared. Would expect if not all entries in column klast are zero, declare singular if the index of the last column number in U being modified klast is not the same as the corresponding location of the modified column in A - this may not be klast due to permutations. A loop is used to find krep in lusol.f90 : https://github.com/nwh/lusol/blob/master/src/lusol.f90#L6333-L6336 Probably something similar is needed such as: jrep = n + 1
600 jrep = jrep - 1 if (iq(jrep) /= klast) go to 600 This differs slightly from the earlier loop to find klast: https://github.com/nwh/lusol/blob/master/src/lusol8b.f#L608-L612 iii) if not declared singular, declare singular if the diagonal entry of U corresponding to the column klast is less than some threshold, otherwise it is not singular.
Description of a portion of the work: https://apps.dtic.mil/sti/tr/pdf/ADA169255.pdf