Annotation of rpl/lapack/lapack/dbdsdc.f, revision 1.1
1.1 ! bertrand 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) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 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 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: * =====================================================================
! 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: IF( INFO.NE.0 ) THEN
! 372: RETURN
! 373: END IF
! 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dbdsdc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack