1*DECK SSVDC
2      SUBROUTINE SSVDC (X, LDX, N, P, S, E, U, LDU, V, LDV, WORK, JOB,
3     +   INFO)
4C***BEGIN PROLOGUE  SSVDC
5C***PURPOSE  Perform the singular value decomposition of a rectangular
6C            matrix.
7C***LIBRARY   SLATEC (LINPACK)
8C***CATEGORY  D6
9C***TYPE      SINGLE PRECISION (SSVDC-S, DSVDC-D, CSVDC-C)
10C***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX,
11C             SINGULAR VALUE DECOMPOSITION
12C***AUTHOR  Stewart, G. W., (U. of Maryland)
13C***DESCRIPTION
14C
15C     SSVDC is a subroutine to reduce a real NxP matrix X by orthogonal
16C     transformations U and V to diagonal form.  The elements S(I) are
17C     the singular values of X.  The columns of U are the corresponding
18C     left singular vectors, and the columns of V the right singular
19C     vectors.
20C
21C     On Entry
22C
23C         X         REAL(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 SSVDC.
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      REAL(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) singular
58C                                  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         REAL(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         REAL(P).
72C                   E ordinarily contains zeros.  However, see the
73C                   discussion of INFO for exceptions.
74C
75C         U         REAL(LDU,K), where LDU .GE. N.  If JOBA .EQ. 1, then
76C                                   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 .EQ. 2, then U may be identified with X
81C                   in the subroutine call.
82C
83C         V         REAL(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 = TRANS(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 (TRANS(U)
98C                   is the transpose of U).  Thus the singular
99C                   values of X and B are the same.
100C
101C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
102C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
103C***ROUTINES CALLED  SAXPY, SDOT, SNRM2, SROT, SROTG, SSCAL, SSWAP
104C***REVISION HISTORY  (YYMMDD)
105C   790319  DATE WRITTEN
106C   890531  Changed all specific intrinsics to generic.  (WRB)
107C   890531  REVISION DATE from Version 3.2
108C   891214  Prologue converted to Version 4.0 format.  (BAB)
109C   900326  Removed duplicate information from DESCRIPTION section.
110C           (WRB)
111C   920501  Reformatted the REFERENCES section.  (WRB)
112C***END PROLOGUE  SSVDC
113      INTEGER LDX,N,P,LDU,LDV,JOB,INFO
114      REAL X(LDX,*),S(*),E(*),U(LDU,*),V(LDV,*),WORK(*)
115C
116C
117      INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
118     1        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
119      REAL SDOT,T
120      REAL B,C,CS,EL,EMM1,F,G,SNRM2,SCALE,SHIFT,SL,SM,SN,SMM1,T1,TEST,
121     1     ZTEST
122      LOGICAL WANTU,WANTV
123C***FIRST EXECUTABLE STATEMENT  SSVDC
124C
125C     SET THE MAXIMUM NUMBER OF ITERATIONS.
126C
127      MAXIT = 30
128C
129C     DETERMINE WHAT IS TO BE COMPUTED.
130C
131      WANTU = .FALSE.
132      WANTV = .FALSE.
133      JOBU = MOD(JOB,100)/10
134      NCU = N
135      IF (JOBU .GT. 1) NCU = MIN(N,P)
136      IF (JOBU .NE. 0) WANTU = .TRUE.
137      IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
138C
139C     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
140C     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
141C
142      INFO = 0
143      NCT = MIN(N-1,P)
144      NRT = MAX(0,MIN(P-2,N))
145      LU = MAX(NCT,NRT)
146      IF (LU .LT. 1) GO TO 170
147      DO 160 L = 1, LU
148         LP1 = L + 1
149         IF (L .GT. NCT) GO TO 20
150C
151C           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
152C           PLACE THE L-TH DIAGONAL IN S(L).
153C
154            S(L) = SNRM2(N-L+1,X(L,L),1)
155            IF (S(L) .EQ. 0.0E0) GO TO 10
156               IF (X(L,L) .NE. 0.0E0) S(L) = SIGN(S(L),X(L,L))
157               CALL SSCAL(N-L+1,1.0E0/S(L),X(L,L),1)
158               X(L,L) = 1.0E0 + X(L,L)
159   10       CONTINUE
160            S(L) = -S(L)
161   20    CONTINUE
162         IF (P .LT. LP1) GO TO 50
163         DO 40 J = LP1, P
164            IF (L .GT. NCT) GO TO 30
165            IF (S(L) .EQ. 0.0E0) GO TO 30
166C
167C              APPLY THE TRANSFORMATION.
168C
169               T = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
170               CALL SAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
171   30       CONTINUE
172C
173C           PLACE THE L-TH ROW OF X INTO  E FOR THE
174C           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
175C
176            E(J) = X(L,J)
177   40    CONTINUE
178   50    CONTINUE
179         IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70
180C
181C           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
182C           MULTIPLICATION.
183C
184            DO 60 I = L, N
185               U(I,L) = X(I,L)
186   60       CONTINUE
187   70    CONTINUE
188         IF (L .GT. NRT) GO TO 150
189C
190C           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
191C           L-TH SUPER-DIAGONAL IN E(L).
192C
193            E(L) = SNRM2(P-L,E(LP1),1)
194            IF (E(L) .EQ. 0.0E0) GO TO 80
195               IF (E(LP1) .NE. 0.0E0) E(L) = SIGN(E(L),E(LP1))
196               CALL SSCAL(P-L,1.0E0/E(L),E(LP1),1)
197               E(LP1) = 1.0E0 + E(LP1)
198   80       CONTINUE
199            E(L) = -E(L)
200            IF (LP1 .GT. N .OR. E(L) .EQ. 0.0E0) GO TO 120
201C
202C              APPLY THE TRANSFORMATION.
203C
204               DO 90 I = LP1, N
205                  WORK(I) = 0.0E0
206   90          CONTINUE
207               DO 100 J = LP1, P
208                  CALL SAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1)
209  100          CONTINUE
210               DO 110 J = LP1, P
211                  CALL SAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1)
212  110          CONTINUE
213  120       CONTINUE
214            IF (.NOT.WANTV) GO TO 140
215C
216C              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
217C              BACK MULTIPLICATION.
218C
219               DO 130 I = LP1, P
220                  V(I,L) = E(I)
221  130          CONTINUE
222  140       CONTINUE
223  150    CONTINUE
224  160 CONTINUE
225  170 CONTINUE
226C
227C     SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
228C
229      M = MIN(P,N+1)
230      NCTP1 = NCT + 1
231      NRTP1 = NRT + 1
232      IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)
233      IF (N .LT. M) S(M) = 0.0E0
234      IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)
235      E(M) = 0.0E0
236C
237C     IF REQUIRED, GENERATE U.
238C
239      IF (.NOT.WANTU) GO TO 300
240         IF (NCU .LT. NCTP1) GO TO 200
241         DO 190 J = NCTP1, NCU
242            DO 180 I = 1, N
243               U(I,J) = 0.0E0
244  180       CONTINUE
245            U(J,J) = 1.0E0
246  190    CONTINUE
247  200    CONTINUE
248         IF (NCT .LT. 1) GO TO 290
249         DO 280 LL = 1, NCT
250            L = NCT - LL + 1
251            IF (S(L) .EQ. 0.0E0) GO TO 250
252               LP1 = L + 1
253               IF (NCU .LT. LP1) GO TO 220
254               DO 210 J = LP1, NCU
255                  T = -SDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)
256                  CALL SAXPY(N-L+1,T,U(L,L),1,U(L,J),1)
257  210          CONTINUE
258  220          CONTINUE
259               CALL SSCAL(N-L+1,-1.0E0,U(L,L),1)
260               U(L,L) = 1.0E0 + U(L,L)
261               LM1 = L - 1
262               IF (LM1 .LT. 1) GO TO 240
263               DO 230 I = 1, LM1
264                  U(I,L) = 0.0E0
265  230          CONTINUE
266  240          CONTINUE
267            GO TO 270
268  250       CONTINUE
269               DO 260 I = 1, N
270                  U(I,L) = 0.0E0
271  260          CONTINUE
272               U(L,L) = 1.0E0
273  270       CONTINUE
274  280    CONTINUE
275  290    CONTINUE
276  300 CONTINUE
277C
278C     IF IT IS REQUIRED, GENERATE V.
279C
280      IF (.NOT.WANTV) GO TO 350
281         DO 340 LL = 1, P
282            L = P - LL + 1
283            LP1 = L + 1
284            IF (L .GT. NRT) GO TO 320
285            IF (E(L) .EQ. 0.0E0) GO TO 320
286               DO 310 J = LP1, P
287                  T = -SDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)
288                  CALL SAXPY(P-L,T,V(LP1,L),1,V(LP1,J),1)
289  310          CONTINUE
290  320       CONTINUE
291            DO 330 I = 1, P
292               V(I,L) = 0.0E0
293  330       CONTINUE
294            V(L,L) = 1.0E0
295  340    CONTINUE
296  350 CONTINUE
297C
298C     MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
299C
300      MM = M
301      ITER = 0
302  360 CONTINUE
303C
304C        QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
305C
306         IF (M .EQ. 0) GO TO 620
307C
308C        IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
309C        FLAG AND RETURN.
310C
311         IF (ITER .LT. MAXIT) GO TO 370
312            INFO = M
313            GO TO 620
314  370    CONTINUE
315C
316C        THIS SECTION OF THE PROGRAM INSPECTS FOR
317C        NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS.  ON
318C        COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS.
319C
320C           KASE = 1     IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M
321C           KASE = 2     IF S(L) IS NEGLIGIBLE AND L.LT.M
322C           KASE = 3     IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND
323C                        S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
324C           KASE = 4     IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).
325C
326         DO 390 LL = 1, M
327            L = M - LL
328            IF (L .EQ. 0) GO TO 400
329            TEST = ABS(S(L)) + ABS(S(L+1))
330            ZTEST = TEST + ABS(E(L))
331            IF (ZTEST .NE. TEST) GO TO 380
332               E(L) = 0.0E0
333               GO TO 400
334  380       CONTINUE
335  390    CONTINUE
336  400    CONTINUE
337         IF (L .NE. M - 1) GO TO 410
338            KASE = 4
339         GO TO 480
340  410    CONTINUE
341            LP1 = L + 1
342            MP1 = M + 1
343            DO 430 LLS = LP1, MP1
344               LS = M - LLS + LP1
345               IF (LS .EQ. L) GO TO 440
346               TEST = 0.0E0
347               IF (LS .NE. M) TEST = TEST + ABS(E(LS))
348               IF (LS .NE. L + 1) TEST = TEST + ABS(E(LS-1))
349               ZTEST = TEST + ABS(S(LS))
350               IF (ZTEST .NE. TEST) GO TO 420
351                  S(LS) = 0.0E0
352                  GO TO 440
353  420          CONTINUE
354  430       CONTINUE
355  440       CONTINUE
356            IF (LS .NE. L) GO TO 450
357               KASE = 3
358            GO TO 470
359  450       CONTINUE
360            IF (LS .NE. M) GO TO 460
361               KASE = 1
362            GO TO 470
363  460       CONTINUE
364               KASE = 2
365               L = LS
366  470       CONTINUE
367  480    CONTINUE
368         L = L + 1
369C
370C        PERFORM THE TASK INDICATED BY KASE.
371C
372         GO TO (490,520,540,570), KASE
373C
374C        DEFLATE NEGLIGIBLE S(M).
375C
376  490    CONTINUE
377            MM1 = M - 1
378            F = E(M-1)
379            E(M-1) = 0.0E0
380            DO 510 KK = L, MM1
381               K = MM1 - KK + L
382               T1 = S(K)
383               CALL SROTG(T1,F,CS,SN)
384               S(K) = T1
385               IF (K .EQ. L) GO TO 500
386                  F = -SN*E(K-1)
387                  E(K-1) = CS*E(K-1)
388  500          CONTINUE
389               IF (WANTV) CALL SROT(P,V(1,K),1,V(1,M),1,CS,SN)
390  510       CONTINUE
391         GO TO 610
392C
393C        SPLIT AT NEGLIGIBLE S(L).
394C
395  520    CONTINUE
396            F = E(L-1)
397            E(L-1) = 0.0E0
398            DO 530 K = L, M
399               T1 = S(K)
400               CALL SROTG(T1,F,CS,SN)
401               S(K) = T1
402               F = -SN*E(K)
403               E(K) = CS*E(K)
404               IF (WANTU) CALL SROT(N,U(1,K),1,U(1,L-1),1,CS,SN)
405  530       CONTINUE
406         GO TO 610
407C
408C        PERFORM ONE QR STEP.
409C
410  540    CONTINUE
411C
412C           CALCULATE THE SHIFT.
413C
414            SCALE = MAX(ABS(S(M)),ABS(S(M-1)),ABS(E(M-1)),ABS(S(L)),
415     1                    ABS(E(L)))
416            SM = S(M)/SCALE
417            SMM1 = S(M-1)/SCALE
418            EMM1 = E(M-1)/SCALE
419            SL = S(L)/SCALE
420            EL = E(L)/SCALE
421            B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0E0
422            C = (SM*EMM1)**2
423            SHIFT = 0.0E0
424            IF (B .EQ. 0.0E0 .AND. C .EQ. 0.0E0) GO TO 550
425               SHIFT = SQRT(B**2+C)
426               IF (B .LT. 0.0E0) SHIFT = -SHIFT
427               SHIFT = C/(B + SHIFT)
428  550       CONTINUE
429            F = (SL + SM)*(SL - SM) - SHIFT
430            G = SL*EL
431C
432C           CHASE ZEROS.
433C
434            MM1 = M - 1
435            DO 560 K = L, MM1
436               CALL SROTG(F,G,CS,SN)
437               IF (K .NE. L) E(K-1) = F
438               F = CS*S(K) + SN*E(K)
439               E(K) = CS*E(K) - SN*S(K)
440               G = SN*S(K+1)
441               S(K+1) = CS*S(K+1)
442               IF (WANTV) CALL SROT(P,V(1,K),1,V(1,K+1),1,CS,SN)
443               CALL SROTG(F,G,CS,SN)
444               S(K) = F
445               F = CS*E(K) + SN*S(K+1)
446               S(K+1) = -SN*E(K) + CS*S(K+1)
447               G = SN*E(K+1)
448               E(K+1) = CS*E(K+1)
449               IF (WANTU .AND. K .LT. N)
450     1            CALL SROT(N,U(1,K),1,U(1,K+1),1,CS,SN)
451  560       CONTINUE
452            E(M-1) = F
453            ITER = ITER + 1
454         GO TO 610
455C
456C        CONVERGENCE.
457C
458  570    CONTINUE
459C
460C           MAKE THE SINGULAR VALUE  POSITIVE.
461C
462            IF (S(L) .GE. 0.0E0) GO TO 580
463               S(L) = -S(L)
464               IF (WANTV) CALL SSCAL(P,-1.0E0,V(1,L),1)
465  580       CONTINUE
466C
467C           ORDER THE SINGULAR VALUE.
468C
469  590       IF (L .EQ. MM) GO TO 600
470               IF (S(L) .GE. S(L+1)) GO TO 600
471               T = S(L)
472               S(L) = S(L+1)
473               S(L+1) = T
474               IF (WANTV .AND. L .LT. P)
475     1            CALL SSWAP(P,V(1,L),1,V(1,L+1),1)
476               IF (WANTU .AND. L .LT. N)
477     1            CALL SSWAP(N,U(1,L),1,U(1,L+1),1)
478               L = L + 1
479            GO TO 590
480  600       CONTINUE
481            ITER = 0
482            M = M - 1
483  610    CONTINUE
484      GO TO 360
485  620 CONTINUE
486      RETURN
487      END
488