SRC/zsteqr.f

   1 *> \brief \b ZSTEQR
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZSTEQR + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsteqr.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsteqr.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsteqr.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          COMPZ
  25 *       INTEGER            INFO, LDZ, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
  29 *       COMPLEX*16         Z( LDZ, * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> ZSTEQR computes all eigenvalues and, optionally, eigenvectors of a
  39 *> symmetric tridiagonal matrix using the implicit QL or QR method.
  40 *> The eigenvectors of a full or band complex Hermitian matrix can also
  41 *> be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this
  42 *> matrix to tridiagonal form.
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in] COMPZ
  49 *> \verbatim
  50 *>          COMPZ is CHARACTER*1
  51 *>          = 'N':  Compute eigenvalues only.
  52 *>          = 'V':  Compute eigenvalues and eigenvectors of the original
  53 *>                  Hermitian matrix.  On entry, Z must contain the
  54 *>                  unitary matrix used to reduce the original matrix
  55 *>                  to tridiagonal form.
  56 *>          = 'I':  Compute eigenvalues and eigenvectors of the
  57 *>                  tridiagonal matrix.  Z is initialized to the identity
  58 *>                  matrix.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] N
  62 *> \verbatim
  63 *>          N is INTEGER
  64 *>          The order of the matrix.  N >= 0.
  65 *> \endverbatim
  66 *>
  67 *> \param[in,out] D
  68 *> \verbatim
  69 *>          D is DOUBLE PRECISION array, dimension (N)
  70 *>          On entry, the diagonal elements of the tridiagonal matrix.
  71 *>          On exit, if INFO = 0, the eigenvalues in ascending order.
  72 *> \endverbatim
  73 *>
  74 *> \param[in,out] E
  75 *> \verbatim
  76 *>          E is DOUBLE PRECISION array, dimension (N-1)
  77 *>          On entry, the (n-1) subdiagonal elements of the tridiagonal
  78 *>          matrix.
  79 *>          On exit, E has been destroyed.
  80 *> \endverbatim
  81 *>
  82 *> \param[in,out] Z
  83 *> \verbatim
  84 *>          Z is COMPLEX*16 array, dimension (LDZ, N)
  85 *>          On entry, if  COMPZ = 'V', then Z contains the unitary
  86 *>          matrix used in the reduction to tridiagonal form.
  87 *>          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
  88 *>          orthonormal eigenvectors of the original Hermitian matrix,
  89 *>          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
  90 *>          of the symmetric tridiagonal matrix.
  91 *>          If COMPZ = 'N', then Z is not referenced.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] LDZ
  95 *> \verbatim
  96 *>          LDZ is INTEGER
  97 *>          The leading dimension of the array Z.  LDZ >= 1, and if
  98 *>          eigenvectors are desired, then  LDZ >= max(1,N).
  99 *> \endverbatim
 100 *>
 101 *> \param[out] WORK
 102 *> \verbatim
 103 *>          WORK is DOUBLE PRECISION array, dimension (max(1,2*N-2))
 104 *>          If COMPZ = 'N', then WORK is not referenced.
 105 *> \endverbatim
 106 *>
 107 *> \param[out] INFO
 108 *> \verbatim
 109 *>          INFO is INTEGER
 110 *>          = 0:  successful exit
 111 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 112 *>          > 0:  the algorithm has failed to find all the eigenvalues in
 113 *>                a total of 30*N iterations; if INFO = i, then i
 114 *>                elements of E have not converged to zero; on exit, D
 115 *>                and E contain the elements of a symmetric tridiagonal
 116 *>                matrix which is unitarily similar to the original
 117 *>                matrix.
 118 *> \endverbatim
 119 *
 120 *  Authors:
 121 *  ========
 122 *
 123 *> \author Univ. of Tennessee
 124 *> \author Univ. of California Berkeley
 125 *> \author Univ. of Colorado Denver
 126 *> \author NAG Ltd.
 127 *
 128 *> \date November 2011
 129 *
 130 *> \ingroup complex16OTHERcomputational
 131 *
 132 *  =====================================================================
 133       SUBROUTINE ZSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
 134 *
 135 *  -- LAPACK computational routine (version 3.4.0) --
 136 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 137 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 138 *     November 2011
 139 *
 140 *     .. Scalar Arguments ..
 141       CHARACTER          COMPZ
 142       INTEGER            INFO, LDZ, N
 143 *     ..
 144 *     .. Array Arguments ..
 145       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
 146       COMPLEX*16         Z( LDZ, * )
 147 *     ..
 148 *
 149 *  =====================================================================
 150 *
 151 *     .. Parameters ..
 152       DOUBLE PRECISION   ZERO, ONE, TWO, THREE
 153       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
 154      $                   THREE = 3.0D0 )
 155       COMPLEX*16         CZERO, CONE
 156       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
 157      $                   CONE = ( 1.0D0, 0.0D0 ) )
 158       INTEGER            MAXIT
 159       PARAMETER          ( MAXIT = 30 )
 160 *     ..
 161 *     .. Local Scalars ..
 162       INTEGER            I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
 163      $                   LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
 164      $                   NM1, NMAXIT
 165       DOUBLE PRECISION   ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
 166      $                   S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
 167 *     ..
 168 *     .. External Functions ..
 169       LOGICAL            LSAME
 170       DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
 171       EXTERNAL           LSAME, DLAMCH, DLANST, DLAPY2
 172 *     ..
 173 *     .. External Subroutines ..
 174       EXTERNAL           DLAE2, DLAEV2, DLARTG, DLASCL, DLASRT, XERBLA,
 175      $                   ZLASET, ZLASR, ZSWAP
 176 *     ..
 177 *     .. Intrinsic Functions ..
 178       INTRINSIC          ABS, MAX, SIGN, SQRT
 179 *     ..
 180 *     .. Executable Statements ..
 181 *
 182 *     Test the input parameters.
 183 *
 184       INFO = 0
 185 *
 186       IF( LSAME( COMPZ, 'N' ) ) THEN
 187          ICOMPZ = 0
 188       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
 189          ICOMPZ = 1
 190       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
 191          ICOMPZ = 2
 192       ELSE
 193          ICOMPZ = -1
 194       END IF
 195       IF( ICOMPZ.LT.0 ) THEN
 196          INFO = -1
 197       ELSE IF( N.LT.0 ) THEN
 198          INFO = -2
 199       ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
 200      $         N ) ) ) THEN
 201          INFO = -6
 202       END IF
 203       IF( INFO.NE.0 ) THEN
 204          CALL XERBLA( 'ZSTEQR', -INFO )
 205          RETURN
 206       END IF
 207 *
 208 *     Quick return if possible
 209 *
 210       IF( N.EQ.0 )
 211      $   RETURN
 212 *
 213       IF( N.EQ.1 ) THEN
 214          IF( ICOMPZ.EQ.2 )
 215      $      Z( 1, 1 ) = CONE
 216          RETURN
 217       END IF
 218 *
 219 *     Determine the unit roundoff and over/underflow thresholds.
 220 *
 221       EPS = DLAMCH( 'E' )
 222       EPS2 = EPS**2
 223       SAFMIN = DLAMCH( 'S' )
 224       SAFMAX = ONE / SAFMIN
 225       SSFMAX = SQRT( SAFMAX ) / THREE
 226       SSFMIN = SQRT( SAFMIN ) / EPS2
 227 *
 228 *     Compute the eigenvalues and eigenvectors of the tridiagonal
 229 *     matrix.
 230 *
 231       IF( ICOMPZ.EQ.2 )
 232      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
 233 *
 234       NMAXIT = N*MAXIT
 235       JTOT = 0
 236 *
 237 *     Determine where the matrix splits and choose QL or QR iteration
 238 *     for each block, according to whether top or bottom diagonal
 239 *     element is smaller.
 240 *
 241       L1 = 1
 242       NM1 = N - 1
 243 *
 244    10 CONTINUE
 245       IF( L1.GT.N )
 246      $   GO TO 160
 247       IF( L1.GT.1 )
 248      $   E( L1-1 ) = ZERO
 249       IF( L1.LE.NM1 ) THEN
 250          DO 20 M = L1, NM1
 251             TST = ABS( E( M ) )
 252             IF( TST.EQ.ZERO )
 253      $         GO TO 30
 254             IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
 255      $          1 ) ) ) )*EPS ) THEN
 256                E( M ) = ZERO
 257                GO TO 30
 258             END IF
 259    20    CONTINUE
 260       END IF
 261       M = N
 262 *
 263    30 CONTINUE
 264       L = L1
 265       LSV = L
 266       LEND = M
 267       LENDSV = LEND
 268       L1 = M + 1
 269       IF( LEND.EQ.L )
 270      $   GO TO 10
 271 *
 272 *     Scale submatrix in rows and columns L to LEND
 273 *
 274       ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
 275       ISCALE = 0
 276       IF( ANORM.EQ.ZERO )
 277      $   GO TO 10
 278       IF( ANORM.GT.SSFMAX ) THEN
 279          ISCALE = 1
 280          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
 281      $                INFO )
 282          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
 283      $                INFO )
 284       ELSE IF( ANORM.LT.SSFMIN ) THEN
 285          ISCALE = 2
 286          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
 287      $                INFO )
 288          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
 289      $                INFO )
 290       END IF
 291 *
 292 *     Choose between QL and QR iteration
 293 *
 294       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
 295          LEND = LSV
 296          L = LENDSV
 297       END IF
 298 *
 299       IF( LEND.GT.L ) THEN
 300 *
 301 *        QL Iteration
 302 *
 303 *        Look for small subdiagonal element.
 304 *
 305    40    CONTINUE
 306          IF( L.NE.LEND ) THEN
 307             LENDM1 = LEND - 1
 308             DO 50 M = L, LENDM1
 309                TST = ABS( E( M ) )**2
 310                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
 311      $             SAFMIN )GO TO 60
 312    50       CONTINUE
 313          END IF
 314 *
 315          M = LEND
 316 *
 317    60    CONTINUE
 318          IF( M.LT.LEND )
 319      $      E( M ) = ZERO
 320          P = D( L )
 321          IF( M.EQ.L )
 322      $      GO TO 80
 323 *
 324 *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
 325 *        to compute its eigensystem.
 326 *
 327          IF( M.EQ.L+1 ) THEN
 328             IF( ICOMPZ.GT.0 ) THEN
 329                CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
 330                WORK( L ) = C
 331                WORK( N-1+L ) = S
 332                CALL ZLASR( 'R', 'V', 'B', N, 2, WORK( L ),
 333      $                     WORK( N-1+L ), Z( 1, L ), LDZ )
 334             ELSE
 335                CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
 336             END IF
 337             D( L ) = RT1
 338             D( L+1 ) = RT2
 339             E( L ) = ZERO
 340             L = L + 2
 341             IF( L.LE.LEND )
 342      $         GO TO 40
 343             GO TO 140
 344          END IF
 345 *
 346          IF( JTOT.EQ.NMAXIT )
 347      $      GO TO 140
 348          JTOT = JTOT + 1
 349 *
 350 *        Form shift.
 351 *
 352          G = ( D( L+1 )-P ) / ( TWO*E( L ) )
 353          R = DLAPY2( G, ONE )
 354          G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
 355 *
 356          S = ONE
 357          C = ONE
 358          P = ZERO
 359 *
 360 *        Inner loop
 361 *
 362          MM1 = M - 1
 363          DO 70 I = MM1, L, -1
 364             F = S*E( I )
 365             B = C*E( I )
 366             CALL DLARTG( G, F, C, S, R )
 367             IF( I.NE.M-1 )
 368      $         E( I+1 ) = R
 369             G = D( I+1 ) - P
 370             R = ( D( I )-G )*S + TWO*C*B
 371             P = S*R
 372             D( I+1 ) = G + P
 373             G = C*R - B
 374 *
 375 *           If eigenvectors are desired, then save rotations.
 376 *
 377             IF( ICOMPZ.GT.0 ) THEN
 378                WORK( I ) = C
 379                WORK( N-1+I ) = -S
 380             END IF
 381 *
 382    70    CONTINUE
 383 *
 384 *        If eigenvectors are desired, then apply saved rotations.
 385 *
 386          IF( ICOMPZ.GT.0 ) THEN
 387             MM = M - L + 1
 388             CALL ZLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ),
 389      $                  Z( 1, L ), LDZ )
 390          END IF
 391 *
 392          D( L ) = D( L ) - P
 393          E( L ) = G
 394          GO TO 40
 395 *
 396 *        Eigenvalue found.
 397 *
 398    80    CONTINUE
 399          D( L ) = P
 400 *
 401          L = L + 1
 402          IF( L.LE.LEND )
 403      $      GO TO 40
 404          GO TO 140
 405 *
 406       ELSE
 407 *
 408 *        QR Iteration
 409 *
 410 *        Look for small superdiagonal element.
 411 *
 412    90    CONTINUE
 413          IF( L.NE.LEND ) THEN
 414             LENDP1 = LEND + 1
 415             DO 100 M = L, LENDP1, -1
 416                TST = ABS( E( M-1 ) )**2
 417                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
 418      $             SAFMIN )GO TO 110
 419   100       CONTINUE
 420          END IF
 421 *
 422          M = LEND
 423 *
 424   110    CONTINUE
 425          IF( M.GT.LEND )
 426      $      E( M-1 ) = ZERO
 427          P = D( L )
 428          IF( M.EQ.L )
 429      $      GO TO 130
 430 *
 431 *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
 432 *        to compute its eigensystem.
 433 *
 434          IF( M.EQ.L-1 ) THEN
 435             IF( ICOMPZ.GT.0 ) THEN
 436                CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
 437                WORK( M ) = C
 438                WORK( N-1+M ) = S
 439                CALL ZLASR( 'R', 'V', 'F', N, 2, WORK( M ),
 440      $                     WORK( N-1+M ), Z( 1, L-1 ), LDZ )
 441             ELSE
 442                CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
 443             END IF
 444             D( L-1 ) = RT1
 445             D( L ) = RT2
 446             E( L-1 ) = ZERO
 447             L = L - 2
 448             IF( L.GE.LEND )
 449      $         GO TO 90
 450             GO TO 140
 451          END IF
 452 *
 453          IF( JTOT.EQ.NMAXIT )
 454      $      GO TO 140
 455          JTOT = JTOT + 1
 456 *
 457 *        Form shift.
 458 *
 459          G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
 460          R = DLAPY2( G, ONE )
 461          G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
 462 *
 463          S = ONE
 464          C = ONE
 465          P = ZERO
 466 *
 467 *        Inner loop
 468 *
 469          LM1 = L - 1
 470          DO 120 I = M, LM1
 471             F = S*E( I )
 472             B = C*E( I )
 473             CALL DLARTG( G, F, C, S, R )
 474             IF( I.NE.M )
 475      $         E( I-1 ) = R
 476             G = D( I ) - P
 477             R = ( D( I+1 )-G )*S + TWO*C*B
 478             P = S*R
 479             D( I ) = G + P
 480             G = C*R - B
 481 *
 482 *           If eigenvectors are desired, then save rotations.
 483 *
 484             IF( ICOMPZ.GT.0 ) THEN
 485                WORK( I ) = C
 486                WORK( N-1+I ) = S
 487             END IF
 488 *
 489   120    CONTINUE
 490 *
 491 *        If eigenvectors are desired, then apply saved rotations.
 492 *
 493          IF( ICOMPZ.GT.0 ) THEN
 494             MM = L - M + 1
 495             CALL ZLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
 496      $                  Z( 1, M ), LDZ )
 497          END IF
 498 *
 499          D( L ) = D( L ) - P
 500          E( LM1 ) = G
 501          GO TO 90
 502 *
 503 *        Eigenvalue found.
 504 *
 505   130    CONTINUE
 506          D( L ) = P
 507 *
 508          L = L - 1
 509          IF( L.GE.LEND )
 510      $      GO TO 90
 511          GO TO 140
 512 *
 513       END IF
 514 *
 515 *     Undo scaling if necessary
 516 *
 517   140 CONTINUE
 518       IF( ISCALE.EQ.1 ) THEN
 519          CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
 520      $                D( LSV ), N, INFO )
 521          CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
 522      $                N, INFO )
 523       ELSE IF( ISCALE.EQ.2 ) THEN
 524          CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
 525      $                D( LSV ), N, INFO )
 526          CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
 527      $                N, INFO )
 528       END IF
 529 *
 530 *     Check for no convergence to an eigenvalue after a total
 531 *     of N*MAXIT iterations.
 532 *
 533       IF( JTOT.EQ.NMAXIT ) THEN
 534          DO 150 I = 1, N - 1
 535             IF( E( I ).NE.ZERO )
 536      $         INFO = INFO + 1
 537   150    CONTINUE
 538          RETURN
 539       END IF
 540       GO TO 10
 541 *
 542 *     Order eigenvalues and eigenvectors.
 543 *
 544   160 CONTINUE
 545       IF( ICOMPZ.EQ.0 ) THEN
 546 *
 547 *        Use Quick Sort
 548 *
 549          CALL DLASRT( 'I', N, D, INFO )
 550 *
 551       ELSE
 552 *
 553 *        Use Selection Sort to minimize swaps of eigenvectors
 554 *
 555          DO 180 II = 2, N
 556             I = II - 1
 557             K = I
 558             P = D( I )
 559             DO 170 J = II, N
 560                IF( D( J ).LT.P ) THEN
 561                   K = J
 562                   P = D( J )
 563                END IF
 564   170       CONTINUE
 565             IF( K.NE.I ) THEN
 566                D( K ) = D( I )
 567                D( I ) = P
 568                CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
 569             END IF
 570   180    CONTINUE
 571       END IF
 572       RETURN
 573 *
 574 *     End of ZSTEQR
 575 *
 576       END