[ublas] bug in lu_substitute

Subject: [ublas] bug in lu_substitute
From: David Applegate (david_at_[hidden])
Date: 2013-02-08 17:52:30

Hello,
There is a bug in
lu_substitute (MV &mv, const M &m, const permutation_matrix<PMT, PMA> &pm)

This is the lu_substitute for solving y' A = b', equivalent to solving
A' y = b.

The problem is demonstrated by the attached test case.  The cause is
that lu_substitute applies the permutation pm in the same way it does
for the transposed substitute:

(from lu.hpp lines 343-348):
    template<class MV, class M, class PMT, class PMA>
    void lu_substitute (MV &mv, const M &m, const permutation_matrix<PMT, PMA> &pm) {
        swap_rows (pm, mv);
        lu_substitute (mv, m);
    }

However, when solving y' A = b' using the LU factorization of A, it
should apply the inverse permutation after the substitution instead.

A corrected version is also contained in the attached test case, in
the function alt_lu_substitute, which uses functions unswap_rows for
the inverse permutation.

Thank you,

David Applegate               AT&T Labs Research
Tel:    [phone]       Email:  david_at_[hidden]
Fax:    [phone]       Recycle yourself -- be an organ donor