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