Decompose a matrix A as A.P = Q.R using Householder transforms.
As suggested in the P. Lascaux and R. Theodor book
Analyse numérique matricielle appliquée à
l'art de l'ingénieur (Masson, 1986), instead of representing
the Householder transforms with uk unit vectors such that:
Hk = I - 2uk.ukt
we use k non-unit vectors such that:
Hk = I - betakvk.vkt
where vk = ak - alphak ek.
The betak coefficients are provided upon exit as recomputing
them from the vk vectors would be costly.
This decomposition handles rank deficient cases since the tranformations
are performed in non-increasing columns norms order thanks to columns
pivoting. The diagonal elements of the R matrix are therefore also in
non-increasing absolute values order.