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
