C....*...1.........2.........3.........4.........5.........6.........7.*.......8 C C This is a conversion of Algorithm 358 from COMPLEX to REAL*8. C C Algorithm 358 C Singular Value Decomposition of a Complex Matrix C Peter A. Businger and Gene H. Golub C Communications of the ACM, Volume 12, Number 10, C October, 1969, 564-565. C C SVD finds the singular values S(1) >= S(2) >= ... >= S(N) of an C M by N matrix A with M >= N. The computed singular values are C stored in the vector S. (B, C, T are work vectors of dimension N C or more). SVD also finds the first NU columns of an M by N C orthogonal matrix U and the first NV columns of an N by N C orthogonal matrix V such that ||A - U*D*V'|| is negligible C relative to ||A||, where D = diag[S(1), S(2), ..., S(N)]. (The C only values permitted for NU are 0, N, or M; those for NV are 0 C or N.) Moreover, the transformation T(U) is applied to the P C vectors given in columns N_1, N+2, N+P of the array A. This C feature can be used as follows to find the least squares solution C of minimal Euclidean length (the pseudoinverse solution) of an C over determined system Ax ~= b. Call SVD with NV = N and with C columns N+1, N+2, ..., N+P of A containing P right hand sides b. C From the computed singular values determine the rank r of C diag(S) and define R = diag[1/S(1), ..., 1/S(r), 0, ..., 0]. C Now x = VRd, where d = U'b is furnished by SVD in place of each C right-hand side b. C C SVD can also be used to solve a homogeneous system of linear C equations. To find an orthonormal basis for all solutions of C the system Ax = 0, call SVD with NV = N. The desired basis C consists of those columns of V which correspond to negligible C singular values. Further applications are mentioned in the C references. C C The constants used in the program for ETA and TOL are machine C dependent. ETA is the relative machine precision, TOL the C smallest normalized positive number divided by ETA. The C assignments made are valid for a GE635 computer (a two's C complement binary machine with a signed 27-bit mantissa and a C signed 7-bit exponent). For this machine, ETA = 2**(-26) = 1.5E-8 C and TOL = 2**(-129)/2**(-26) = 1.E-31. C C References C C Golub, G. Least squares singular values and matrix approximations. C Aplickace Matematiky 13 (1968), 44-51. C C Golub, G., and Kahan, W. Calculating the singular value and C pseudoinverse of a matrix. J. SIAM Numer. Anal. (1965), 205-224. C C Golub, G., and Reinsch, C. Singular value decomposition and least C squares solutions. Numer. Math. (to appear). SUBROUTINE DSVD(A,MMAX,NMAX,M,N,P,NU,NV,S,U,V,B,C,T) C IMPLICIT REAL*8 (A-H,O-Z) save REAL*8 A(MMAX,1), U(MMAX,1), V(NMAX,1), S(1), B(1), C(1), T(1) INTEGER M,N,P,NU,NV DATA ETA,TOL/0.222D-15,0.4D-289/ NP=N+P N1=N+1 C C HOUSEHOLDER REDUCTION C(1)=0.D0 K=1 10 K1=K+1 C C ELIMINATION OF A(I,K), I=K+1....,M Z=0.D0 DO 20 I=K,M 20 Z=Z+A(I,K)**2 B(K)=0.D0 IF(Z.LE.TOL)GO TO 70 Z=DSQRT(Z) B(K)=Z W=DABS(A(K,K)) Q=1.D0 IF(W.NE.0.D0)Q=A(K,K)/W A(K,K)=Q*(Z+W) IF(K.EQ.NP)GO TO 70 DO 50 J=K1,NP Q=0.D0 DO 30 I=K,M 30 Q=Q+A(I,K)*A(I,J) Q=Q/(Z*(Z+W)) DO 40 I=K,M 40 A(I,J)=A(I,J)-Q*A(I,K) 50 CONTINUE C C PHASE TRANSFORMATION Q=-A(K,K)/DABS(A(K,K)) DO 60 J=K1,NP 60 A(K,J)=Q*A(K,J) C C ELIMINATION OF A(K,J), J=K+Z,...,N 70 IF(K.EQ.N)GO TO 140 Z=0.D0 DO 80 J=K1,N 80 Z=Z+A(K,J)**2 C(K1)=0.D0 IF(Z.LE.TOL)GO TO 130 Z=DSQRT(Z) C(K1)=Z W=DABS(A(K,K1)) Q=1.D0 IF(W.NE.0.D0)Q=A(K,K1)/W A(K,K1)=Q*(Z+W) DO 110 I=K1,M Q=0.D0 DO 90 J=K1,N 90 Q=Q+A(K,J)*A(I,J) Q=Q/(Z*(Z+W)) DO 100 J=K1,N 100 A(I,J)=A(I,J)-Q*A(K,J) 110 CONTINUE C C PHASE TRANSFORMATION Q=-A(K,K1)/DABS(A(K,K1)) DO 120 I=K1,M 120 A(I,K1)=A(I,K1)*Q 130 K=K1 GO TO 10 C C TOLERANCE FOR NEGLIGIBLE ELEMENTS 140 EPS=0.D0 DO 150 K=1,N S(K)=B(K) T(K)=C(K) 150 EPS=DMAX1(EPS,S(K)+T(K)) EPS=EPS*ETA C C INITIALIZATION OF U AND V IF(NU.EQ.0)GO TO 180 DO 170 J=1,NU DO 160 I=1,M 160 U(I,J)=0.D0 170 U(J,J)=1.D0 180 IF(NV.EQ.0)GO TO 210 DO 200 J=1,NV DO 190 I=1,N 190 V(I,J)=0.D0 200 V(J,J)=1.D0 C C QR DIAGONALIZATION 210 DO 380 KK=1,N K=N1-KK C C TEST FOR SPLIT 220 DO 230 LL=1,K L=K+1-LL IF(DABS(T(L)).LE.EPS)GO TO 290 IF(DABS(S(L-1)).LE.EPS)GO TO 240 230 CONTINUE C C CANCELLATION OF E(L) 240 CS=0.D0 SN=1.D0 L1=L-1 DO 280 I=L,K F=SN*T(I) T(I)=CS*T(I) IF(DABS(F).LE.EPS)GO TO 290 H=S(I) W=DSQRT(F*F+H*H) S(I)=W CS=H/W SN=-F/W IF(NU.EQ.0)GO TO 260 DO 250 J=1,N X=U(J,L1) Y=U(J,I) U(J,L1)=X*CS+Y*SN 250 U(J,I)=Y*CS-X*SN 260 IF(NP.EQ.N)GO TO 280 DO 270 J=N1,NP Q=A(L1,J) R=A(I,J) A(L1,J)=Q*CS+R*SN 270 A(I,J)=R*CS-Q*SN 280 CONTINUE C C TEST FOR CONVERGENCE 290 W=S(K) IF(L.EQ.K)GO TO 360 C C ORIGIN SHIFT X=S(L) Y=S(K-1) G=T(K-1) H=T(K) F=((Y-W)*(Y+W)+(G-H)*(G+H))/(2.D0*H*Y) G=DSQRT(F*F+1.D0) IF(F.LT.0.D0)G=-G F=((X-W)*(X+W)+(Y/(F+G)-H)*H)/X C C QR STEP CS=1.D0 SN=1.D0 L1=L+1 DO 350 I=L1,K G=T(I) Y=S(I) H=SN*G G=CS*G W=DSQRT(H*H+F*F) T(I-1)=W CS=F/W SN=H/W F=X*CS+G*SN G=G*CS-X*SN H=Y*SN Y=Y*CS IF(NV.EQ.0)GO TO 310 DO 300 J=1,N X=V(J,I-1) W=V(J,I) V(J,I-1)=X*CS+W*SN 300 V(J,I)=W*CS-X*SN 310 W=DSQRT(H*H+F*F) S(I-1)=W CS=F/W SN=H/W F=CS*G+SN*Y X=CS*Y-SN*G IF(NU.EQ.0)GO TO 330 DO 320 J=1,N Y=U(J,I-1) W=U(J,I) U(J,I-1)=Y*CS+W*SN 320 U(J,I)=W*CS-Y*SN 330 IF(N.EQ.NP)GO TO 350 DO 340 J=N1,NP Q=A(I-1,J) R=A(I,J) A(I-1,J)=Q*CS+R*SN 340 A(I,J)=R*CS-Q*SN 350 CONTINUE T(L)=0.D0 T(K)=F S(K)=X GO TO 220 C C CONVERGENCE 360 IF(W.GE.0.D0)GO TO 380 S(K)=-W IF(NV.EQ.0)GO TO 380 DO 370 J=1,N 370 V(J,K)=-V(J,K) 380 CONTINUE C C SORT SINGULAR VALUES DO 450 K=1,N G=-1.D0 J=K DO 390 I=K,N IF(S(I).LE.G)GO TO 390 G=S(I) J=I 390 CONTINUE IF(J.EQ.K)GO TO 450 S(J)=S(K) S(K)=G IF(NV.EQ.0)GO TO 410 DO 400 I=1,N Q=V(I,J) V(I,J)=V(I,K) 400 V(I,K)=Q 410 IF(NU.EQ.0)GO TO 430 DO 420 I=1,N Q=U(I,J) U(I,J)=U(I,K) 420 U(I,K)=Q 430 IF(N.EQ.NP)GO TO 450 DO 440 I=N1,NP Q=A(J,I) A(J,I)=A(K,I) 440 A(K,I)=Q 450 CONTINUE C C BACK TRANSFORMATION IF(NU.EQ.0)GO TO 510 DO 500 KK=1,N K=N1-KK IF(B(K).EQ.0.D0)GO TO 500 Q=-A(K,K)/DABS(A(K,K)) DO 460 J=1,NU 460 U(K,J)=Q*U(K,J) DO 490 J=1,NU Q=0.D0 DO 470 I=K,M 470 Q=Q+A(I,K)*U(I,J) Q=Q/(DABS(A(K,K))*B(K)) DO 480 I=K,M 480 U(I,J)=U(I,J)-Q*A(I,K) 490 CONTINUE 500 CONTINUE 510 IF(NV.EQ.0)GO TO 570 IF(N.LT.2)GO TO 570 DO 560 KK=2,N K=N1-KK K1=K+1 IF(C(K1).EQ.0.D0)GO TO 560 Q=-A(K,K1)/DABS(A(K,K1)) DO 520 J=1,NV 520 V(K1,J)=Q*V(K1,J) DO 550 J=1,NV Q=0.D0 DO 530 I=K1,N 530 Q=Q+A(K,I)*V(I,J) Q=Q/(DABS(A(K,K1))*C(K1)) DO 540 I=K1,N 540 V(I,J)=V(I,J)-Q*A(K,I) 550 CONTINUE 560 CONTINUE 570 RETURN END