SRC/zlaed8.f

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