[SOLVED] What does lu_factorize return?

What does lu_factorize return?

boost::number::ublas contains the M::size_type lu_factorize(M& m) function. Its name suggests that it performs the LU decomposition of a given matrix m, i.e. should produce two matrices that m = L*U. There seems to be no documentation provided for this function.

It is easy to deduce that it returns 0 to indicate successful decomposition, and a non-zero value when the matrix is singular. However, it is completely unclear where is the result. Taking the matrix by reference suggests that it works in-place, however it should produce two matrices (L and U) not one. So what does it do?

Solution

There is no documentation in boost, but looking at the documentation of SciPy's lu_factor one can see, that it's not uncommon to return one result for the LU decomposition.

This is enough, because in a typical approach to LU decomposition, L's diagonal consists of ones only, as presented in this answer from Mathematics, for example.

So, it is possible to fit both L and U into one matrix, putting L in result's lower part, omitting the diagonal (which is assumed to contain only ones), and U in the upper part. For example, for a 3x3 problem the result is:

    u11 u12 u13
m = l21 u22 u23
    l31 l32 u33

which implies:

     1    0   0
L =  l21  1   0
     l31  l32 1

and

    u11 u12 u13
U = 0   u22 u23
    0   0   u33

Inspecting boost's void lu_substitute(const M& m, vector_expression<E>& e) function, from the same namespace seems to confirm this. It solves the equation LUx = e, where both L and U are contained in its m argument in two steps.

First solve Lz = e for z, where z = Ux, using lower part of m:

inplace_solve(m, e, unit_lower_tag ());

then, having computed z = Ux (with e modified in place), Ux = e can be solved, using upper part of m:

inplace_solve(m, e, upper_tag ());

inplace_solve is mentioned in the documentation, and it:

Solves a system of linear equations with triangular form, i.e. A is triangular.

So everything seems to make sense.

Also as mentioned by @Catriel in the comments, beside the output argument, the function returns non-zero value if the matrix is not singular, and zero otherwise.