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