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