SRC/dlaed8.f

   1 *> \brief \b DLAED8 used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense.
   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 *
 231 *> \date September 2012
 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 *  =====================================================================
 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 *
 246 *  -- LAPACK computational routine (version 3.4.2) --
 247 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 248 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 249 *     September 2012
 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 *
 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 *
 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