Return to zlaed8.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack vers la version 3.3.0.

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