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