bug in lu_substitute

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:    +1 973 360 7127       Email:  [email protected]
Fax:    +1 973 360 8178       Recycle yourself -- be an organ donor

_______________________________________________
ublas mailing list
[email protected]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [email protected]