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?

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.