https://bugzilla.redhat.com/show_bug.cgi?id=2253356 --- Comment #7 from Benson Muite <benson_muite@xxxxxxxxxxxxx> --- 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 -- You are receiving this mail because: You are always notified about changes to this product and component You are on the CC list for the bug. https://bugzilla.redhat.com/show_bug.cgi?id=2253356 Report this comment as SPAM: https://bugzilla.redhat.com/enter_bug.cgi?product=Bugzilla&format=report-spam&short_desc=Report%20of%20Bug%202253356%23c7 -- _______________________________________________ package-review mailing list -- package-review@xxxxxxxxxxxxxxxxxxxxxxx To unsubscribe send an email to package-review-leave@xxxxxxxxxxxxxxxxxxxxxxx Fedora Code of Conduct: https://docs.fedoraproject.org/en-US/project/code-of-conduct/ List Guidelines: https://fedoraproject.org/wiki/Mailing_list_guidelines List Archives: https://lists.fedoraproject.org/archives/list/package-review@xxxxxxxxxxxxxxxxxxxxxxx Do not reply to spam, report it: https://pagure.io/fedora-infrastructure/new_issue
