Singular value decomposition in extended double precision arithmetic

A well-known and successful algorithm to compute the singular value decomposition (SVD) of a matrix was published by Golub and Reinsch ( Numer. Math. 14:403–420, 1970 ), together with an implementation in Algol. We give an updated implementation in extended double precision arithmetic in the C programming language. Extended double precision is native for Intel x86 processors and provides improved accuracy at full hardware speed. The complete program for computing the SVD is listed. Additionally, a comprehensive explanation of the original algorithm of Golub and Reinsch ( Numer. Math. 14:403–420, 1970 ) is given at an elementary level without referring to the more general results of Francis ( Comput. J. 4:265–271, 1961 , 1962 ).