1      SUBROUTINE CSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)
2C***BEGIN PROLOGUE  CSVDC
3C***DATE WRITTEN   790319   (YYMMDD)
4C***REVISION DATE  820801   (YYMMDD)
5C***REVISION HISTORY  (YYMMDD)
6C   000330  Modified array declarations.  (JEC)
7C***CATEGORY NO.  D6
8C***KEYWORDS  COMPLEX,LINEAR ALGEBRA,LINPACK,MATRIX,
9C             SINGULAR VALUE DECOMPOSITION
10C***AUTHOR  STEWART, G. W., (U. OF MARYLAND)
11C***PURPOSE  Perform the singular value decomposition of a COMPLEX NXP
12C            matrix.
13C***DESCRIPTION
14C
15C     CSVDC is a subroutine to reduce a complex NxP matrix X by
16C     unitary transformations U and V to diagonal form.  The
17C     diagonal elements S(I) are the singular values of X.  The
18C     columns of U are the corresponding left singular vectors,
19C     and the columns of V the right singular vectors.
20C
21C     On Entry
22C
23C         X         COMPLEX(LDX,P), where LDX .GE. N.
24C                   X contains the matrix whose singular value
25C                   decomposition is to be computed.  X is
26C                   destroyed by CSVDC.
27C
28C         LDX       INTEGER.
29C                   LDX is the leading dimension of the array X.
30C
31C         N         INTEGER.
32C                   N is the number of rows of the matrix X.
33C
34C         P         INTEGER.
35C                   P is the number of columns of the matrix X.
36C
37C         LDU       INTEGER.
38C                   LDU is the leading dimension of the array U
39C                   (see below).
40C
41C         LDV       INTEGER.
42C                   LDV is the leading dimension of the array V
43C                   (see below).
44C
45C         WORK      COMPLEX(N).
46C                   WORK is a scratch array.
47C
48C         JOB       INTEGER.
49C                   JOB controls the computation of the singular
50C                   vectors.  It has the decimal expansion AB
51C                   with the following meaning
52C
53C                        A .EQ. 0    Do not compute the left singular
54C                                    vectors.
55C                        A .EQ. 1    Return the N left singular vectors
56C                                    in U.
57C                        A .GE. 2    Return the first MIN(N,P)
58C                                    left singular vectors in U.
59C                        B .EQ. 0    Do not compute the right singular
60C                                    vectors.
61C                        B .EQ. 1    Return the right singular vectors
62C                                    in V.
63C
64C     On Return
65C
66C         S         COMPLEX(MM), where MM = MIN(N+1,P).
67C                   The first MIN(N,P) entries of S contain the
68C                   singular values of X arranged in descending
69C                   order of magnitude.
70C
71C         E         COMPLEX(P).
72C                   E ordinarily contains zeros.  However see the
73C                   discussion of INFO for exceptions.
74C
75C         U         COMPLEX(LDU,K), where LDU .GE. N.  If JOBA .EQ. 1
76C                                   then K .EQ. N.  If JOBA .GE. 2 then
77C                                   K .EQ. MIN(N,P).
78C                   U contains the matrix of right singular vectors.
79C                   U is not referenced if JOBA .EQ. 0.  If N .LE. P
80C                   or if JOBA .GT. 2, then U may be identified with X
81C                   in the subroutine call.
82C
83C         V         COMPLEX(LDV,P), where LDV .GE. P.
84C                   V contains the matrix of right singular vectors.
85C                   V is not referenced if JOB .EQ. 0.  If P .LE. N,
86C                   then V may be identified with X in the
87C                   subroutine call.
88C
89C         INFO      INTEGER.
90C                   The singular values (and their corresponding
91C                   singular vectors) S(INFO+1),S(INFO+2),...,S(M)
92C                   are correct (here M=MIN(N,P)).  Thus if
93C                   INFO.EQ. 0, all the singular values and their
94C                   vectors are correct.  In any event, the matrix
95C                   B = CTRANS(U)*X*V is the bidiagonal matrix
96C                   with the elements of S on its diagonal and the
97C                   elements of E on its super-diagonal (CTRANS(U)
98C                   is the conjugate-transpose of U).  Thus the
99C                   singular values of X and B are the same.
100C
101C     LINPACK.  This version dated 03/19/79 .
102C     Stewart, G. W., University of Maryland, Argonne National Lab.
103C
104C     CSVDC uses the following functions and subprograms.
105C
106C     External CSROT
107C     BLAS CAXPY,CDOTC,CSCAL,CSWAP,SCNRM2,SROTG
108C     Fortran ABS,AIMAG,AMAX1,CABS,CMPLX
109C     Fortran CONJG,MAX0,MIN0,MOD,REAL,SQRT
110C***REFERENCES  DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,
111C                 *LINPACK USERS  GUIDE*, SIAM, 1979.
112C***ROUTINES CALLED  CAXPY,CDOTC,CSCAL,CSROT,CSWAP,SCNRM2,SROTG
113C***END PROLOGUE  CSVDC
114      INTEGER LDX,N,P,LDU,LDV,JOB,INFO
115      COMPLEX X(LDX,*),S(*),E(*),U(LDU,*),V(LDV,*),WORK(*)
116C
117C
118      INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
119     1        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
120      COMPLEX CDOTC,T,R
121      REAL B,C,CS,EL,EMM1,F,G,SCNRM2,SCALE,SHIFT,SL,SM,SN,SMM1,T1,TEST,
122     1     ZTEST
123      LOGICAL WANTU,WANTV
124      COMPLEX CSIGN,ZDUM,ZDUM1,ZDUM2
125      REAL CABS1
126      CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM))
127      CSIGN(ZDUM1,ZDUM2) = CABS(ZDUM1)*(ZDUM2/CABS(ZDUM2))
128C
129C     SET THE MAXIMUM NUMBER OF ITERATIONS.
130C
131C***FIRST EXECUTABLE STATEMENT  CSVDC
132      MAXIT = 30
133C
134C     DETERMINE WHAT IS TO BE COMPUTED.
135C
136      WANTU = .FALSE.
137      WANTV = .FALSE.
138      JOBU = MOD(JOB,100)/10
139      NCU = N
140      IF (JOBU .GT. 1) NCU = MIN0(N,P)
141      IF (JOBU .NE. 0) WANTU = .TRUE.
142      IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
143C
144C     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
145C     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
146C
147      INFO = 0
148      NCT = MIN0(N-1,P)
149      NRT = MAX0(0,MIN0(P-2,N))
150      LU = MAX0(NCT,NRT)
151      IF (LU .LT. 1) GO TO 170
152      DO 160 L = 1, LU
153         LP1 = L + 1
154         IF (L .GT. NCT) GO TO 20
155C
156C           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
157C           PLACE THE L-TH DIAGONAL IN S(L).
158C
159            S(L) = CMPLX(SCNRM2(N-L+1,X(L,L),1),0.0E0)
160            IF (CABS1(S(L)) .EQ. 0.0E0) GO TO 10
161               IF (CABS1(X(L,L)) .NE. 0.0E0) S(L) = CSIGN(S(L),X(L,L))
162               CALL CSCAL(N-L+1,1.0E0/S(L),X(L,L),1)
163               X(L,L) = (1.0E0,0.0E0) + X(L,L)
164   10       CONTINUE
165            S(L) = -S(L)
166   20    CONTINUE
167         IF (P .LT. LP1) GO TO 50
168         DO 40 J = LP1, P
169            IF (L .GT. NCT) GO TO 30
170            IF (CABS1(S(L)) .EQ. 0.0E0) GO TO 30
171C
172C              APPLY THE TRANSFORMATION.
173C
174               T = -CDOTC(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
175               CALL CAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
176   30       CONTINUE
177C
178C           PLACE THE L-TH ROW OF X INTO  E FOR THE
179C           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
180C
181            E(J) = CONJG(X(L,J))
182   40    CONTINUE
183   50    CONTINUE
184         IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70
185C
186C           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
187C           MULTIPLICATION.
188C
189            DO 60 I = L, N
190               U(I,L) = X(I,L)
191   60       CONTINUE
192   70    CONTINUE
193         IF (L .GT. NRT) GO TO 150
194C
195C           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
196C           L-TH SUPER-DIAGONAL IN E(L).
197C
198            E(L) = CMPLX(SCNRM2(P-L,E(LP1),1),0.0E0)
199            IF (CABS1(E(L)) .EQ. 0.0E0) GO TO 80
200               IF (CABS1(E(LP1)) .NE. 0.0E0) E(L) = CSIGN(E(L),E(LP1))
201               CALL CSCAL(P-L,1.0E0/E(L),E(LP1),1)
202               E(LP1) = (1.0E0,0.0E0) + E(LP1)
203   80       CONTINUE
204            E(L) = -CONJG(E(L))
205            IF (LP1 .GT. N .OR. CABS1(E(L)) .EQ. 0.0E0) GO TO 120
206C
207C              APPLY THE TRANSFORMATION.
208C
209               DO 90 I = LP1, N
210                  WORK(I) = (0.0E0,0.0E0)
211   90          CONTINUE
212               DO 100 J = LP1, P
213                  CALL CAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1)
214  100          CONTINUE
215               DO 110 J = LP1, P
216                  CALL CAXPY(N-L,CONJG(-E(J)/E(LP1)),WORK(LP1),1,
217     1                       X(LP1,J),1)
218  110          CONTINUE
219  120       CONTINUE
220            IF (.NOT.WANTV) GO TO 140
221C
222C              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
223C              BACK MULTIPLICATION.
224C
225               DO 130 I = LP1, P
226                  V(I,L) = E(I)
227  130          CONTINUE
228  140       CONTINUE
229  150    CONTINUE
230  160 CONTINUE
231  170 CONTINUE
232C
233C     SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
234C
235      M = MIN0(P,N+1)
236      NCTP1 = NCT + 1
237      NRTP1 = NRT + 1
238      IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)
239      IF (N .LT. M) S(M) = (0.0E0,0.0E0)
240      IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)
241      E(M) = (0.0E0,0.0E0)
242C
243C     IF REQUIRED, GENERATE U.
244C
245      IF (.NOT.WANTU) GO TO 300
246         IF (NCU .LT. NCTP1) GO TO 200
247         DO 190 J = NCTP1, NCU
248            DO 180 I = 1, N
249               U(I,J) = (0.0E0,0.0E0)
250  180       CONTINUE
251            U(J,J) = (1.0E0,0.0E0)
252  190    CONTINUE
253  200    CONTINUE
254         IF (NCT .LT. 1) GO TO 290
255         DO 280 LL = 1, NCT
256            L = NCT - LL + 1
257            IF (CABS1(S(L)) .EQ. 0.0E0) GO TO 250
258               LP1 = L + 1
259               IF (NCU .LT. LP1) GO TO 220
260               DO 210 J = LP1, NCU
261                  T = -CDOTC(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)
262                  CALL CAXPY(N-L+1,T,U(L,L),1,U(L,J),1)
263  210          CONTINUE
264  220          CONTINUE
265               CALL CSCAL(N-L+1,(-1.0E0,0.0E0),U(L,L),1)
266               U(L,L) = (1.0E0,0.0E0) + U(L,L)
267               LM1 = L - 1
268               IF (LM1 .LT. 1) GO TO 240
269               DO 230 I = 1, LM1
270                  U(I,L) = (0.0E0,0.0E0)
271  230          CONTINUE
272  240          CONTINUE
273            GO TO 270
274  250       CONTINUE
275               DO 260 I = 1, N
276                  U(I,L) = (0.0E0,0.0E0)
277  260          CONTINUE
278               U(L,L) = (1.0E0,0.0E0)
279  270       CONTINUE
280  280    CONTINUE
281  290    CONTINUE
282  300 CONTINUE
283C
284C     IF IT IS REQUIRED, GENERATE V.
285C
286      IF (.NOT.WANTV) GO TO 350
287         DO 340 LL = 1, P
288            L = P - LL + 1
289            LP1 = L + 1
290            IF (L .GT. NRT) GO TO 320
291            IF (CABS1(E(L)) .EQ. 0.0E0) GO TO 320
292               DO 310 J = LP1, P
293                  T = -CDOTC(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)
294                  CALL CAXPY(P-L,T,V(LP1,L),1,V(LP1,J),1)
295  310          CONTINUE
296  320       CONTINUE
297            DO 330 I = 1, P
298               V(I,L) = (0.0E0,0.0E0)
299  330       CONTINUE
300            V(L,L) = (1.0E0,0.0E0)
301  340    CONTINUE
302  350 CONTINUE
303C
304C     TRANSFORM S AND E SO THAT THEY ARE REAL.
305C
306      DO 380 I = 1, M
307         IF (CABS1(S(I)) .EQ. 0.0E0) GO TO 360
308            T = CMPLX(CABS(S(I)),0.0E0)
309            R = S(I)/T
310            S(I) = T
311            IF (I .LT. M) E(I) = E(I)/R
312            IF (WANTU) CALL CSCAL(N,R,U(1,I),1)
313  360    CONTINUE
314C     ...EXIT
315         IF (I .EQ. M) GO TO 390
316         IF (CABS1(E(I)) .EQ. 0.0E0) GO TO 370
317            T = CMPLX(CABS(E(I)),0.0E0)
318            R = T/E(I)
319            E(I) = T
320            S(I+1) = S(I+1)*R
321            IF (WANTV) CALL CSCAL(P,R,V(1,I+1),1)
322  370    CONTINUE
323  380 CONTINUE
324  390 CONTINUE
325C
326C     MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
327C
328      MM = M
329      ITER = 0
330  400 CONTINUE
331C
332C        QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
333C
334C     ...EXIT
335         IF (M .EQ. 0) GO TO 660
336C
337C        IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
338C        FLAG AND RETURN.
339C
340         IF (ITER .LT. MAXIT) GO TO 410
341            INFO = M
342C     ......EXIT
343            GO TO 660
344  410    CONTINUE
345C
346C        THIS SECTION OF THE PROGRAM INSPECTS FOR
347C        NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS.  ON
348C        COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS.
349C
350C           KASE = 1     IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M
351C           KASE = 2     IF S(L) IS NEGLIGIBLE AND L.LT.M
352C           KASE = 3     IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND
353C                        S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
354C           KASE = 4     IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).
355C
356         DO 430 LL = 1, M
357            L = M - LL
358C        ...EXIT
359            IF (L .EQ. 0) GO TO 440
360            TEST = CABS(S(L)) + CABS(S(L+1))
361            ZTEST = TEST + CABS(E(L))
362            IF (ZTEST .NE. TEST) GO TO 420
363               E(L) = (0.0E0,0.0E0)
364C        ......EXIT
365               GO TO 440
366  420       CONTINUE
367  430    CONTINUE
368  440    CONTINUE
369         IF (L .NE. M - 1) GO TO 450
370            KASE = 4
371         GO TO 520
372  450    CONTINUE
373            LP1 = L + 1
374            MP1 = M + 1
375            DO 470 LLS = LP1, MP1
376               LS = M - LLS + LP1
377C           ...EXIT
378               IF (LS .EQ. L) GO TO 480
379               TEST = 0.0E0
380               IF (LS .NE. M) TEST = TEST + CABS(E(LS))
381               IF (LS .NE. L + 1) TEST = TEST + CABS(E(LS-1))
382               ZTEST = TEST + CABS(S(LS))
383               IF (ZTEST .NE. TEST) GO TO 460
384                  S(LS) = (0.0E0,0.0E0)
385C           ......EXIT
386                  GO TO 480
387  460          CONTINUE
388  470       CONTINUE
389  480       CONTINUE
390            IF (LS .NE. L) GO TO 490
391               KASE = 3
392            GO TO 510
393  490       CONTINUE
394            IF (LS .NE. M) GO TO 500
395               KASE = 1
396            GO TO 510
397  500       CONTINUE
398               KASE = 2
399               L = LS
400  510       CONTINUE
401  520    CONTINUE
402         L = L + 1
403C
404C        PERFORM THE TASK INDICATED BY KASE.
405C
406         GO TO (530, 560, 580, 610), KASE
407C
408C        DEFLATE NEGLIGIBLE S(M).
409C
410  530    CONTINUE
411            MM1 = M - 1
412            F = REAL(E(M-1))
413            E(M-1) = (0.0E0,0.0E0)
414            DO 550 KK = L, MM1
415               K = MM1 - KK + L
416               T1 = REAL(S(K))
417               CALL SROTG(T1,F,CS,SN)
418               S(K) = CMPLX(T1,0.0E0)
419               IF (K .EQ. L) GO TO 540
420                  F = -SN*REAL(E(K-1))
421                  E(K-1) = CS*E(K-1)
422  540          CONTINUE
423               IF (WANTV) CALL CSROT(P,V(1,K),1,V(1,M),1,CS,SN)
424  550       CONTINUE
425         GO TO 650
426C
427C        SPLIT AT NEGLIGIBLE S(L).
428C
429  560    CONTINUE
430            F = REAL(E(L-1))
431            E(L-1) = (0.0E0,0.0E0)
432            DO 570 K = L, M
433               T1 = REAL(S(K))
434               CALL SROTG(T1,F,CS,SN)
435               S(K) = CMPLX(T1,0.0E0)
436               F = -SN*REAL(E(K))
437               E(K) = CS*E(K)
438               IF (WANTU) CALL CSROT(N,U(1,K),1,U(1,L-1),1,CS,SN)
439  570       CONTINUE
440         GO TO 650
441C
442C        PERFORM ONE QR STEP.
443C
444  580    CONTINUE
445C
446C           CALCULATE THE SHIFT.
447C
448            SCALE = AMAX1(CABS(S(M)),CABS(S(M-1)),CABS(E(M-1)),
449     1                    CABS(S(L)),CABS(E(L)))
450            SM = REAL(S(M))/SCALE
451            SMM1 = REAL(S(M-1))/SCALE
452            EMM1 = REAL(E(M-1))/SCALE
453            SL = REAL(S(L))/SCALE
454            EL = REAL(E(L))/SCALE
455            B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0E0
456            C = (SM*EMM1)**2
457            SHIFT = 0.0E0
458            IF (B .EQ. 0.0E0 .AND. C .EQ. 0.0E0) GO TO 590
459               SHIFT = SQRT(B**2+C)
460               IF (B .LT. 0.0E0) SHIFT = -SHIFT
461               SHIFT = C/(B + SHIFT)
462  590       CONTINUE
463            F = (SL + SM)*(SL - SM) - SHIFT
464            G = SL*EL
465C
466C           CHASE ZEROS.
467C
468            MM1 = M - 1
469            DO 600 K = L, MM1
470               CALL SROTG(F,G,CS,SN)
471               IF (K .NE. L) E(K-1) = CMPLX(F,0.0E0)
472               F = CS*REAL(S(K)) + SN*REAL(E(K))
473               E(K) = CS*E(K) - SN*S(K)
474               G = SN*REAL(S(K+1))
475               S(K+1) = CS*S(K+1)
476               IF (WANTV) CALL CSROT(P,V(1,K),1,V(1,K+1),1,CS,SN)
477               CALL SROTG(F,G,CS,SN)
478               S(K) = CMPLX(F,0.0E0)
479               F = CS*REAL(E(K)) + SN*REAL(S(K+1))
480               S(K+1) = -SN*E(K) + CS*S(K+1)
481               G = SN*REAL(E(K+1))
482               E(K+1) = CS*E(K+1)
483               IF (WANTU .AND. K .LT. N)
484     1            CALL CSROT(N,U(1,K),1,U(1,K+1),1,CS,SN)
485  600       CONTINUE
486            E(M-1) = CMPLX(F,0.0E0)
487            ITER = ITER + 1
488         GO TO 650
489C
490C        CONVERGENCE.
491C
492  610    CONTINUE
493C
494C           MAKE THE SINGULAR VALUE  POSITIVE
495C
496            IF (REAL(S(L)) .GE. 0.0E0) GO TO 620
497               S(L) = -S(L)
498               IF (WANTV) CALL CSCAL(P,(-1.0E0,0.0E0),V(1,L),1)
499  620       CONTINUE
500C
501C           ORDER THE SINGULAR VALUE.
502C
503  630       IF (L .EQ. MM) GO TO 640
504C           ...EXIT
505               IF (REAL(S(L)) .GE. REAL(S(L+1))) GO TO 640
506               T = S(L)
507               S(L) = S(L+1)
508               S(L+1) = T
509               IF (WANTV .AND. L .LT. P)
510     1            CALL CSWAP(P,V(1,L),1,V(1,L+1),1)
511               IF (WANTU .AND. L .LT. N)
512     1            CALL CSWAP(N,U(1,L),1,U(1,L+1),1)
513               L = L + 1
514            GO TO 630
515  640       CONTINUE
516            ITER = 0
517            M = M - 1
518  650    CONTINUE
519      GO TO 400
520  660 CONTINUE
521      RETURN
522      END
523