1 SUBROUTINE SBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, 2 $ WORK, IWORK, INFO ) 3* 4* -- LAPACK routine (instrumented to count ops, version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* October 31, 1999 8* 9* .. Scalar Arguments .. 10 CHARACTER COMPQ, UPLO 11 INTEGER INFO, LDU, LDVT, N 12* .. 13* .. Array Arguments .. 14 INTEGER IQ( * ), IWORK( * ) 15 REAL D( * ), E( * ), Q( * ), U( LDU, * ), 16 $ VT( LDVT, * ), WORK( * ) 17* .. 18* .. Common blocks .. 19 COMMON / LATIME / OPS, ITCNT 20* .. 21* .. Scalars in Common .. 22 REAL ITCNT, OPS 23* .. 24* 25* Purpose 26* ======= 27* 28* SBDSDC computes the singular value decomposition (SVD) of a real 29* N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT, 30* using a divide and conquer method, where S is a diagonal matrix 31* with non-negative diagonal elements (the singular values of B), and 32* U and VT are orthogonal matrices of left and right singular vectors, 33* respectively. SBDSDC can be used to compute all singular values, 34* and optionally, singular vectors or singular vectors in compact form. 35* 36* This code makes very mild assumptions about floating point 37* arithmetic. It will work on machines with a guard digit in 38* add/subtract, or on those binary machines without guard digits 39* which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. 40* It could conceivably fail on hexadecimal or decimal machines 41* without guard digits, but we know of none. See SLASD3 for details. 42* 43* The code currently call SLASDQ if singular values only are desired. 44* However, it can be slightly modified to compute singular values 45* using the divide and conquer method. 46* 47* Arguments 48* ========= 49* 50* UPLO (input) CHARACTER*1 51* = 'U': B is upper bidiagonal. 52* = 'L': B is lower bidiagonal. 53* 54* COMPQ (input) CHARACTER*1 55* Specifies whether singular vectors are to be computed 56* as follows: 57* = 'N': Compute singular values only; 58* = 'P': Compute singular values and compute singular 59* vectors in compact form; 60* = 'I': Compute singular values and singular vectors. 61* 62* N (input) INTEGER 63* The order of the matrix B. N >= 0. 64* 65* D (input/output) REAL array, dimension (N) 66* On entry, the n diagonal elements of the bidiagonal matrix B. 67* On exit, if INFO=0, the singular values of B. 68* 69* E (input/output) REAL array, dimension (N) 70* On entry, the elements of E contain the offdiagonal 71* elements of the bidiagonal matrix whose SVD is desired. 72* On exit, E has been destroyed. 73* 74* U (output) REAL array, dimension (LDU,N) 75* If COMPQ = 'I', then: 76* On exit, if INFO = 0, U contains the left singular vectors 77* of the bidiagonal matrix. 78* For other values of COMPQ, U is not referenced. 79* 80* LDU (input) INTEGER 81* The leading dimension of the array U. LDU >= 1. 82* If singular vectors are desired, then LDU >= max( 1, N ). 83* 84* VT (output) REAL array, dimension (LDVT,N) 85* If COMPQ = 'I', then: 86* On exit, if INFO = 0, VT' contains the right singular 87* vectors of the bidiagonal matrix. 88* For other values of COMPQ, VT is not referenced. 89* 90* LDVT (input) INTEGER 91* The leading dimension of the array VT. LDVT >= 1. 92* If singular vectors are desired, then LDVT >= max( 1, N ). 93* 94* Q (output) REAL array, dimension (LDQ) 95* If COMPQ = 'P', then: 96* On exit, if INFO = 0, Q and IQ contain the left 97* and right singular vectors in a compact form, 98* requiring O(N log N) space instead of 2*N**2. 99* In particular, Q contains all the REAL data in 100* LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1)))) 101* words of memory, where SMLSIZ is returned by ILAENV and 102* is equal to the maximum size of the subproblems at the 103* bottom of the computation tree (usually about 25). 104* For other values of COMPQ, Q is not referenced. 105* 106* IQ (output) INTEGER array, dimension (LDIQ) 107* If COMPQ = 'P', then: 108* On exit, if INFO = 0, Q and IQ contain the left 109* and right singular vectors in a compact form, 110* requiring O(N log N) space instead of 2*N**2. 111* In particular, IQ contains all INTEGER data in 112* LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1)))) 113* words of memory, where SMLSIZ is returned by ILAENV and 114* is equal to the maximum size of the subproblems at the 115* bottom of the computation tree (usually about 25). 116* For other values of COMPQ, IQ is not referenced. 117* 118* WORK (workspace) REAL array, dimension (LWORK) 119* If COMPQ = 'N' then LWORK >= (4 * N). 120* If COMPQ = 'P' then LWORK >= (6 * N). 121* If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N). 122* 123* IWORK (workspace) INTEGER array, dimension (7*N) 124* 125* INFO (output) INTEGER 126* = 0: successful exit. 127* < 0: if INFO = -i, the i-th argument had an illegal value. 128* > 0: The algorithm failed to compute an singular value. 129* The update process of divide and conquer failed. 130* 131* Further Details 132* =============== 133* 134* Based on contributions by 135* Ming Gu and Huan Ren, Computer Science Division, University of 136* California at Berkeley, USA 137* 138* ===================================================================== 139* 140* .. Parameters .. 141 REAL ZERO, ONE, TWO 142 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) 143* .. 144* .. Local Scalars .. 145 INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC, 146 $ ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK, 147 $ MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ, 148 $ SMLSZP, SQRE, START, WSTART, Z 149 REAL CS, EPS, ORGNRM, P, R, SN 150* .. 151* .. External Functions .. 152 LOGICAL LSAME 153 INTEGER ILAENV 154 REAL SLAMCH, SLANST 155 EXTERNAL SLAMCH, SLANST, ILAENV, LSAME 156* .. 157* .. External Subroutines .. 158 EXTERNAL SCOPY, SLARTG, SLASCL, SLASD0, SLASDA, SLASDQ, 159 $ SLASET, SLASR, SSWAP, XERBLA 160* .. 161* .. Intrinsic Functions .. 162 INTRINSIC REAL, ABS, INT, LOG, SIGN 163* .. 164* .. Executable Statements .. 165* 166* Test the input parameters. 167* 168 INFO = 0 169* 170 IUPLO = 0 171 IF( LSAME( UPLO, 'U' ) ) 172 $ IUPLO = 1 173 IF( LSAME( UPLO, 'L' ) ) 174 $ IUPLO = 2 175 IF( LSAME( COMPQ, 'N' ) ) THEN 176 ICOMPQ = 0 177 ELSE IF( LSAME( COMPQ, 'P' ) ) THEN 178 ICOMPQ = 1 179 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 180 ICOMPQ = 2 181 ELSE 182 ICOMPQ = -1 183 END IF 184 IF( IUPLO.EQ.0 ) THEN 185 INFO = -1 186 ELSE IF( ICOMPQ.LT.0 ) THEN 187 INFO = -2 188 ELSE IF( N.LT.0 ) THEN 189 INFO = -3 190 ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT. 191 $ N ) ) ) THEN 192 INFO = -7 193 ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT. 194 $ N ) ) ) THEN 195 INFO = -9 196 END IF 197 IF( INFO.NE.0 ) THEN 198 CALL XERBLA( 'SBDSDC', -INFO ) 199 RETURN 200 END IF 201* 202* Quick return if possible 203* 204 IF( N.EQ.0 ) 205 $ RETURN 206 SMLSIZ = ILAENV( 9, 'SBDSDC', ' ', 0, 0, 0, 0 ) 207 IF( N.EQ.1 ) THEN 208 IF( ICOMPQ.EQ.1 ) THEN 209 Q( 1 ) = SIGN( ONE, D( 1 ) ) 210 Q( 1+SMLSIZ*N ) = ONE 211 ELSE IF( ICOMPQ.EQ.2 ) THEN 212 U( 1, 1 ) = SIGN( ONE, D( 1 ) ) 213 VT( 1, 1 ) = ONE 214 END IF 215 D( 1 ) = ABS( D( 1 ) ) 216 RETURN 217 END IF 218 NM1 = N - 1 219* 220* If matrix lower bidiagonal, rotate to be upper bidiagonal 221* by applying Givens rotations on the left 222* 223 WSTART = 1 224 QSTART = 3 225 IF( ICOMPQ.EQ.1 ) THEN 226 CALL SCOPY( N, D, 1, Q( 1 ), 1 ) 227 CALL SCOPY( N-1, E, 1, Q( N+1 ), 1 ) 228 END IF 229 IF( IUPLO.EQ.2 ) THEN 230 QSTART = 5 231 WSTART = 2*N - 1 232 OPS = OPS + REAL( 8*( N-1 ) ) 233 DO 10 I = 1, N - 1 234 CALL SLARTG( D( I ), E( I ), CS, SN, R ) 235 D( I ) = R 236 E( I ) = SN*D( I+1 ) 237 D( I+1 ) = CS*D( I+1 ) 238 IF( ICOMPQ.EQ.1 ) THEN 239 Q( I+2*N ) = CS 240 Q( I+3*N ) = SN 241 ELSE IF( ICOMPQ.EQ.2 ) THEN 242 WORK( I ) = CS 243 WORK( NM1+I ) = -SN 244 END IF 245 10 CONTINUE 246 END IF 247* 248* If ICOMPQ = 0, use SLASDQ to compute the singular values. 249* 250 IF( ICOMPQ.EQ.0 ) THEN 251 CALL SLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U, 252 $ LDU, WORK( WSTART ), INFO ) 253 GO TO 40 254 END IF 255* 256* If N is smaller than the minimum divide size SMLSIZ, then solve 257* the problem with another solver. 258* 259 IF( N.LE.SMLSIZ ) THEN 260 IF( ICOMPQ.EQ.2 ) THEN 261 CALL SLASET( 'A', N, N, ZERO, ONE, U, LDU ) 262 CALL SLASET( 'A', N, N, ZERO, ONE, VT, LDVT ) 263 CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U, 264 $ LDU, WORK( WSTART ), INFO ) 265 ELSE IF( ICOMPQ.EQ.1 ) THEN 266 IU = 1 267 IVT = IU + N 268 CALL SLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ), 269 $ N ) 270 CALL SLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ), 271 $ N ) 272 CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, 273 $ Q( IVT+( QSTART-1 )*N ), N, 274 $ Q( IU+( QSTART-1 )*N ), N, 275 $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ), 276 $ INFO ) 277 END IF 278 GO TO 40 279 END IF 280* 281 IF( ICOMPQ.EQ.2 ) THEN 282 CALL SLASET( 'A', N, N, ZERO, ONE, U, LDU ) 283 CALL SLASET( 'A', N, N, ZERO, ONE, VT, LDVT ) 284 END IF 285* 286* Scale. 287* 288 ORGNRM = SLANST( 'M', N, D, E ) 289 IF( ORGNRM.EQ.ZERO ) 290 $ RETURN 291 OPS = OPS + REAL( N + NM1 ) 292 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR ) 293 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR ) 294* 295 EPS = SLAMCH( 'Epsilon' ) 296* 297 MLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 298 SMLSZP = SMLSIZ + 1 299* 300 IF( ICOMPQ.EQ.1 ) THEN 301 IU = 1 302 IVT = 1 + SMLSIZ 303 DIFL = IVT + SMLSZP 304 DIFR = DIFL + MLVL 305 Z = DIFR + MLVL*2 306 IC = Z + MLVL 307 IS = IC + 1 308 POLES = IS + 1 309 GIVNUM = POLES + 2*MLVL 310* 311 K = 1 312 GIVPTR = 2 313 PERM = 3 314 GIVCOL = PERM + MLVL 315 END IF 316* 317 DO 20 I = 1, N 318 IF( ABS( D( I ) ).LT.EPS ) THEN 319 D( I ) = SIGN( EPS, D( I ) ) 320 END IF 321 20 CONTINUE 322* 323 START = 1 324 SQRE = 0 325* 326 DO 30 I = 1, NM1 327 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN 328* 329* Subproblem found. First determine its size and then 330* apply divide and conquer on it. 331* 332 IF( I.LT.NM1 ) THEN 333* 334* A subproblem with E(I) small for I < NM1. 335* 336 NSIZE = I - START + 1 337 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN 338* 339* A subproblem with E(NM1) not too small but I = NM1. 340* 341 NSIZE = N - START + 1 342 ELSE 343* 344* A subproblem with E(NM1) small. This implies an 345* 1-by-1 subproblem at D(N). Solve this 1-by-1 problem 346* first. 347* 348 NSIZE = I - START + 1 349 IF( ICOMPQ.EQ.2 ) THEN 350 U( N, N ) = SIGN( ONE, D( N ) ) 351 VT( N, N ) = ONE 352 ELSE IF( ICOMPQ.EQ.1 ) THEN 353 Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) ) 354 Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE 355 END IF 356 D( N ) = ABS( D( N ) ) 357 END IF 358 IF( ICOMPQ.EQ.2 ) THEN 359 CALL SLASD0( NSIZE, SQRE, D( START ), E( START ), 360 $ U( START, START ), LDU, VT( START, START ), 361 $ LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO ) 362 ELSE 363 CALL SLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ), 364 $ E( START ), Q( START+( IU+QSTART-2 )*N ), N, 365 $ Q( START+( IVT+QSTART-2 )*N ), 366 $ IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )* 367 $ N ), Q( START+( DIFR+QSTART-2 )*N ), 368 $ Q( START+( Z+QSTART-2 )*N ), 369 $ Q( START+( POLES+QSTART-2 )*N ), 370 $ IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ), 371 $ N, IQ( START+PERM*N ), 372 $ Q( START+( GIVNUM+QSTART-2 )*N ), 373 $ Q( START+( IC+QSTART-2 )*N ), 374 $ Q( START+( IS+QSTART-2 )*N ), 375 $ WORK( WSTART ), IWORK, INFO ) 376 IF( INFO.NE.0 ) THEN 377 RETURN 378 END IF 379 END IF 380 START = I + 1 381 END IF 382 30 CONTINUE 383* 384* Unscale 385* 386 OPS = OPS + REAL( N ) 387 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR ) 388 40 CONTINUE 389* 390* Use Selection Sort to minimize swaps of singular vectors 391* 392 DO 60 II = 2, N 393 I = II - 1 394 KK = I 395 P = D( I ) 396 DO 50 J = II, N 397 IF( D( J ).GT.P ) THEN 398 KK = J 399 P = D( J ) 400 END IF 401 50 CONTINUE 402 IF( KK.NE.I ) THEN 403 D( KK ) = D( I ) 404 D( I ) = P 405 IF( ICOMPQ.EQ.1 ) THEN 406 IQ( I ) = KK 407 ELSE IF( ICOMPQ.EQ.2 ) THEN 408 CALL SSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 ) 409 CALL SSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT ) 410 END IF 411 ELSE IF( ICOMPQ.EQ.1 ) THEN 412 IQ( I ) = I 413 END IF 414 60 CONTINUE 415* 416* If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO 417* 418 IF( ICOMPQ.EQ.1 ) THEN 419 IF( IUPLO.EQ.1 ) THEN 420 IQ( N ) = 1 421 ELSE 422 IQ( N ) = 0 423 END IF 424 END IF 425* 426* If B is lower bidiagonal, update U by those Givens rotations 427* which rotated B to be upper bidiagonal 428* 429 IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) ) THEN 430 OPS = OPS + REAL( 6*( N-1 )*N ) 431 CALL SLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU ) 432 END IF 433* 434 RETURN 435* 436* End of SBDSDC 437* 438 END 439