**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.