001: SUBROUTINE SBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, 002: $ WORK, IWORK, INFO ) 003: * 004: * -- LAPACK routine (version 3.2) -- 005: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 006: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 007: * November 2006 008: * 009: * .. Scalar Arguments .. 010: CHARACTER COMPQ, UPLO 011: INTEGER INFO, LDU, LDVT, N 012: * .. 013: * .. Array Arguments .. 014: INTEGER IQ( * ), IWORK( * ) 015: REAL D( * ), E( * ), Q( * ), U( LDU, * ), 016: $ VT( LDVT, * ), WORK( * ) 017: * .. 018: * 019: * Purpose 020: * ======= 021: * 022: * SBDSDC computes the singular value decomposition (SVD) of a real 023: * N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT, 024: * using a divide and conquer method, where S is a diagonal matrix 025: * with non-negative diagonal elements (the singular values of B), and 026: * U and VT are orthogonal matrices of left and right singular vectors, 027: * respectively. SBDSDC can be used to compute all singular values, 028: * and optionally, singular vectors or singular vectors in compact form. 029: * 030: * This code makes very mild assumptions about floating point 031: * arithmetic. It will work on machines with a guard digit in 032: * add/subtract, or on those binary machines without guard digits 033: * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. 034: * It could conceivably fail on hexadecimal or decimal machines 035: * without guard digits, but we know of none. See SLASD3 for details. 036: * 037: * The code currently calls SLASDQ if singular values only are desired. 038: * However, it can be slightly modified to compute singular values 039: * using the divide and conquer method. 040: * 041: * Arguments 042: * ========= 043: * 044: * UPLO (input) CHARACTER*1 045: * = 'U': B is upper bidiagonal. 046: * = 'L': B is lower bidiagonal. 047: * 048: * COMPQ (input) CHARACTER*1 049: * Specifies whether singular vectors are to be computed 050: * as follows: 051: * = 'N': Compute singular values only; 052: * = 'P': Compute singular values and compute singular 053: * vectors in compact form; 054: * = 'I': Compute singular values and singular vectors. 055: * 056: * N (input) INTEGER 057: * The order of the matrix B. N >= 0. 058: * 059: * D (input/output) REAL array, dimension (N) 060: * On entry, the n diagonal elements of the bidiagonal matrix B. 061: * On exit, if INFO=0, the singular values of B. 062: * 063: * E (input/output) REAL array, dimension (N-1) 064: * On entry, the elements of E contain the offdiagonal 065: * elements of the bidiagonal matrix whose SVD is desired. 066: * On exit, E has been destroyed. 067: * 068: * U (output) REAL array, dimension (LDU,N) 069: * If COMPQ = 'I', then: 070: * On exit, if INFO = 0, U contains the left singular vectors 071: * of the bidiagonal matrix. 072: * For other values of COMPQ, U is not referenced. 073: * 074: * LDU (input) INTEGER 075: * The leading dimension of the array U. LDU >= 1. 076: * If singular vectors are desired, then LDU >= max( 1, N ). 077: * 078: * VT (output) REAL array, dimension (LDVT,N) 079: * If COMPQ = 'I', then: 080: * On exit, if INFO = 0, VT' contains the right singular 081: * vectors of the bidiagonal matrix. 082: * For other values of COMPQ, VT is not referenced. 083: * 084: * LDVT (input) INTEGER 085: * The leading dimension of the array VT. LDVT >= 1. 086: * If singular vectors are desired, then LDVT >= max( 1, N ). 087: * 088: * Q (output) REAL array, dimension (LDQ) 089: * If COMPQ = 'P', then: 090: * On exit, if INFO = 0, Q and IQ contain the left 091: * and right singular vectors in a compact form, 092: * requiring O(N log N) space instead of 2*N**2. 093: * In particular, Q contains all the REAL data in 094: * LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1)))) 095: * words of memory, where SMLSIZ is returned by ILAENV and 096: * is equal to the maximum size of the subproblems at the 097: * bottom of the computation tree (usually about 25). 098: * For other values of COMPQ, Q is not referenced. 099: * 100: * IQ (output) INTEGER array, dimension (LDIQ) 101: * If COMPQ = 'P', then: 102: * On exit, if INFO = 0, Q and IQ contain the left 103: * and right singular vectors in a compact form, 104: * requiring O(N log N) space instead of 2*N**2. 105: * In particular, IQ contains all INTEGER data in 106: * LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1)))) 107: * words of memory, where SMLSIZ is returned by ILAENV and 108: * is equal to the maximum size of the subproblems at the 109: * bottom of the computation tree (usually about 25). 110: * For other values of COMPQ, IQ is not referenced. 111: * 112: * WORK (workspace) REAL array, dimension (MAX(1,LWORK)) 113: * If COMPQ = 'N' then LWORK >= (4 * N). 114: * If COMPQ = 'P' then LWORK >= (6 * N). 115: * If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N). 116: * 117: * IWORK (workspace) INTEGER array, dimension (8*N) 118: * 119: * INFO (output) INTEGER 120: * = 0: successful exit. 121: * < 0: if INFO = -i, the i-th argument had an illegal value. 122: * > 0: The algorithm failed to compute an singular value. 123: * The update process of divide and conquer failed. 124: * 125: * Further Details 126: * =============== 127: * 128: * Based on contributions by 129: * Ming Gu and Huan Ren, Computer Science Division, University of 130: * California at Berkeley, USA 131: * ===================================================================== 132: * Changed dimension statement in comment describing E from (N) to 133: * (N-1). Sven, 17 Feb 05. 134: * ===================================================================== 135: * 136: * .. Parameters .. 137: REAL ZERO, ONE, TWO 138: PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) 139: * .. 140: * .. Local Scalars .. 141: INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC, 142: $ ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK, 143: $ MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ, 144: $ SMLSZP, SQRE, START, WSTART, Z 145: REAL CS, EPS, ORGNRM, P, R, SN 146: * .. 147: * .. External Functions .. 148: LOGICAL LSAME 149: INTEGER ILAENV 150: REAL SLAMCH, SLANST 151: EXTERNAL SLAMCH, SLANST, ILAENV, LSAME 152: * .. 153: * .. External Subroutines .. 154: EXTERNAL SCOPY, SLARTG, SLASCL, SLASD0, SLASDA, SLASDQ, 155: $ SLASET, SLASR, SSWAP, XERBLA 156: * .. 157: * .. Intrinsic Functions .. 158: INTRINSIC REAL, ABS, INT, LOG, SIGN 159: * .. 160: * .. Executable Statements .. 161: * 162: * Test the input parameters. 163: * 164: INFO = 0 165: * 166: IUPLO = 0 167: IF( LSAME( UPLO, 'U' ) ) 168: $ IUPLO = 1 169: IF( LSAME( UPLO, 'L' ) ) 170: $ IUPLO = 2 171: IF( LSAME( COMPQ, 'N' ) ) THEN 172: ICOMPQ = 0 173: ELSE IF( LSAME( COMPQ, 'P' ) ) THEN 174: ICOMPQ = 1 175: ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 176: ICOMPQ = 2 177: ELSE 178: ICOMPQ = -1 179: END IF 180: IF( IUPLO.EQ.0 ) THEN 181: INFO = -1 182: ELSE IF( ICOMPQ.LT.0 ) THEN 183: INFO = -2 184: ELSE IF( N.LT.0 ) THEN 185: INFO = -3 186: ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT. 187: $ N ) ) ) THEN 188: INFO = -7 189: ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT. 190: $ N ) ) ) THEN 191: INFO = -9 192: END IF 193: IF( INFO.NE.0 ) THEN 194: CALL XERBLA( 'SBDSDC', -INFO ) 195: RETURN 196: END IF 197: * 198: * Quick return if possible 199: * 200: IF( N.EQ.0 ) 201: $ RETURN 202: SMLSIZ = ILAENV( 9, 'SBDSDC', ' ', 0, 0, 0, 0 ) 203: IF( N.EQ.1 ) THEN 204: IF( ICOMPQ.EQ.1 ) THEN 205: Q( 1 ) = SIGN( ONE, D( 1 ) ) 206: Q( 1+SMLSIZ*N ) = ONE 207: ELSE IF( ICOMPQ.EQ.2 ) THEN 208: U( 1, 1 ) = SIGN( ONE, D( 1 ) ) 209: VT( 1, 1 ) = ONE 210: END IF 211: D( 1 ) = ABS( D( 1 ) ) 212: RETURN 213: END IF 214: NM1 = N - 1 215: * 216: * If matrix lower bidiagonal, rotate to be upper bidiagonal 217: * by applying Givens rotations on the left 218: * 219: WSTART = 1 220: QSTART = 3 221: IF( ICOMPQ.EQ.1 ) THEN 222: CALL SCOPY( N, D, 1, Q( 1 ), 1 ) 223: CALL SCOPY( N-1, E, 1, Q( N+1 ), 1 ) 224: END IF 225: IF( IUPLO.EQ.2 ) THEN 226: QSTART = 5 227: WSTART = 2*N - 1 228: DO 10 I = 1, N - 1 229: CALL SLARTG( D( I ), E( I ), CS, SN, R ) 230: D( I ) = R 231: E( I ) = SN*D( I+1 ) 232: D( I+1 ) = CS*D( I+1 ) 233: IF( ICOMPQ.EQ.1 ) THEN 234: Q( I+2*N ) = CS 235: Q( I+3*N ) = SN 236: ELSE IF( ICOMPQ.EQ.2 ) THEN 237: WORK( I ) = CS 238: WORK( NM1+I ) = -SN 239: END IF 240: 10 CONTINUE 241: END IF 242: * 243: * If ICOMPQ = 0, use SLASDQ to compute the singular values. 244: * 245: IF( ICOMPQ.EQ.0 ) THEN 246: CALL SLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U, 247: $ LDU, WORK( WSTART ), INFO ) 248: GO TO 40 249: END IF 250: * 251: * If N is smaller than the minimum divide size SMLSIZ, then solve 252: * the problem with another solver. 253: * 254: IF( N.LE.SMLSIZ ) THEN 255: IF( ICOMPQ.EQ.2 ) THEN 256: CALL SLASET( 'A', N, N, ZERO, ONE, U, LDU ) 257: CALL SLASET( 'A', N, N, ZERO, ONE, VT, LDVT ) 258: CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U, 259: $ LDU, WORK( WSTART ), INFO ) 260: ELSE IF( ICOMPQ.EQ.1 ) THEN 261: IU = 1 262: IVT = IU + N 263: CALL SLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ), 264: $ N ) 265: CALL SLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ), 266: $ N ) 267: CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, 268: $ Q( IVT+( QSTART-1 )*N ), N, 269: $ Q( IU+( QSTART-1 )*N ), N, 270: $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ), 271: $ INFO ) 272: END IF 273: GO TO 40 274: END IF 275: * 276: IF( ICOMPQ.EQ.2 ) THEN 277: CALL SLASET( 'A', N, N, ZERO, ONE, U, LDU ) 278: CALL SLASET( 'A', N, N, ZERO, ONE, VT, LDVT ) 279: END IF 280: * 281: * Scale. 282: * 283: ORGNRM = SLANST( 'M', N, D, E ) 284: IF( ORGNRM.EQ.ZERO ) 285: $ RETURN 286: CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR ) 287: CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR ) 288: * 289: EPS = SLAMCH( 'Epsilon' ) 290: * 291: MLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 292: SMLSZP = SMLSIZ + 1 293: * 294: IF( ICOMPQ.EQ.1 ) THEN 295: IU = 1 296: IVT = 1 + SMLSIZ 297: DIFL = IVT + SMLSZP 298: DIFR = DIFL + MLVL 299: Z = DIFR + MLVL*2 300: IC = Z + MLVL 301: IS = IC + 1 302: POLES = IS + 1 303: GIVNUM = POLES + 2*MLVL 304: * 305: K = 1 306: GIVPTR = 2 307: PERM = 3 308: GIVCOL = PERM + MLVL 309: END IF 310: * 311: DO 20 I = 1, N 312: IF( ABS( D( I ) ).LT.EPS ) THEN 313: D( I ) = SIGN( EPS, D( I ) ) 314: END IF 315: 20 CONTINUE 316: * 317: START = 1 318: SQRE = 0 319: * 320: DO 30 I = 1, NM1 321: IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN 322: * 323: * Subproblem found. First determine its size and then 324: * apply divide and conquer on it. 325: * 326: IF( I.LT.NM1 ) THEN 327: * 328: * A subproblem with E(I) small for I < NM1. 329: * 330: NSIZE = I - START + 1 331: ELSE IF( ABS( E( I ) ).GE.EPS ) THEN 332: * 333: * A subproblem with E(NM1) not too small but I = NM1. 334: * 335: NSIZE = N - START + 1 336: ELSE 337: * 338: * A subproblem with E(NM1) small. This implies an 339: * 1-by-1 subproblem at D(N). Solve this 1-by-1 problem 340: * first. 341: * 342: NSIZE = I - START + 1 343: IF( ICOMPQ.EQ.2 ) THEN 344: U( N, N ) = SIGN( ONE, D( N ) ) 345: VT( N, N ) = ONE 346: ELSE IF( ICOMPQ.EQ.1 ) THEN 347: Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) ) 348: Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE 349: END IF 350: D( N ) = ABS( D( N ) ) 351: END IF 352: IF( ICOMPQ.EQ.2 ) THEN 353: CALL SLASD0( NSIZE, SQRE, D( START ), E( START ), 354: $ U( START, START ), LDU, VT( START, START ), 355: $ LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO ) 356: ELSE 357: CALL SLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ), 358: $ E( START ), Q( START+( IU+QSTART-2 )*N ), N, 359: $ Q( START+( IVT+QSTART-2 )*N ), 360: $ IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )* 361: $ N ), Q( START+( DIFR+QSTART-2 )*N ), 362: $ Q( START+( Z+QSTART-2 )*N ), 363: $ Q( START+( POLES+QSTART-2 )*N ), 364: $ IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ), 365: $ N, IQ( START+PERM*N ), 366: $ Q( START+( GIVNUM+QSTART-2 )*N ), 367: $ Q( START+( IC+QSTART-2 )*N ), 368: $ Q( START+( IS+QSTART-2 )*N ), 369: $ WORK( WSTART ), IWORK, INFO ) 370: IF( INFO.NE.0 ) THEN 371: RETURN 372: END IF 373: END IF 374: START = I + 1 375: END IF 376: 30 CONTINUE 377: * 378: * Unscale 379: * 380: CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR ) 381: 40 CONTINUE 382: * 383: * Use Selection Sort to minimize swaps of singular vectors 384: * 385: DO 60 II = 2, N 386: I = II - 1 387: KK = I 388: P = D( I ) 389: DO 50 J = II, N 390: IF( D( J ).GT.P ) THEN 391: KK = J 392: P = D( J ) 393: END IF 394: 50 CONTINUE 395: IF( KK.NE.I ) THEN 396: D( KK ) = D( I ) 397: D( I ) = P 398: IF( ICOMPQ.EQ.1 ) THEN 399: IQ( I ) = KK 400: ELSE IF( ICOMPQ.EQ.2 ) THEN 401: CALL SSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 ) 402: CALL SSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT ) 403: END IF 404: ELSE IF( ICOMPQ.EQ.1 ) THEN 405: IQ( I ) = I 406: END IF 407: 60 CONTINUE 408: * 409: * If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO 410: * 411: IF( ICOMPQ.EQ.1 ) THEN 412: IF( IUPLO.EQ.1 ) THEN 413: IQ( N ) = 1 414: ELSE 415: IQ( N ) = 0 416: END IF 417: END IF 418: * 419: * If B is lower bidiagonal, update U by those Givens rotations 420: * which rotated B to be upper bidiagonal 421: * 422: IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) ) 423: $ CALL SLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU ) 424: * 425: RETURN 426: * 427: * End of SBDSDC 428: * 429: END 430: