Annotation of rpl/lapack/lapack/dlaed8.f, revision 1.4

1.1       bertrand    1:       SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO,
                      2:      $                   CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR,
                      3:      $                   GIVCOL, GIVNUM, INDXP, INDX, INFO )
                      4: *
                      5: *  -- LAPACK 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            CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N,
                     12:      $                   QSIZ
                     13:       DOUBLE PRECISION   RHO
                     14: *     ..
                     15: *     .. Array Arguments ..
                     16:       INTEGER            GIVCOL( 2, * ), INDX( * ), INDXP( * ),
                     17:      $                   INDXQ( * ), PERM( * )
                     18:       DOUBLE PRECISION   D( * ), DLAMDA( * ), GIVNUM( 2, * ),
                     19:      $                   Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
                     20: *     ..
                     21: *
                     22: *  Purpose
                     23: *  =======
                     24: *
                     25: *  DLAED8 merges the two sets of eigenvalues 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: *  eigenvalues are close together or if there is a tiny element in the
                     29: *  Z vector.  For each such occurrence the order of the related secular
                     30: *  equation problem is reduced by one.
                     31: *
                     32: *  Arguments
                     33: *  =========
                     34: *
                     35: *  ICOMPQ  (input) INTEGER
                     36: *          = 0:  Compute eigenvalues only.
                     37: *          = 1:  Compute eigenvectors of original dense symmetric matrix
                     38: *                also.  On entry, Q contains the orthogonal matrix used
                     39: *                to reduce the original matrix to tridiagonal form.
                     40: *
                     41: *  K      (output) INTEGER
                     42: *         The number of non-deflated eigenvalues, and the order of the
                     43: *         related secular equation.
                     44: *
                     45: *  N      (input) INTEGER
                     46: *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
                     47: *
                     48: *  QSIZ   (input) INTEGER
                     49: *         The dimension of the orthogonal matrix used to reduce
                     50: *         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
                     51: *
                     52: *  D      (input/output) DOUBLE PRECISION array, dimension (N)
                     53: *         On entry, the eigenvalues of the two submatrices to be
                     54: *         combined.  On exit, the trailing (N-K) updated eigenvalues
                     55: *         (those which were deflated) sorted into increasing order.
                     56: *
                     57: *  Q      (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
                     58: *         If ICOMPQ = 0, Q is not referenced.  Otherwise,
                     59: *         on entry, Q contains the eigenvectors of the partially solved
                     60: *         system which has been previously updated in matrix
                     61: *         multiplies with other partially solved eigensystems.
                     62: *         On exit, Q contains the trailing (N-K) updated eigenvectors
                     63: *         (those which were deflated) in its last N-K columns.
                     64: *
                     65: *  LDQ    (input) INTEGER
                     66: *         The leading dimension of the array Q.  LDQ >= max(1,N).
                     67: *
                     68: *  INDXQ  (input) INTEGER array, dimension (N)
                     69: *         The permutation which separately sorts the two sub-problems
                     70: *         in D into ascending order.  Note that elements in the second
                     71: *         half of this permutation must first have CUTPNT added to
                     72: *         their values in order to be accurate.
                     73: *
                     74: *  RHO    (input/output) DOUBLE PRECISION
                     75: *         On entry, the off-diagonal element associated with the rank-1
                     76: *         cut which originally split the two submatrices which are now
                     77: *         being recombined.
                     78: *         On exit, RHO has been modified to the value required by
                     79: *         DLAED3.
                     80: *
                     81: *  CUTPNT (input) INTEGER
                     82: *         The location of the last eigenvalue in the leading
                     83: *         sub-matrix.  min(1,N) <= CUTPNT <= N.
                     84: *
                     85: *  Z      (input) DOUBLE PRECISION array, dimension (N)
                     86: *         On entry, Z contains the updating vector (the last row of
                     87: *         the first sub-eigenvector matrix and the first row of the
                     88: *         second sub-eigenvector matrix).
                     89: *         On exit, the contents of Z are destroyed by the updating
                     90: *         process.
                     91: *
                     92: *  DLAMDA (output) DOUBLE PRECISION array, dimension (N)
                     93: *         A copy of the first K eigenvalues which will be used by
                     94: *         DLAED3 to form the secular equation.
                     95: *
                     96: *  Q2     (output) DOUBLE PRECISION array, dimension (LDQ2,N)
                     97: *         If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
                     98: *         a copy of the first K eigenvectors which will be used by
                     99: *         DLAED7 in a matrix multiply (DGEMM) to update the new
                    100: *         eigenvectors.
                    101: *
                    102: *  LDQ2   (input) INTEGER
                    103: *         The leading dimension of the array Q2.  LDQ2 >= max(1,N).
                    104: *
                    105: *  W      (output) DOUBLE PRECISION array, dimension (N)
                    106: *         The first k values of the final deflation-altered z-vector and
                    107: *         will be passed to DLAED3.
                    108: *
                    109: *  PERM   (output) INTEGER array, dimension (N)
                    110: *         The permutations (from deflation and sorting) to be applied
                    111: *         to each eigenblock.
                    112: *
                    113: *  GIVPTR (output) INTEGER
                    114: *         The number of Givens rotations which took place in this
                    115: *         subproblem.
                    116: *
                    117: *  GIVCOL (output) INTEGER array, dimension (2, N)
                    118: *         Each pair of numbers indicates a pair of columns to take place
                    119: *         in a Givens rotation.
                    120: *
                    121: *  GIVNUM (output) DOUBLE PRECISION array, dimension (2, N)
                    122: *         Each number indicates the S value to be used in the
                    123: *         corresponding Givens rotation.
                    124: *
                    125: *  INDXP  (workspace) INTEGER array, dimension (N)
                    126: *         The permutation used to place deflated values of D at the end
                    127: *         of the array.  INDXP(1:K) points to the nondeflated D-values
                    128: *         and INDXP(K+1:N) points to the deflated eigenvalues.
                    129: *
                    130: *  INDX   (workspace) INTEGER array, dimension (N)
                    131: *         The permutation used to sort the contents of D into ascending
                    132: *         order.
                    133: *
                    134: *  INFO   (output) INTEGER
                    135: *          = 0:  successful exit.
                    136: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
                    137: *
                    138: *  Further Details
                    139: *  ===============
                    140: *
                    141: *  Based on contributions by
                    142: *     Jeff Rutter, Computer Science Division, University of California
                    143: *     at Berkeley, USA
                    144: *
                    145: *  =====================================================================
                    146: *
                    147: *     .. Parameters ..
                    148:       DOUBLE PRECISION   MONE, ZERO, ONE, TWO, EIGHT
                    149:       PARAMETER          ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0,
                    150:      $                   TWO = 2.0D0, EIGHT = 8.0D0 )
                    151: *     ..
                    152: *     .. Local Scalars ..
                    153: *
                    154:       INTEGER            I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
                    155:       DOUBLE PRECISION   C, EPS, S, T, TAU, TOL
                    156: *     ..
                    157: *     .. External Functions ..
                    158:       INTEGER            IDAMAX
                    159:       DOUBLE PRECISION   DLAMCH, DLAPY2
                    160:       EXTERNAL           IDAMAX, DLAMCH, DLAPY2
                    161: *     ..
                    162: *     .. External Subroutines ..
                    163:       EXTERNAL           DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA
                    164: *     ..
                    165: *     .. Intrinsic Functions ..
                    166:       INTRINSIC          ABS, MAX, MIN, SQRT
                    167: *     ..
                    168: *     .. Executable Statements ..
                    169: *
                    170: *     Test the input parameters.
                    171: *
                    172:       INFO = 0
                    173: *
                    174:       IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
                    175:          INFO = -1
                    176:       ELSE IF( N.LT.0 ) THEN
                    177:          INFO = -3
                    178:       ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
                    179:          INFO = -4
                    180:       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
                    181:          INFO = -7
                    182:       ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN
                    183:          INFO = -10
                    184:       ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN
                    185:          INFO = -14
                    186:       END IF
                    187:       IF( INFO.NE.0 ) THEN
                    188:          CALL XERBLA( 'DLAED8', -INFO )
                    189:          RETURN
                    190:       END IF
                    191: *
                    192: *     Quick return if possible
                    193: *
                    194:       IF( N.EQ.0 )
                    195:      $   RETURN
                    196: *
                    197:       N1 = CUTPNT
                    198:       N2 = N - N1
                    199:       N1P1 = N1 + 1
                    200: *
                    201:       IF( RHO.LT.ZERO ) THEN
                    202:          CALL DSCAL( N2, MONE, Z( N1P1 ), 1 )
                    203:       END IF
                    204: *
                    205: *     Normalize z so that norm(z) = 1
                    206: *
                    207:       T = ONE / SQRT( TWO )
                    208:       DO 10 J = 1, N
                    209:          INDX( J ) = J
                    210:    10 CONTINUE
                    211:       CALL DSCAL( N, T, Z, 1 )
                    212:       RHO = ABS( TWO*RHO )
                    213: *
                    214: *     Sort the eigenvalues into increasing order
                    215: *
                    216:       DO 20 I = CUTPNT + 1, N
                    217:          INDXQ( I ) = INDXQ( I ) + CUTPNT
                    218:    20 CONTINUE
                    219:       DO 30 I = 1, N
                    220:          DLAMDA( I ) = D( INDXQ( I ) )
                    221:          W( I ) = Z( INDXQ( I ) )
                    222:    30 CONTINUE
                    223:       I = 1
                    224:       J = CUTPNT + 1
                    225:       CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
                    226:       DO 40 I = 1, N
                    227:          D( I ) = DLAMDA( INDX( I ) )
                    228:          Z( I ) = W( INDX( I ) )
                    229:    40 CONTINUE
                    230: *
                    231: *     Calculate the allowable deflation tolerence
                    232: *
                    233:       IMAX = IDAMAX( N, Z, 1 )
                    234:       JMAX = IDAMAX( N, D, 1 )
                    235:       EPS = DLAMCH( 'Epsilon' )
                    236:       TOL = EIGHT*EPS*ABS( D( JMAX ) )
                    237: *
                    238: *     If the rank-1 modifier is small enough, no more needs to be done
                    239: *     except to reorganize Q so that its columns correspond with the
                    240: *     elements in D.
                    241: *
                    242:       IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
                    243:          K = 0
                    244:          IF( ICOMPQ.EQ.0 ) THEN
                    245:             DO 50 J = 1, N
                    246:                PERM( J ) = INDXQ( INDX( J ) )
                    247:    50       CONTINUE
                    248:          ELSE
                    249:             DO 60 J = 1, N
                    250:                PERM( J ) = INDXQ( INDX( J ) )
                    251:                CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
                    252:    60       CONTINUE
                    253:             CALL DLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ),
                    254:      $                   LDQ )
                    255:          END IF
                    256:          RETURN
                    257:       END IF
                    258: *
                    259: *     If there are multiple eigenvalues then the problem deflates.  Here
                    260: *     the number of equal eigenvalues are found.  As each equal
                    261: *     eigenvalue is found, an elementary reflector is computed to rotate
                    262: *     the corresponding eigensubspace so that the corresponding
                    263: *     components of Z are zero in this new basis.
                    264: *
                    265:       K = 0
                    266:       GIVPTR = 0
                    267:       K2 = N + 1
                    268:       DO 70 J = 1, N
                    269:          IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
                    270: *
                    271: *           Deflate due to small z component.
                    272: *
                    273:             K2 = K2 - 1
                    274:             INDXP( K2 ) = J
                    275:             IF( J.EQ.N )
                    276:      $         GO TO 110
                    277:          ELSE
                    278:             JLAM = J
                    279:             GO TO 80
                    280:          END IF
                    281:    70 CONTINUE
                    282:    80 CONTINUE
                    283:       J = J + 1
                    284:       IF( J.GT.N )
                    285:      $   GO TO 100
                    286:       IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
                    287: *
                    288: *        Deflate due to small z component.
                    289: *
                    290:          K2 = K2 - 1
                    291:          INDXP( K2 ) = J
                    292:       ELSE
                    293: *
                    294: *        Check if eigenvalues are close enough to allow deflation.
                    295: *
                    296:          S = Z( JLAM )
                    297:          C = Z( J )
                    298: *
                    299: *        Find sqrt(a**2+b**2) without overflow or
                    300: *        destructive underflow.
                    301: *
                    302:          TAU = DLAPY2( C, S )
                    303:          T = D( J ) - D( JLAM )
                    304:          C = C / TAU
                    305:          S = -S / TAU
                    306:          IF( ABS( T*C*S ).LE.TOL ) THEN
                    307: *
                    308: *           Deflation is possible.
                    309: *
                    310:             Z( J ) = TAU
                    311:             Z( JLAM ) = ZERO
                    312: *
                    313: *           Record the appropriate Givens rotation
                    314: *
                    315:             GIVPTR = GIVPTR + 1
                    316:             GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
                    317:             GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
                    318:             GIVNUM( 1, GIVPTR ) = C
                    319:             GIVNUM( 2, GIVPTR ) = S
                    320:             IF( ICOMPQ.EQ.1 ) THEN
                    321:                CALL DROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
                    322:      $                    Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
                    323:             END IF
                    324:             T = D( JLAM )*C*C + D( J )*S*S
                    325:             D( J ) = D( JLAM )*S*S + D( J )*C*C
                    326:             D( JLAM ) = T
                    327:             K2 = K2 - 1
                    328:             I = 1
                    329:    90       CONTINUE
                    330:             IF( K2+I.LE.N ) THEN
                    331:                IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
                    332:                   INDXP( K2+I-1 ) = INDXP( K2+I )
                    333:                   INDXP( K2+I ) = JLAM
                    334:                   I = I + 1
                    335:                   GO TO 90
                    336:                ELSE
                    337:                   INDXP( K2+I-1 ) = JLAM
                    338:                END IF
                    339:             ELSE
                    340:                INDXP( K2+I-1 ) = JLAM
                    341:             END IF
                    342:             JLAM = J
                    343:          ELSE
                    344:             K = K + 1
                    345:             W( K ) = Z( JLAM )
                    346:             DLAMDA( K ) = D( JLAM )
                    347:             INDXP( K ) = JLAM
                    348:             JLAM = J
                    349:          END IF
                    350:       END IF
                    351:       GO TO 80
                    352:   100 CONTINUE
                    353: *
                    354: *     Record the last eigenvalue.
                    355: *
                    356:       K = K + 1
                    357:       W( K ) = Z( JLAM )
                    358:       DLAMDA( K ) = D( JLAM )
                    359:       INDXP( K ) = JLAM
                    360: *
                    361:   110 CONTINUE
                    362: *
                    363: *     Sort the eigenvalues and corresponding eigenvectors into DLAMDA
                    364: *     and Q2 respectively.  The eigenvalues/vectors which were not
                    365: *     deflated go into the first K slots of DLAMDA and Q2 respectively,
                    366: *     while those which were deflated go into the last N - K slots.
                    367: *
                    368:       IF( ICOMPQ.EQ.0 ) THEN
                    369:          DO 120 J = 1, N
                    370:             JP = INDXP( J )
                    371:             DLAMDA( J ) = D( JP )
                    372:             PERM( J ) = INDXQ( INDX( JP ) )
                    373:   120    CONTINUE
                    374:       ELSE
                    375:          DO 130 J = 1, N
                    376:             JP = INDXP( J )
                    377:             DLAMDA( J ) = D( JP )
                    378:             PERM( J ) = INDXQ( INDX( JP ) )
                    379:             CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
                    380:   130    CONTINUE
                    381:       END IF
                    382: *
                    383: *     The deflated eigenvalues and their corresponding vectors go back
                    384: *     into the last N - K slots of D and Q respectively.
                    385: *
                    386:       IF( K.LT.N ) THEN
                    387:          IF( ICOMPQ.EQ.0 ) THEN
                    388:             CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
                    389:          ELSE
                    390:             CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
                    391:             CALL DLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2,
                    392:      $                   Q( 1, K+1 ), LDQ )
                    393:          END IF
                    394:       END IF
                    395: *
                    396:       RETURN
                    397: *
                    398: *     End of DLAED8
                    399: *
                    400:       END

CVSweb interface <joel.bertrand@systella.fr>