Khaletsky method (LU decomposition)

So given a matrix I want to find the its LU decomposition.

So I could use

 > A: matrix( [1, 2, 3], [4, 5, 6], [7, 8, 9] );
 > lu_factor(A)$
 > get_lu_factors(%);

and that would give me L and U in %[2], and [%3], respectively.

However, the LU decomposition, but I am looking for it's:
     L - lower triangular matrix (not lower triangular matrix with 
unitary diagonal, which is what get_lu_factors returns)
   U - as upper triangular matrix with unitary diagonal. (not upper 
triangular without unitary diagonal, which is what get_lu_factors returns)

This is the Khaletsky method. It's not very known. A search in google 
yields very few interesting results, but yet, I would like to be able to 
compute it that way.

You can find a reference to this method here 
http://simlab.mas.bg.ac.rs/ZaSkidanje/LinAlgEq.pdf at slide 17.

So, do you have any idea how I can achieve that?

Thanks in advance.