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

1.12    ! bertrand    1: *> \brief \b DLAED8 used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense.
1.9       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
                      5: * Online html documentation available at 
                      6: *            http://www.netlib.org/lapack/explore-html/ 
                      7: *
                      8: *> \htmlonly
                      9: *> Download DLAED8 + dependencies 
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed8.f"> 
                     11: *> [TGZ]</a> 
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed8.f"> 
                     13: *> [ZIP]</a> 
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed8.f"> 
                     15: *> [TXT]</a>
                     16: *> \endhtmlonly 
                     17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO,
                     22: *                          CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR,
                     23: *                          GIVCOL, GIVNUM, INDXP, INDX, INFO )
                     24: * 
                     25: *       .. Scalar Arguments ..
                     26: *       INTEGER            CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N,
                     27: *      $                   QSIZ
                     28: *       DOUBLE PRECISION   RHO
                     29: *       ..
                     30: *       .. Array Arguments ..
                     31: *       INTEGER            GIVCOL( 2, * ), INDX( * ), INDXP( * ),
                     32: *      $                   INDXQ( * ), PERM( * )
                     33: *       DOUBLE PRECISION   D( * ), DLAMDA( * ), GIVNUM( 2, * ),
                     34: *      $                   Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
                     35: *       ..
                     36: *  
                     37: *
                     38: *> \par Purpose:
                     39: *  =============
                     40: *>
                     41: *> \verbatim
                     42: *>
                     43: *> DLAED8 merges the two sets of eigenvalues together into a single
                     44: *> sorted set.  Then it tries to deflate the size of the problem.
                     45: *> There are two ways in which deflation can occur:  when two or more
                     46: *> eigenvalues are close together or if there is a tiny element in the
                     47: *> Z vector.  For each such occurrence the order of the related secular
                     48: *> equation problem is reduced by one.
                     49: *> \endverbatim
                     50: *
                     51: *  Arguments:
                     52: *  ==========
                     53: *
                     54: *> \param[in] ICOMPQ
                     55: *> \verbatim
                     56: *>          ICOMPQ is INTEGER
                     57: *>          = 0:  Compute eigenvalues only.
                     58: *>          = 1:  Compute eigenvectors of original dense symmetric matrix
                     59: *>                also.  On entry, Q contains the orthogonal matrix used
                     60: *>                to reduce the original matrix to tridiagonal form.
                     61: *> \endverbatim
                     62: *>
                     63: *> \param[out] K
                     64: *> \verbatim
                     65: *>          K is INTEGER
                     66: *>         The number of non-deflated eigenvalues, and the order of the
                     67: *>         related secular equation.
                     68: *> \endverbatim
                     69: *>
                     70: *> \param[in] N
                     71: *> \verbatim
                     72: *>          N is INTEGER
                     73: *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
                     74: *> \endverbatim
                     75: *>
                     76: *> \param[in] QSIZ
                     77: *> \verbatim
                     78: *>          QSIZ is INTEGER
                     79: *>         The dimension of the orthogonal matrix used to reduce
                     80: *>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
                     81: *> \endverbatim
                     82: *>
                     83: *> \param[in,out] D
                     84: *> \verbatim
                     85: *>          D is DOUBLE PRECISION array, dimension (N)
                     86: *>         On entry, the eigenvalues of the two submatrices to be
                     87: *>         combined.  On exit, the trailing (N-K) updated eigenvalues
                     88: *>         (those which were deflated) sorted into increasing order.
                     89: *> \endverbatim
                     90: *>
                     91: *> \param[in,out] Q
                     92: *> \verbatim
                     93: *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
                     94: *>         If ICOMPQ = 0, Q is not referenced.  Otherwise,
                     95: *>         on entry, Q contains the eigenvectors of the partially solved
                     96: *>         system which has been previously updated in matrix
                     97: *>         multiplies with other partially solved eigensystems.
                     98: *>         On exit, Q contains the trailing (N-K) updated eigenvectors
                     99: *>         (those which were deflated) in its last N-K columns.
                    100: *> \endverbatim
                    101: *>
                    102: *> \param[in] LDQ
                    103: *> \verbatim
                    104: *>          LDQ is INTEGER
                    105: *>         The leading dimension of the array Q.  LDQ >= max(1,N).
                    106: *> \endverbatim
                    107: *>
                    108: *> \param[in] INDXQ
                    109: *> \verbatim
                    110: *>          INDXQ is INTEGER array, dimension (N)
                    111: *>         The permutation which separately sorts the two sub-problems
                    112: *>         in D into ascending order.  Note that elements in the second
                    113: *>         half of this permutation must first have CUTPNT added to
                    114: *>         their values in order to be accurate.
                    115: *> \endverbatim
                    116: *>
                    117: *> \param[in,out] RHO
                    118: *> \verbatim
                    119: *>          RHO is DOUBLE PRECISION
                    120: *>         On entry, the off-diagonal element associated with the rank-1
                    121: *>         cut which originally split the two submatrices which are now
                    122: *>         being recombined.
                    123: *>         On exit, RHO has been modified to the value required by
                    124: *>         DLAED3.
                    125: *> \endverbatim
                    126: *>
                    127: *> \param[in] CUTPNT
                    128: *> \verbatim
                    129: *>          CUTPNT is INTEGER
                    130: *>         The location of the last eigenvalue in the leading
                    131: *>         sub-matrix.  min(1,N) <= CUTPNT <= N.
                    132: *> \endverbatim
                    133: *>
                    134: *> \param[in] Z
                    135: *> \verbatim
                    136: *>          Z is DOUBLE PRECISION array, dimension (N)
                    137: *>         On entry, Z contains the updating vector (the last row of
                    138: *>         the first sub-eigenvector matrix and the first row of the
                    139: *>         second sub-eigenvector matrix).
                    140: *>         On exit, the contents of Z are destroyed by the updating
                    141: *>         process.
                    142: *> \endverbatim
                    143: *>
                    144: *> \param[out] DLAMDA
                    145: *> \verbatim
                    146: *>          DLAMDA is DOUBLE PRECISION array, dimension (N)
                    147: *>         A copy of the first K eigenvalues which will be used by
                    148: *>         DLAED3 to form the secular equation.
                    149: *> \endverbatim
                    150: *>
                    151: *> \param[out] Q2
                    152: *> \verbatim
                    153: *>          Q2 is DOUBLE PRECISION array, dimension (LDQ2,N)
                    154: *>         If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
                    155: *>         a copy of the first K eigenvectors which will be used by
                    156: *>         DLAED7 in a matrix multiply (DGEMM) to update the new
                    157: *>         eigenvectors.
                    158: *> \endverbatim
                    159: *>
                    160: *> \param[in] LDQ2
                    161: *> \verbatim
                    162: *>          LDQ2 is INTEGER
                    163: *>         The leading dimension of the array Q2.  LDQ2 >= max(1,N).
                    164: *> \endverbatim
                    165: *>
                    166: *> \param[out] W
                    167: *> \verbatim
                    168: *>          W is DOUBLE PRECISION array, dimension (N)
                    169: *>         The first k values of the final deflation-altered z-vector and
                    170: *>         will be passed to DLAED3.
                    171: *> \endverbatim
                    172: *>
                    173: *> \param[out] PERM
                    174: *> \verbatim
                    175: *>          PERM is INTEGER array, dimension (N)
                    176: *>         The permutations (from deflation and sorting) to be applied
                    177: *>         to each eigenblock.
                    178: *> \endverbatim
                    179: *>
                    180: *> \param[out] GIVPTR
                    181: *> \verbatim
                    182: *>          GIVPTR is INTEGER
                    183: *>         The number of Givens rotations which took place in this
                    184: *>         subproblem.
                    185: *> \endverbatim
                    186: *>
                    187: *> \param[out] GIVCOL
                    188: *> \verbatim
                    189: *>          GIVCOL is INTEGER array, dimension (2, N)
                    190: *>         Each pair of numbers indicates a pair of columns to take place
                    191: *>         in a Givens rotation.
                    192: *> \endverbatim
                    193: *>
                    194: *> \param[out] GIVNUM
                    195: *> \verbatim
                    196: *>          GIVNUM is DOUBLE PRECISION array, dimension (2, N)
                    197: *>         Each number indicates the S value to be used in the
                    198: *>         corresponding Givens rotation.
                    199: *> \endverbatim
                    200: *>
                    201: *> \param[out] INDXP
                    202: *> \verbatim
                    203: *>          INDXP is INTEGER array, dimension (N)
                    204: *>         The permutation used to place deflated values of D at the end
                    205: *>         of the array.  INDXP(1:K) points to the nondeflated D-values
                    206: *>         and INDXP(K+1:N) points to the deflated eigenvalues.
                    207: *> \endverbatim
                    208: *>
                    209: *> \param[out] INDX
                    210: *> \verbatim
                    211: *>          INDX is INTEGER array, dimension (N)
                    212: *>         The permutation used to sort the contents of D into ascending
                    213: *>         order.
                    214: *> \endverbatim
                    215: *>
                    216: *> \param[out] INFO
                    217: *> \verbatim
                    218: *>          INFO is INTEGER
                    219: *>          = 0:  successful exit.
                    220: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
                    221: *> \endverbatim
                    222: *
                    223: *  Authors:
                    224: *  ========
                    225: *
                    226: *> \author Univ. of Tennessee 
                    227: *> \author Univ. of California Berkeley 
                    228: *> \author Univ. of Colorado Denver 
                    229: *> \author NAG Ltd. 
                    230: *
1.12    ! bertrand  231: *> \date September 2012
1.9       bertrand  232: *
                    233: *> \ingroup auxOTHERcomputational
                    234: *
                    235: *> \par Contributors:
                    236: *  ==================
                    237: *>
                    238: *> Jeff Rutter, Computer Science Division, University of California
                    239: *> at Berkeley, USA
                    240: *
                    241: *  =====================================================================
1.1       bertrand  242:       SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO,
                    243:      $                   CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR,
                    244:      $                   GIVCOL, GIVNUM, INDXP, INDX, INFO )
                    245: *
1.12    ! bertrand  246: *  -- LAPACK computational routine (version 3.4.2) --
1.1       bertrand  247: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    248: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.12    ! bertrand  249: *     September 2012
1.1       bertrand  250: *
                    251: *     .. Scalar Arguments ..
                    252:       INTEGER            CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N,
                    253:      $                   QSIZ
                    254:       DOUBLE PRECISION   RHO
                    255: *     ..
                    256: *     .. Array Arguments ..
                    257:       INTEGER            GIVCOL( 2, * ), INDX( * ), INDXP( * ),
                    258:      $                   INDXQ( * ), PERM( * )
                    259:       DOUBLE PRECISION   D( * ), DLAMDA( * ), GIVNUM( 2, * ),
                    260:      $                   Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
                    261: *     ..
                    262: *
                    263: *  =====================================================================
                    264: *
                    265: *     .. Parameters ..
                    266:       DOUBLE PRECISION   MONE, ZERO, ONE, TWO, EIGHT
                    267:       PARAMETER          ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0,
                    268:      $                   TWO = 2.0D0, EIGHT = 8.0D0 )
                    269: *     ..
                    270: *     .. Local Scalars ..
                    271: *
                    272:       INTEGER            I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
                    273:       DOUBLE PRECISION   C, EPS, S, T, TAU, TOL
                    274: *     ..
                    275: *     .. External Functions ..
                    276:       INTEGER            IDAMAX
                    277:       DOUBLE PRECISION   DLAMCH, DLAPY2
                    278:       EXTERNAL           IDAMAX, DLAMCH, DLAPY2
                    279: *     ..
                    280: *     .. External Subroutines ..
                    281:       EXTERNAL           DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA
                    282: *     ..
                    283: *     .. Intrinsic Functions ..
                    284:       INTRINSIC          ABS, MAX, MIN, SQRT
                    285: *     ..
                    286: *     .. Executable Statements ..
                    287: *
                    288: *     Test the input parameters.
                    289: *
                    290:       INFO = 0
                    291: *
                    292:       IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
                    293:          INFO = -1
                    294:       ELSE IF( N.LT.0 ) THEN
                    295:          INFO = -3
                    296:       ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
                    297:          INFO = -4
                    298:       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
                    299:          INFO = -7
                    300:       ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN
                    301:          INFO = -10
                    302:       ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN
                    303:          INFO = -14
                    304:       END IF
                    305:       IF( INFO.NE.0 ) THEN
                    306:          CALL XERBLA( 'DLAED8', -INFO )
                    307:          RETURN
                    308:       END IF
                    309: *
1.5       bertrand  310: *     Need to initialize GIVPTR to O here in case of quick exit
                    311: *     to prevent an unspecified code behavior (usually sigfault) 
                    312: *     when IWORK array on entry to *stedc is not zeroed 
                    313: *     (or at least some IWORK entries which used in *laed7 for GIVPTR).
                    314: *
                    315:       GIVPTR = 0
                    316: *
1.1       bertrand  317: *     Quick return if possible
                    318: *
                    319:       IF( N.EQ.0 )
                    320:      $   RETURN
                    321: *
                    322:       N1 = CUTPNT
                    323:       N2 = N - N1
                    324:       N1P1 = N1 + 1
                    325: *
                    326:       IF( RHO.LT.ZERO ) THEN
                    327:          CALL DSCAL( N2, MONE, Z( N1P1 ), 1 )
                    328:       END IF
                    329: *
                    330: *     Normalize z so that norm(z) = 1
                    331: *
                    332:       T = ONE / SQRT( TWO )
                    333:       DO 10 J = 1, N
                    334:          INDX( J ) = J
                    335:    10 CONTINUE
                    336:       CALL DSCAL( N, T, Z, 1 )
                    337:       RHO = ABS( TWO*RHO )
                    338: *
                    339: *     Sort the eigenvalues into increasing order
                    340: *
                    341:       DO 20 I = CUTPNT + 1, N
                    342:          INDXQ( I ) = INDXQ( I ) + CUTPNT
                    343:    20 CONTINUE
                    344:       DO 30 I = 1, N
                    345:          DLAMDA( I ) = D( INDXQ( I ) )
                    346:          W( I ) = Z( INDXQ( I ) )
                    347:    30 CONTINUE
                    348:       I = 1
                    349:       J = CUTPNT + 1
                    350:       CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
                    351:       DO 40 I = 1, N
                    352:          D( I ) = DLAMDA( INDX( I ) )
                    353:          Z( I ) = W( INDX( I ) )
                    354:    40 CONTINUE
                    355: *
                    356: *     Calculate the allowable deflation tolerence
                    357: *
                    358:       IMAX = IDAMAX( N, Z, 1 )
                    359:       JMAX = IDAMAX( N, D, 1 )
                    360:       EPS = DLAMCH( 'Epsilon' )
                    361:       TOL = EIGHT*EPS*ABS( D( JMAX ) )
                    362: *
                    363: *     If the rank-1 modifier is small enough, no more needs to be done
                    364: *     except to reorganize Q so that its columns correspond with the
                    365: *     elements in D.
                    366: *
                    367:       IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
                    368:          K = 0
                    369:          IF( ICOMPQ.EQ.0 ) THEN
                    370:             DO 50 J = 1, N
                    371:                PERM( J ) = INDXQ( INDX( J ) )
                    372:    50       CONTINUE
                    373:          ELSE
                    374:             DO 60 J = 1, N
                    375:                PERM( J ) = INDXQ( INDX( J ) )
                    376:                CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
                    377:    60       CONTINUE
                    378:             CALL DLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ),
                    379:      $                   LDQ )
                    380:          END IF
                    381:          RETURN
                    382:       END IF
                    383: *
                    384: *     If there are multiple eigenvalues then the problem deflates.  Here
                    385: *     the number of equal eigenvalues are found.  As each equal
                    386: *     eigenvalue is found, an elementary reflector is computed to rotate
                    387: *     the corresponding eigensubspace so that the corresponding
                    388: *     components of Z are zero in this new basis.
                    389: *
                    390:       K = 0
                    391:       K2 = N + 1
                    392:       DO 70 J = 1, N
                    393:          IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
                    394: *
                    395: *           Deflate due to small z component.
                    396: *
                    397:             K2 = K2 - 1
                    398:             INDXP( K2 ) = J
                    399:             IF( J.EQ.N )
                    400:      $         GO TO 110
                    401:          ELSE
                    402:             JLAM = J
                    403:             GO TO 80
                    404:          END IF
                    405:    70 CONTINUE
                    406:    80 CONTINUE
                    407:       J = J + 1
                    408:       IF( J.GT.N )
                    409:      $   GO TO 100
                    410:       IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
                    411: *
                    412: *        Deflate due to small z component.
                    413: *
                    414:          K2 = K2 - 1
                    415:          INDXP( K2 ) = J
                    416:       ELSE
                    417: *
                    418: *        Check if eigenvalues are close enough to allow deflation.
                    419: *
                    420:          S = Z( JLAM )
                    421:          C = Z( J )
                    422: *
                    423: *        Find sqrt(a**2+b**2) without overflow or
                    424: *        destructive underflow.
                    425: *
                    426:          TAU = DLAPY2( C, S )
                    427:          T = D( J ) - D( JLAM )
                    428:          C = C / TAU
                    429:          S = -S / TAU
                    430:          IF( ABS( T*C*S ).LE.TOL ) THEN
                    431: *
                    432: *           Deflation is possible.
                    433: *
                    434:             Z( J ) = TAU
                    435:             Z( JLAM ) = ZERO
                    436: *
                    437: *           Record the appropriate Givens rotation
                    438: *
                    439:             GIVPTR = GIVPTR + 1
                    440:             GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
                    441:             GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
                    442:             GIVNUM( 1, GIVPTR ) = C
                    443:             GIVNUM( 2, GIVPTR ) = S
                    444:             IF( ICOMPQ.EQ.1 ) THEN
                    445:                CALL DROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
                    446:      $                    Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
                    447:             END IF
                    448:             T = D( JLAM )*C*C + D( J )*S*S
                    449:             D( J ) = D( JLAM )*S*S + D( J )*C*C
                    450:             D( JLAM ) = T
                    451:             K2 = K2 - 1
                    452:             I = 1
                    453:    90       CONTINUE
                    454:             IF( K2+I.LE.N ) THEN
                    455:                IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
                    456:                   INDXP( K2+I-1 ) = INDXP( K2+I )
                    457:                   INDXP( K2+I ) = JLAM
                    458:                   I = I + 1
                    459:                   GO TO 90
                    460:                ELSE
                    461:                   INDXP( K2+I-1 ) = JLAM
                    462:                END IF
                    463:             ELSE
                    464:                INDXP( K2+I-1 ) = JLAM
                    465:             END IF
                    466:             JLAM = J
                    467:          ELSE
                    468:             K = K + 1
                    469:             W( K ) = Z( JLAM )
                    470:             DLAMDA( K ) = D( JLAM )
                    471:             INDXP( K ) = JLAM
                    472:             JLAM = J
                    473:          END IF
                    474:       END IF
                    475:       GO TO 80
                    476:   100 CONTINUE
                    477: *
                    478: *     Record the last eigenvalue.
                    479: *
                    480:       K = K + 1
                    481:       W( K ) = Z( JLAM )
                    482:       DLAMDA( K ) = D( JLAM )
                    483:       INDXP( K ) = JLAM
                    484: *
                    485:   110 CONTINUE
                    486: *
                    487: *     Sort the eigenvalues and corresponding eigenvectors into DLAMDA
                    488: *     and Q2 respectively.  The eigenvalues/vectors which were not
                    489: *     deflated go into the first K slots of DLAMDA and Q2 respectively,
                    490: *     while those which were deflated go into the last N - K slots.
                    491: *
                    492:       IF( ICOMPQ.EQ.0 ) THEN
                    493:          DO 120 J = 1, N
                    494:             JP = INDXP( J )
                    495:             DLAMDA( J ) = D( JP )
                    496:             PERM( J ) = INDXQ( INDX( JP ) )
                    497:   120    CONTINUE
                    498:       ELSE
                    499:          DO 130 J = 1, N
                    500:             JP = INDXP( J )
                    501:             DLAMDA( J ) = D( JP )
                    502:             PERM( J ) = INDXQ( INDX( JP ) )
                    503:             CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
                    504:   130    CONTINUE
                    505:       END IF
                    506: *
                    507: *     The deflated eigenvalues and their corresponding vectors go back
                    508: *     into the last N - K slots of D and Q respectively.
                    509: *
                    510:       IF( K.LT.N ) THEN
                    511:          IF( ICOMPQ.EQ.0 ) THEN
                    512:             CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
                    513:          ELSE
                    514:             CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
                    515:             CALL DLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2,
                    516:      $                   Q( 1, K+1 ), LDQ )
                    517:          END IF
                    518:       END IF
                    519: *
                    520:       RETURN
                    521: *
                    522: *     End of DLAED8
                    523: *
                    524:       END

CVSweb interface <joel.bertrand@systella.fr>