Bidiagonalization and SVD of matrix I can't find a single solid explanation of how to implement this -- whitepapers too detailed/confusing. Closest I came to an answer was this:
http://www.hep.ucl.ac.uk/~bino/libbpm/doc/pro/html/gsl__linalg_8c-source.html
see:
[1] gsl_linalg_bidiag_decomp
[2] gsl_linalg_SV_decomp (which calls [1])
Which perform it, but really, really obfusticate the process underneath.
 A: Golub-Kahan Bidiagonalisation
In this process householder reflectors are applied alternatively on the left and then the right. The $i^{\text{th}}$ left reflector introduces zeros below the diagonal in the $i^{\text{th}}$ column. The $i^{\text{th}}$ right reflector introduces zeros to the right of the first super-diagonal in the $i^{\text{th}}$ row.
In software packages I suspect they use a mixture of this Golub-Kahan bidiagonalisation and a process called Lawson-Hanson-Chan (LHC) bidiagonalisation depending on the size of the matrix.
Computing the SVD
The first phase of computing the SVD is bidiagonalising the matrix.  Then the SVD of the bidiagonal matrix is determined by a process very similar to the QR algorithm. This process is described in Golub and Kahn, "Calculating the singular values and pseudo-inverse of a matrix" (1960).
Since this paper there have been some alterations to provide better accuracy when the singular values are small, see Demmel and Kahan: "Accurate singular values of bidiagonal matrices" (1990)
Having been lectured by N. Trefethen on this very subject he briefly mentioned a divide-and-conquer type algorithm was now state-of-the-art though I don't know much about the details.
