Annotation of rpl/lapack/lapack/dlasd2.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLASD2( NL, NR, SQRE, K, D, Z, ALPHA, BETA, U, LDU, VT,
! 2: $ LDVT, DSIGMA, U2, LDU2, VT2, LDVT2, IDXP, IDX,
! 3: $ IDXC, IDXQ, COLTYP, INFO )
! 4: *
! 5: * -- LAPACK auxiliary routine (version 3.2) --
! 6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 8: * November 2006
! 9: *
! 10: * .. Scalar Arguments ..
! 11: INTEGER INFO, K, LDU, LDU2, LDVT, LDVT2, NL, NR, SQRE
! 12: DOUBLE PRECISION ALPHA, BETA
! 13: * ..
! 14: * .. Array Arguments ..
! 15: INTEGER COLTYP( * ), IDX( * ), IDXC( * ), IDXP( * ),
! 16: $ IDXQ( * )
! 17: DOUBLE PRECISION D( * ), DSIGMA( * ), U( LDU, * ),
! 18: $ U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ),
! 19: $ Z( * )
! 20: * ..
! 21: *
! 22: * Purpose
! 23: * =======
! 24: *
! 25: * DLASD2 merges the two sets of singular values together into a single
! 26: * sorted set. Then it tries to deflate the size of the problem.
! 27: * There are two ways in which deflation can occur: when two or more
! 28: * singular values are close together or if there is a tiny entry in the
! 29: * Z vector. For each such occurrence the order of the related secular
! 30: * equation problem is reduced by one.
! 31: *
! 32: * DLASD2 is called from DLASD1.
! 33: *
! 34: * Arguments
! 35: * =========
! 36: *
! 37: * NL (input) INTEGER
! 38: * The row dimension of the upper block. NL >= 1.
! 39: *
! 40: * NR (input) INTEGER
! 41: * The row dimension of the lower block. NR >= 1.
! 42: *
! 43: * SQRE (input) INTEGER
! 44: * = 0: the lower block is an NR-by-NR square matrix.
! 45: * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
! 46: *
! 47: * The bidiagonal matrix has N = NL + NR + 1 rows and
! 48: * M = N + SQRE >= N columns.
! 49: *
! 50: * K (output) INTEGER
! 51: * Contains the dimension of the non-deflated matrix,
! 52: * This is the order of the related secular equation. 1 <= K <=N.
! 53: *
! 54: * D (input/output) DOUBLE PRECISION array, dimension(N)
! 55: * On entry D contains the singular values of the two submatrices
! 56: * to be combined. On exit D contains the trailing (N-K) updated
! 57: * singular values (those which were deflated) sorted into
! 58: * increasing order.
! 59: *
! 60: * Z (output) DOUBLE PRECISION array, dimension(N)
! 61: * On exit Z contains the updating row vector in the secular
! 62: * equation.
! 63: *
! 64: * ALPHA (input) DOUBLE PRECISION
! 65: * Contains the diagonal element associated with the added row.
! 66: *
! 67: * BETA (input) DOUBLE PRECISION
! 68: * Contains the off-diagonal element associated with the added
! 69: * row.
! 70: *
! 71: * U (input/output) DOUBLE PRECISION array, dimension(LDU,N)
! 72: * On entry U contains the left singular vectors of two
! 73: * submatrices in the two square blocks with corners at (1,1),
! 74: * (NL, NL), and (NL+2, NL+2), (N,N).
! 75: * On exit U contains the trailing (N-K) updated left singular
! 76: * vectors (those which were deflated) in its last N-K columns.
! 77: *
! 78: * LDU (input) INTEGER
! 79: * The leading dimension of the array U. LDU >= N.
! 80: *
! 81: * VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
! 82: * On entry VT' contains the right singular vectors of two
! 83: * submatrices in the two square blocks with corners at (1,1),
! 84: * (NL+1, NL+1), and (NL+2, NL+2), (M,M).
! 85: * On exit VT' contains the trailing (N-K) updated right singular
! 86: * vectors (those which were deflated) in its last N-K columns.
! 87: * In case SQRE =1, the last row of VT spans the right null
! 88: * space.
! 89: *
! 90: * LDVT (input) INTEGER
! 91: * The leading dimension of the array VT. LDVT >= M.
! 92: *
! 93: * DSIGMA (output) DOUBLE PRECISION array, dimension (N)
! 94: * Contains a copy of the diagonal elements (K-1 singular values
! 95: * and one zero) in the secular equation.
! 96: *
! 97: * U2 (output) DOUBLE PRECISION array, dimension(LDU2,N)
! 98: * Contains a copy of the first K-1 left singular vectors which
! 99: * will be used by DLASD3 in a matrix multiply (DGEMM) to solve
! 100: * for the new left singular vectors. U2 is arranged into four
! 101: * blocks. The first block contains a column with 1 at NL+1 and
! 102: * zero everywhere else; the second block contains non-zero
! 103: * entries only at and above NL; the third contains non-zero
! 104: * entries only below NL+1; and the fourth is dense.
! 105: *
! 106: * LDU2 (input) INTEGER
! 107: * The leading dimension of the array U2. LDU2 >= N.
! 108: *
! 109: * VT2 (output) DOUBLE PRECISION array, dimension(LDVT2,N)
! 110: * VT2' contains a copy of the first K right singular vectors
! 111: * which will be used by DLASD3 in a matrix multiply (DGEMM) to
! 112: * solve for the new right singular vectors. VT2 is arranged into
! 113: * three blocks. The first block contains a row that corresponds
! 114: * to the special 0 diagonal element in SIGMA; the second block
! 115: * contains non-zeros only at and before NL +1; the third block
! 116: * contains non-zeros only at and after NL +2.
! 117: *
! 118: * LDVT2 (input) INTEGER
! 119: * The leading dimension of the array VT2. LDVT2 >= M.
! 120: *
! 121: * IDXP (workspace) INTEGER array dimension(N)
! 122: * This will contain the permutation used to place deflated
! 123: * values of D at the end of the array. On output IDXP(2:K)
! 124: * points to the nondeflated D-values and IDXP(K+1:N)
! 125: * points to the deflated singular values.
! 126: *
! 127: * IDX (workspace) INTEGER array dimension(N)
! 128: * This will contain the permutation used to sort the contents of
! 129: * D into ascending order.
! 130: *
! 131: * IDXC (output) INTEGER array dimension(N)
! 132: * This will contain the permutation used to arrange the columns
! 133: * of the deflated U matrix into three groups: the first group
! 134: * contains non-zero entries only at and above NL, the second
! 135: * contains non-zero entries only below NL+2, and the third is
! 136: * dense.
! 137: *
! 138: * IDXQ (input/output) INTEGER array dimension(N)
! 139: * This contains the permutation which separately sorts the two
! 140: * sub-problems in D into ascending order. Note that entries in
! 141: * the first hlaf of this permutation must first be moved one
! 142: * position backward; and entries in the second half
! 143: * must first have NL+1 added to their values.
! 144: *
! 145: * COLTYP (workspace/output) INTEGER array dimension(N)
! 146: * As workspace, this will contain a label which will indicate
! 147: * which of the following types a column in the U2 matrix or a
! 148: * row in the VT2 matrix is:
! 149: * 1 : non-zero in the upper half only
! 150: * 2 : non-zero in the lower half only
! 151: * 3 : dense
! 152: * 4 : deflated
! 153: *
! 154: * On exit, it is an array of dimension 4, with COLTYP(I) being
! 155: * the dimension of the I-th type columns.
! 156: *
! 157: * INFO (output) INTEGER
! 158: * = 0: successful exit.
! 159: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 160: *
! 161: * Further Details
! 162: * ===============
! 163: *
! 164: * Based on contributions by
! 165: * Ming Gu and Huan Ren, Computer Science Division, University of
! 166: * California at Berkeley, USA
! 167: *
! 168: * =====================================================================
! 169: *
! 170: * .. Parameters ..
! 171: DOUBLE PRECISION ZERO, ONE, TWO, EIGHT
! 172: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
! 173: $ EIGHT = 8.0D+0 )
! 174: * ..
! 175: * .. Local Arrays ..
! 176: INTEGER CTOT( 4 ), PSM( 4 )
! 177: * ..
! 178: * .. Local Scalars ..
! 179: INTEGER CT, I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M,
! 180: $ N, NLP1, NLP2
! 181: DOUBLE PRECISION C, EPS, HLFTOL, S, TAU, TOL, Z1
! 182: * ..
! 183: * .. External Functions ..
! 184: DOUBLE PRECISION DLAMCH, DLAPY2
! 185: EXTERNAL DLAMCH, DLAPY2
! 186: * ..
! 187: * .. External Subroutines ..
! 188: EXTERNAL DCOPY, DLACPY, DLAMRG, DLASET, DROT, XERBLA
! 189: * ..
! 190: * .. Intrinsic Functions ..
! 191: INTRINSIC ABS, MAX
! 192: * ..
! 193: * .. Executable Statements ..
! 194: *
! 195: * Test the input parameters.
! 196: *
! 197: INFO = 0
! 198: *
! 199: IF( NL.LT.1 ) THEN
! 200: INFO = -1
! 201: ELSE IF( NR.LT.1 ) THEN
! 202: INFO = -2
! 203: ELSE IF( ( SQRE.NE.1 ) .AND. ( SQRE.NE.0 ) ) THEN
! 204: INFO = -3
! 205: END IF
! 206: *
! 207: N = NL + NR + 1
! 208: M = N + SQRE
! 209: *
! 210: IF( LDU.LT.N ) THEN
! 211: INFO = -10
! 212: ELSE IF( LDVT.LT.M ) THEN
! 213: INFO = -12
! 214: ELSE IF( LDU2.LT.N ) THEN
! 215: INFO = -15
! 216: ELSE IF( LDVT2.LT.M ) THEN
! 217: INFO = -17
! 218: END IF
! 219: IF( INFO.NE.0 ) THEN
! 220: CALL XERBLA( 'DLASD2', -INFO )
! 221: RETURN
! 222: END IF
! 223: *
! 224: NLP1 = NL + 1
! 225: NLP2 = NL + 2
! 226: *
! 227: * Generate the first part of the vector Z; and move the singular
! 228: * values in the first part of D one position backward.
! 229: *
! 230: Z1 = ALPHA*VT( NLP1, NLP1 )
! 231: Z( 1 ) = Z1
! 232: DO 10 I = NL, 1, -1
! 233: Z( I+1 ) = ALPHA*VT( I, NLP1 )
! 234: D( I+1 ) = D( I )
! 235: IDXQ( I+1 ) = IDXQ( I ) + 1
! 236: 10 CONTINUE
! 237: *
! 238: * Generate the second part of the vector Z.
! 239: *
! 240: DO 20 I = NLP2, M
! 241: Z( I ) = BETA*VT( I, NLP2 )
! 242: 20 CONTINUE
! 243: *
! 244: * Initialize some reference arrays.
! 245: *
! 246: DO 30 I = 2, NLP1
! 247: COLTYP( I ) = 1
! 248: 30 CONTINUE
! 249: DO 40 I = NLP2, N
! 250: COLTYP( I ) = 2
! 251: 40 CONTINUE
! 252: *
! 253: * Sort the singular values into increasing order
! 254: *
! 255: DO 50 I = NLP2, N
! 256: IDXQ( I ) = IDXQ( I ) + NLP1
! 257: 50 CONTINUE
! 258: *
! 259: * DSIGMA, IDXC, IDXC, and the first column of U2
! 260: * are used as storage space.
! 261: *
! 262: DO 60 I = 2, N
! 263: DSIGMA( I ) = D( IDXQ( I ) )
! 264: U2( I, 1 ) = Z( IDXQ( I ) )
! 265: IDXC( I ) = COLTYP( IDXQ( I ) )
! 266: 60 CONTINUE
! 267: *
! 268: CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
! 269: *
! 270: DO 70 I = 2, N
! 271: IDXI = 1 + IDX( I )
! 272: D( I ) = DSIGMA( IDXI )
! 273: Z( I ) = U2( IDXI, 1 )
! 274: COLTYP( I ) = IDXC( IDXI )
! 275: 70 CONTINUE
! 276: *
! 277: * Calculate the allowable deflation tolerance
! 278: *
! 279: EPS = DLAMCH( 'Epsilon' )
! 280: TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
! 281: TOL = EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
! 282: *
! 283: * There are 2 kinds of deflation -- first a value in the z-vector
! 284: * is small, second two (or more) singular values are very close
! 285: * together (their difference is small).
! 286: *
! 287: * If the value in the z-vector is small, we simply permute the
! 288: * array so that the corresponding singular value is moved to the
! 289: * end.
! 290: *
! 291: * If two values in the D-vector are close, we perform a two-sided
! 292: * rotation designed to make one of the corresponding z-vector
! 293: * entries zero, and then permute the array so that the deflated
! 294: * singular value is moved to the end.
! 295: *
! 296: * If there are multiple singular values then the problem deflates.
! 297: * Here the number of equal singular values are found. As each equal
! 298: * singular value is found, an elementary reflector is computed to
! 299: * rotate the corresponding singular subspace so that the
! 300: * corresponding components of Z are zero in this new basis.
! 301: *
! 302: K = 1
! 303: K2 = N + 1
! 304: DO 80 J = 2, N
! 305: IF( ABS( Z( J ) ).LE.TOL ) THEN
! 306: *
! 307: * Deflate due to small z component.
! 308: *
! 309: K2 = K2 - 1
! 310: IDXP( K2 ) = J
! 311: COLTYP( J ) = 4
! 312: IF( J.EQ.N )
! 313: $ GO TO 120
! 314: ELSE
! 315: JPREV = J
! 316: GO TO 90
! 317: END IF
! 318: 80 CONTINUE
! 319: 90 CONTINUE
! 320: J = JPREV
! 321: 100 CONTINUE
! 322: J = J + 1
! 323: IF( J.GT.N )
! 324: $ GO TO 110
! 325: IF( ABS( Z( J ) ).LE.TOL ) THEN
! 326: *
! 327: * Deflate due to small z component.
! 328: *
! 329: K2 = K2 - 1
! 330: IDXP( K2 ) = J
! 331: COLTYP( J ) = 4
! 332: ELSE
! 333: *
! 334: * Check if singular values are close enough to allow deflation.
! 335: *
! 336: IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
! 337: *
! 338: * Deflation is possible.
! 339: *
! 340: S = Z( JPREV )
! 341: C = Z( J )
! 342: *
! 343: * Find sqrt(a**2+b**2) without overflow or
! 344: * destructive underflow.
! 345: *
! 346: TAU = DLAPY2( C, S )
! 347: C = C / TAU
! 348: S = -S / TAU
! 349: Z( J ) = TAU
! 350: Z( JPREV ) = ZERO
! 351: *
! 352: * Apply back the Givens rotation to the left and right
! 353: * singular vector matrices.
! 354: *
! 355: IDXJP = IDXQ( IDX( JPREV )+1 )
! 356: IDXJ = IDXQ( IDX( J )+1 )
! 357: IF( IDXJP.LE.NLP1 ) THEN
! 358: IDXJP = IDXJP - 1
! 359: END IF
! 360: IF( IDXJ.LE.NLP1 ) THEN
! 361: IDXJ = IDXJ - 1
! 362: END IF
! 363: CALL DROT( N, U( 1, IDXJP ), 1, U( 1, IDXJ ), 1, C, S )
! 364: CALL DROT( M, VT( IDXJP, 1 ), LDVT, VT( IDXJ, 1 ), LDVT, C,
! 365: $ S )
! 366: IF( COLTYP( J ).NE.COLTYP( JPREV ) ) THEN
! 367: COLTYP( J ) = 3
! 368: END IF
! 369: COLTYP( JPREV ) = 4
! 370: K2 = K2 - 1
! 371: IDXP( K2 ) = JPREV
! 372: JPREV = J
! 373: ELSE
! 374: K = K + 1
! 375: U2( K, 1 ) = Z( JPREV )
! 376: DSIGMA( K ) = D( JPREV )
! 377: IDXP( K ) = JPREV
! 378: JPREV = J
! 379: END IF
! 380: END IF
! 381: GO TO 100
! 382: 110 CONTINUE
! 383: *
! 384: * Record the last singular value.
! 385: *
! 386: K = K + 1
! 387: U2( K, 1 ) = Z( JPREV )
! 388: DSIGMA( K ) = D( JPREV )
! 389: IDXP( K ) = JPREV
! 390: *
! 391: 120 CONTINUE
! 392: *
! 393: * Count up the total number of the various types of columns, then
! 394: * form a permutation which positions the four column types into
! 395: * four groups of uniform structure (although one or more of these
! 396: * groups may be empty).
! 397: *
! 398: DO 130 J = 1, 4
! 399: CTOT( J ) = 0
! 400: 130 CONTINUE
! 401: DO 140 J = 2, N
! 402: CT = COLTYP( J )
! 403: CTOT( CT ) = CTOT( CT ) + 1
! 404: 140 CONTINUE
! 405: *
! 406: * PSM(*) = Position in SubMatrix (of types 1 through 4)
! 407: *
! 408: PSM( 1 ) = 2
! 409: PSM( 2 ) = 2 + CTOT( 1 )
! 410: PSM( 3 ) = PSM( 2 ) + CTOT( 2 )
! 411: PSM( 4 ) = PSM( 3 ) + CTOT( 3 )
! 412: *
! 413: * Fill out the IDXC array so that the permutation which it induces
! 414: * will place all type-1 columns first, all type-2 columns next,
! 415: * then all type-3's, and finally all type-4's, starting from the
! 416: * second column. This applies similarly to the rows of VT.
! 417: *
! 418: DO 150 J = 2, N
! 419: JP = IDXP( J )
! 420: CT = COLTYP( JP )
! 421: IDXC( PSM( CT ) ) = J
! 422: PSM( CT ) = PSM( CT ) + 1
! 423: 150 CONTINUE
! 424: *
! 425: * Sort the singular values and corresponding singular vectors into
! 426: * DSIGMA, U2, and VT2 respectively. The singular values/vectors
! 427: * which were not deflated go into the first K slots of DSIGMA, U2,
! 428: * and VT2 respectively, while those which were deflated go into the
! 429: * last N - K slots, except that the first column/row will be treated
! 430: * separately.
! 431: *
! 432: DO 160 J = 2, N
! 433: JP = IDXP( J )
! 434: DSIGMA( J ) = D( JP )
! 435: IDXJ = IDXQ( IDX( IDXP( IDXC( J ) ) )+1 )
! 436: IF( IDXJ.LE.NLP1 ) THEN
! 437: IDXJ = IDXJ - 1
! 438: END IF
! 439: CALL DCOPY( N, U( 1, IDXJ ), 1, U2( 1, J ), 1 )
! 440: CALL DCOPY( M, VT( IDXJ, 1 ), LDVT, VT2( J, 1 ), LDVT2 )
! 441: 160 CONTINUE
! 442: *
! 443: * Determine DSIGMA(1), DSIGMA(2) and Z(1)
! 444: *
! 445: DSIGMA( 1 ) = ZERO
! 446: HLFTOL = TOL / TWO
! 447: IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
! 448: $ DSIGMA( 2 ) = HLFTOL
! 449: IF( M.GT.N ) THEN
! 450: Z( 1 ) = DLAPY2( Z1, Z( M ) )
! 451: IF( Z( 1 ).LE.TOL ) THEN
! 452: C = ONE
! 453: S = ZERO
! 454: Z( 1 ) = TOL
! 455: ELSE
! 456: C = Z1 / Z( 1 )
! 457: S = Z( M ) / Z( 1 )
! 458: END IF
! 459: ELSE
! 460: IF( ABS( Z1 ).LE.TOL ) THEN
! 461: Z( 1 ) = TOL
! 462: ELSE
! 463: Z( 1 ) = Z1
! 464: END IF
! 465: END IF
! 466: *
! 467: * Move the rest of the updating row to Z.
! 468: *
! 469: CALL DCOPY( K-1, U2( 2, 1 ), 1, Z( 2 ), 1 )
! 470: *
! 471: * Determine the first column of U2, the first row of VT2 and the
! 472: * last row of VT.
! 473: *
! 474: CALL DLASET( 'A', N, 1, ZERO, ZERO, U2, LDU2 )
! 475: U2( NLP1, 1 ) = ONE
! 476: IF( M.GT.N ) THEN
! 477: DO 170 I = 1, NLP1
! 478: VT( M, I ) = -S*VT( NLP1, I )
! 479: VT2( 1, I ) = C*VT( NLP1, I )
! 480: 170 CONTINUE
! 481: DO 180 I = NLP2, M
! 482: VT2( 1, I ) = S*VT( M, I )
! 483: VT( M, I ) = C*VT( M, I )
! 484: 180 CONTINUE
! 485: ELSE
! 486: CALL DCOPY( M, VT( NLP1, 1 ), LDVT, VT2( 1, 1 ), LDVT2 )
! 487: END IF
! 488: IF( M.GT.N ) THEN
! 489: CALL DCOPY( M, VT( M, 1 ), LDVT, VT2( M, 1 ), LDVT2 )
! 490: END IF
! 491: *
! 492: * The deflated singular values and their corresponding vectors go
! 493: * into the back of D, U, and V respectively.
! 494: *
! 495: IF( N.GT.K ) THEN
! 496: CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
! 497: CALL DLACPY( 'A', N, N-K, U2( 1, K+1 ), LDU2, U( 1, K+1 ),
! 498: $ LDU )
! 499: CALL DLACPY( 'A', N-K, M, VT2( K+1, 1 ), LDVT2, VT( K+1, 1 ),
! 500: $ LDVT )
! 501: END IF
! 502: *
! 503: * Copy CTOT into COLTYP for referencing in DLASD3.
! 504: *
! 505: DO 190 J = 1, 4
! 506: COLTYP( J ) = CTOT( J )
! 507: 190 CONTINUE
! 508: *
! 509: RETURN
! 510: *
! 511: * End of DLASD2
! 512: *
! 513: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack