Annotation of rpl/lapack/lapack/dlasd7.f, revision 1.6

1.1       bertrand    1:       SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
                      2:      $                   VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
                      3:      $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
                      4:      $                   C, S, INFO )
                      5: *
                      6: *  -- LAPACK auxiliary routine (version 3.2) --
                      7: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      8: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                      9: *     November 2006
                     10: *
                     11: *     .. Scalar Arguments ..
                     12:       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
                     13:      $                   NR, SQRE
                     14:       DOUBLE PRECISION   ALPHA, BETA, C, S
                     15: *     ..
                     16: *     .. Array Arguments ..
                     17:       INTEGER            GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
                     18:      $                   IDXQ( * ), PERM( * )
                     19:       DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
                     20:      $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
                     21:      $                   ZW( * )
                     22: *     ..
                     23: *
                     24: *  Purpose
                     25: *  =======
                     26: *
                     27: *  DLASD7 merges the two sets of singular values together into a single
                     28: *  sorted set. Then it tries to deflate the size of the problem. There
                     29: *  are two ways in which deflation can occur:  when two or more singular
                     30: *  values are close together or if there is a tiny entry in the Z
                     31: *  vector. For each such occurrence the order of the related
                     32: *  secular equation problem is reduced by one.
                     33: *
                     34: *  DLASD7 is called from DLASD6.
                     35: *
                     36: *  Arguments
                     37: *  =========
                     38: *
                     39: *  ICOMPQ  (input) INTEGER
                     40: *          Specifies whether singular vectors are to be computed
                     41: *          in compact form, as follows:
                     42: *          = 0: Compute singular values only.
                     43: *          = 1: Compute singular vectors of upper
                     44: *               bidiagonal matrix in compact form.
                     45: *
                     46: *  NL     (input) INTEGER
                     47: *         The row dimension of the upper block. NL >= 1.
                     48: *
                     49: *  NR     (input) INTEGER
                     50: *         The row dimension of the lower block. NR >= 1.
                     51: *
                     52: *  SQRE   (input) INTEGER
                     53: *         = 0: the lower block is an NR-by-NR square matrix.
                     54: *         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
                     55: *
                     56: *         The bidiagonal matrix has
                     57: *         N = NL + NR + 1 rows and
                     58: *         M = N + SQRE >= N columns.
                     59: *
                     60: *  K      (output) INTEGER
                     61: *         Contains the dimension of the non-deflated matrix, this is
                     62: *         the order of the related secular equation. 1 <= K <=N.
                     63: *
                     64: *  D      (input/output) DOUBLE PRECISION array, dimension ( N )
                     65: *         On entry D contains the singular values of the two submatrices
                     66: *         to be combined. On exit D contains the trailing (N-K) updated
                     67: *         singular values (those which were deflated) sorted into
                     68: *         increasing order.
                     69: *
                     70: *  Z      (output) DOUBLE PRECISION array, dimension ( M )
                     71: *         On exit Z contains the updating row vector in the secular
                     72: *         equation.
                     73: *
                     74: *  ZW     (workspace) DOUBLE PRECISION array, dimension ( M )
                     75: *         Workspace for Z.
                     76: *
                     77: *  VF     (input/output) DOUBLE PRECISION array, dimension ( M )
                     78: *         On entry, VF(1:NL+1) contains the first components of all
                     79: *         right singular vectors of the upper block; and VF(NL+2:M)
                     80: *         contains the first components of all right singular vectors
                     81: *         of the lower block. On exit, VF contains the first components
                     82: *         of all right singular vectors of the bidiagonal matrix.
                     83: *
                     84: *  VFW    (workspace) DOUBLE PRECISION array, dimension ( M )
                     85: *         Workspace for VF.
                     86: *
                     87: *  VL     (input/output) DOUBLE PRECISION array, dimension ( M )
                     88: *         On entry, VL(1:NL+1) contains the  last components of all
                     89: *         right singular vectors of the upper block; and VL(NL+2:M)
                     90: *         contains the last components of all right singular vectors
                     91: *         of the lower block. On exit, VL contains the last components
                     92: *         of all right singular vectors of the bidiagonal matrix.
                     93: *
                     94: *  VLW    (workspace) DOUBLE PRECISION array, dimension ( M )
                     95: *         Workspace for VL.
                     96: *
                     97: *  ALPHA  (input) DOUBLE PRECISION
                     98: *         Contains the diagonal element associated with the added row.
                     99: *
                    100: *  BETA   (input) DOUBLE PRECISION
                    101: *         Contains the off-diagonal element associated with the added
                    102: *         row.
                    103: *
                    104: *  DSIGMA (output) DOUBLE PRECISION array, dimension ( N )
                    105: *         Contains a copy of the diagonal elements (K-1 singular values
                    106: *         and one zero) in the secular equation.
                    107: *
                    108: *  IDX    (workspace) INTEGER array, dimension ( N )
                    109: *         This will contain the permutation used to sort the contents of
                    110: *         D into ascending order.
                    111: *
                    112: *  IDXP   (workspace) INTEGER array, dimension ( N )
                    113: *         This will contain the permutation used to place deflated
                    114: *         values of D at the end of the array. On output IDXP(2:K)
                    115: *         points to the nondeflated D-values and IDXP(K+1:N)
                    116: *         points to the deflated singular values.
                    117: *
                    118: *  IDXQ   (input) INTEGER array, dimension ( N )
                    119: *         This contains the permutation which separately sorts the two
                    120: *         sub-problems in D into ascending order.  Note that entries in
                    121: *         the first half of this permutation must first be moved one
                    122: *         position backward; and entries in the second half
                    123: *         must first have NL+1 added to their values.
                    124: *
                    125: *  PERM   (output) INTEGER array, dimension ( N )
                    126: *         The permutations (from deflation and sorting) to be applied
                    127: *         to each singular block. Not referenced if ICOMPQ = 0.
                    128: *
                    129: *  GIVPTR (output) INTEGER
                    130: *         The number of Givens rotations which took place in this
                    131: *         subproblem. Not referenced if ICOMPQ = 0.
                    132: *
                    133: *  GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 )
                    134: *         Each pair of numbers indicates a pair of columns to take place
                    135: *         in a Givens rotation. Not referenced if ICOMPQ = 0.
                    136: *
                    137: *  LDGCOL (input) INTEGER
                    138: *         The leading dimension of GIVCOL, must be at least N.
                    139: *
                    140: *  GIVNUM (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
                    141: *         Each number indicates the C or S value to be used in the
                    142: *         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
                    143: *
                    144: *  LDGNUM (input) INTEGER
                    145: *         The leading dimension of GIVNUM, must be at least N.
                    146: *
                    147: *  C      (output) DOUBLE PRECISION
                    148: *         C contains garbage if SQRE =0 and the C-value of a Givens
                    149: *         rotation related to the right null space if SQRE = 1.
                    150: *
                    151: *  S      (output) DOUBLE PRECISION
                    152: *         S contains garbage if SQRE =0 and the S-value of a Givens
                    153: *         rotation related to the right null space if SQRE = 1.
                    154: *
                    155: *  INFO   (output) INTEGER
                    156: *         = 0:  successful exit.
                    157: *         < 0:  if INFO = -i, the i-th argument had an illegal value.
                    158: *
                    159: *  Further Details
                    160: *  ===============
                    161: *
                    162: *  Based on contributions by
                    163: *     Ming Gu and Huan Ren, Computer Science Division, University of
                    164: *     California at Berkeley, USA
                    165: *
                    166: *  =====================================================================
                    167: *
                    168: *     .. Parameters ..
                    169:       DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT
                    170:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
                    171:      $                   EIGHT = 8.0D+0 )
                    172: *     ..
                    173: *     .. Local Scalars ..
                    174: *
                    175:       INTEGER            I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
                    176:      $                   NLP1, NLP2
                    177:       DOUBLE PRECISION   EPS, HLFTOL, TAU, TOL, Z1
                    178: *     ..
                    179: *     .. External Subroutines ..
                    180:       EXTERNAL           DCOPY, DLAMRG, DROT, XERBLA
                    181: *     ..
                    182: *     .. External Functions ..
                    183:       DOUBLE PRECISION   DLAMCH, DLAPY2
                    184:       EXTERNAL           DLAMCH, DLAPY2
                    185: *     ..
                    186: *     .. Intrinsic Functions ..
                    187:       INTRINSIC          ABS, MAX
                    188: *     ..
                    189: *     .. Executable Statements ..
                    190: *
                    191: *     Test the input parameters.
                    192: *
                    193:       INFO = 0
                    194:       N = NL + NR + 1
                    195:       M = N + SQRE
                    196: *
                    197:       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
                    198:          INFO = -1
                    199:       ELSE IF( NL.LT.1 ) THEN
                    200:          INFO = -2
                    201:       ELSE IF( NR.LT.1 ) THEN
                    202:          INFO = -3
                    203:       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
                    204:          INFO = -4
                    205:       ELSE IF( LDGCOL.LT.N ) THEN
                    206:          INFO = -22
                    207:       ELSE IF( LDGNUM.LT.N ) THEN
                    208:          INFO = -24
                    209:       END IF
                    210:       IF( INFO.NE.0 ) THEN
                    211:          CALL XERBLA( 'DLASD7', -INFO )
                    212:          RETURN
                    213:       END IF
                    214: *
                    215:       NLP1 = NL + 1
                    216:       NLP2 = NL + 2
                    217:       IF( ICOMPQ.EQ.1 ) THEN
                    218:          GIVPTR = 0
                    219:       END IF
                    220: *
                    221: *     Generate the first part of the vector Z and move the singular
                    222: *     values in the first part of D one position backward.
                    223: *
                    224:       Z1 = ALPHA*VL( NLP1 )
                    225:       VL( NLP1 ) = ZERO
                    226:       TAU = VF( NLP1 )
                    227:       DO 10 I = NL, 1, -1
                    228:          Z( I+1 ) = ALPHA*VL( I )
                    229:          VL( I ) = ZERO
                    230:          VF( I+1 ) = VF( I )
                    231:          D( I+1 ) = D( I )
                    232:          IDXQ( I+1 ) = IDXQ( I ) + 1
                    233:    10 CONTINUE
                    234:       VF( 1 ) = TAU
                    235: *
                    236: *     Generate the second part of the vector Z.
                    237: *
                    238:       DO 20 I = NLP2, M
                    239:          Z( I ) = BETA*VF( I )
                    240:          VF( I ) = ZERO
                    241:    20 CONTINUE
                    242: *
                    243: *     Sort the singular values into increasing order
                    244: *
                    245:       DO 30 I = NLP2, N
                    246:          IDXQ( I ) = IDXQ( I ) + NLP1
                    247:    30 CONTINUE
                    248: *
                    249: *     DSIGMA, IDXC, IDXC, and ZW are used as storage space.
                    250: *
                    251:       DO 40 I = 2, N
                    252:          DSIGMA( I ) = D( IDXQ( I ) )
                    253:          ZW( I ) = Z( IDXQ( I ) )
                    254:          VFW( I ) = VF( IDXQ( I ) )
                    255:          VLW( I ) = VL( IDXQ( I ) )
                    256:    40 CONTINUE
                    257: *
                    258:       CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
                    259: *
                    260:       DO 50 I = 2, N
                    261:          IDXI = 1 + IDX( I )
                    262:          D( I ) = DSIGMA( IDXI )
                    263:          Z( I ) = ZW( IDXI )
                    264:          VF( I ) = VFW( IDXI )
                    265:          VL( I ) = VLW( IDXI )
                    266:    50 CONTINUE
                    267: *
                    268: *     Calculate the allowable deflation tolerence
                    269: *
                    270:       EPS = DLAMCH( 'Epsilon' )
                    271:       TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
                    272:       TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
                    273: *
                    274: *     There are 2 kinds of deflation -- first a value in the z-vector
                    275: *     is small, second two (or more) singular values are very close
                    276: *     together (their difference is small).
                    277: *
                    278: *     If the value in the z-vector is small, we simply permute the
                    279: *     array so that the corresponding singular value is moved to the
                    280: *     end.
                    281: *
                    282: *     If two values in the D-vector are close, we perform a two-sided
                    283: *     rotation designed to make one of the corresponding z-vector
                    284: *     entries zero, and then permute the array so that the deflated
                    285: *     singular value is moved to the end.
                    286: *
                    287: *     If there are multiple singular values then the problem deflates.
                    288: *     Here the number of equal singular values are found.  As each equal
                    289: *     singular value is found, an elementary reflector is computed to
                    290: *     rotate the corresponding singular subspace so that the
                    291: *     corresponding components of Z are zero in this new basis.
                    292: *
                    293:       K = 1
                    294:       K2 = N + 1
                    295:       DO 60 J = 2, N
                    296:          IF( ABS( Z( J ) ).LE.TOL ) THEN
                    297: *
                    298: *           Deflate due to small z component.
                    299: *
                    300:             K2 = K2 - 1
                    301:             IDXP( K2 ) = J
                    302:             IF( J.EQ.N )
                    303:      $         GO TO 100
                    304:          ELSE
                    305:             JPREV = J
                    306:             GO TO 70
                    307:          END IF
                    308:    60 CONTINUE
                    309:    70 CONTINUE
                    310:       J = JPREV
                    311:    80 CONTINUE
                    312:       J = J + 1
                    313:       IF( J.GT.N )
                    314:      $   GO TO 90
                    315:       IF( ABS( Z( J ) ).LE.TOL ) THEN
                    316: *
                    317: *        Deflate due to small z component.
                    318: *
                    319:          K2 = K2 - 1
                    320:          IDXP( K2 ) = J
                    321:       ELSE
                    322: *
                    323: *        Check if singular values are close enough to allow deflation.
                    324: *
                    325:          IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
                    326: *
                    327: *           Deflation is possible.
                    328: *
                    329:             S = Z( JPREV )
                    330:             C = Z( J )
                    331: *
                    332: *           Find sqrt(a**2+b**2) without overflow or
                    333: *           destructive underflow.
                    334: *
                    335:             TAU = DLAPY2( C, S )
                    336:             Z( J ) = TAU
                    337:             Z( JPREV ) = ZERO
                    338:             C = C / TAU
                    339:             S = -S / TAU
                    340: *
                    341: *           Record the appropriate Givens rotation
                    342: *
                    343:             IF( ICOMPQ.EQ.1 ) THEN
                    344:                GIVPTR = GIVPTR + 1
                    345:                IDXJP = IDXQ( IDX( JPREV )+1 )
                    346:                IDXJ = IDXQ( IDX( J )+1 )
                    347:                IF( IDXJP.LE.NLP1 ) THEN
                    348:                   IDXJP = IDXJP - 1
                    349:                END IF
                    350:                IF( IDXJ.LE.NLP1 ) THEN
                    351:                   IDXJ = IDXJ - 1
                    352:                END IF
                    353:                GIVCOL( GIVPTR, 2 ) = IDXJP
                    354:                GIVCOL( GIVPTR, 1 ) = IDXJ
                    355:                GIVNUM( GIVPTR, 2 ) = C
                    356:                GIVNUM( GIVPTR, 1 ) = S
                    357:             END IF
                    358:             CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
                    359:             CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
                    360:             K2 = K2 - 1
                    361:             IDXP( K2 ) = JPREV
                    362:             JPREV = J
                    363:          ELSE
                    364:             K = K + 1
                    365:             ZW( K ) = Z( JPREV )
                    366:             DSIGMA( K ) = D( JPREV )
                    367:             IDXP( K ) = JPREV
                    368:             JPREV = J
                    369:          END IF
                    370:       END IF
                    371:       GO TO 80
                    372:    90 CONTINUE
                    373: *
                    374: *     Record the last singular value.
                    375: *
                    376:       K = K + 1
                    377:       ZW( K ) = Z( JPREV )
                    378:       DSIGMA( K ) = D( JPREV )
                    379:       IDXP( K ) = JPREV
                    380: *
                    381:   100 CONTINUE
                    382: *
                    383: *     Sort the singular values into DSIGMA. The singular values which
                    384: *     were not deflated go into the first K slots of DSIGMA, except
                    385: *     that DSIGMA(1) is treated separately.
                    386: *
                    387:       DO 110 J = 2, N
                    388:          JP = IDXP( J )
                    389:          DSIGMA( J ) = D( JP )
                    390:          VFW( J ) = VF( JP )
                    391:          VLW( J ) = VL( JP )
                    392:   110 CONTINUE
                    393:       IF( ICOMPQ.EQ.1 ) THEN
                    394:          DO 120 J = 2, N
                    395:             JP = IDXP( J )
                    396:             PERM( J ) = IDXQ( IDX( JP )+1 )
                    397:             IF( PERM( J ).LE.NLP1 ) THEN
                    398:                PERM( J ) = PERM( J ) - 1
                    399:             END IF
                    400:   120    CONTINUE
                    401:       END IF
                    402: *
                    403: *     The deflated singular values go back into the last N - K slots of
                    404: *     D.
                    405: *
                    406:       CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
                    407: *
                    408: *     Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
                    409: *     VL(M).
                    410: *
                    411:       DSIGMA( 1 ) = ZERO
                    412:       HLFTOL = TOL / TWO
                    413:       IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
                    414:      $   DSIGMA( 2 ) = HLFTOL
                    415:       IF( M.GT.N ) THEN
                    416:          Z( 1 ) = DLAPY2( Z1, Z( M ) )
                    417:          IF( Z( 1 ).LE.TOL ) THEN
                    418:             C = ONE
                    419:             S = ZERO
                    420:             Z( 1 ) = TOL
                    421:          ELSE
                    422:             C = Z1 / Z( 1 )
                    423:             S = -Z( M ) / Z( 1 )
                    424:          END IF
                    425:          CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
                    426:          CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
                    427:       ELSE
                    428:          IF( ABS( Z1 ).LE.TOL ) THEN
                    429:             Z( 1 ) = TOL
                    430:          ELSE
                    431:             Z( 1 ) = Z1
                    432:          END IF
                    433:       END IF
                    434: *
                    435: *     Restore Z, VF, and VL.
                    436: *
                    437:       CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
                    438:       CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
                    439:       CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
                    440: *
                    441:       RETURN
                    442: *
                    443: *     End of DLASD7
                    444: *
                    445:       END

CVSweb interface <joel.bertrand@systella.fr>