1 2 SUBROUTINE SVD(NM,M,N,A,W,MATU,U,MATV,V,IERR,RV1) 3C***BEGIN PROLOGUE SVD 4C***REFER TO EISDOC 5C 6C This subroutine is a translation of the ALGOL procedure SVD, 7C NUM. MATH. 14, 403-420(1970) by Golub and Reinsch. 8C HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971). 9C 10C This subroutine determines the singular value decomposition 11C T 12C A=USV of a REAL M by N rectangular matrix. Householder 13C bidiagonalization and a variant of the QR algorithm are used. 14C 15C On Input 16C 17C NM must be set to the row dimension of two-dimensional 18C array parameters as declared in the calling program 19C dimension statement. Note that NM must be at least 20C as large as the maximum of M and N. 21C 22C M is the number of rows of A (and U). 23C 24C N is the number of columns of A (and U) and the order of V. 25C 26C A contains the rectangular input matrix to be decomposed. 27C 28C MATU should be set to .TRUE. if the U matrix in the 29C decomposition is desired, and to .FALSE. otherwise. 30C 31C MATV should be set to .TRUE. if the V matrix in the 32C decomposition is desired, and to .FALSE. otherwise. 33C 34C On Output 35C 36C A is unaltered (unless overwritten by U or V). 37C 38C W contains the N (non-negative) singular values of A (the 39C diagonal elements of S). They are UNORDERED. If an 40C error exit is made, the singular values should be correct 41C for indices IERR+1,IERR+2,...,N. 42C 43C U contains the matrix U (orthogonal column vectors) of the 44C decomposition if MATU has been set to .TRUE. Otherwise 45C U is used as a temporary array. U may coincide with A. 46C If an error exit is made, the columns of U corresponding 47C to indices of correct singular values should be correct. 48C 49C V contains the matrix V (orthogonal) of the decomposition if 50C MATV has been set to .TRUE. Otherwise V is not referenced. 51C V may also coincide with A if U is not needed. If an error 52C exit is made, the columns of V corresponding to indices of 53C correct singular values should be correct. 54C 55C IERR is set to 56C Zero for normal return, 57C K if the K-th singular value has not been 58C determined after 30 iterations. 59C 60C RV1 is a temporary storage array. 61C 62C CALLS PYTHAG(A,B) for sqrt(A**2 + B**2). 63C 64C Questions and comments should be directed to B. S. Garbow, 65C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 66C ------------------------------------------------------------------ 67C***ROUTINES CALLED PYTHAG 68C***END PROLOGUE SVD 69 70 71