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.