SRC/zhsein.f

   1 *> \brief \b ZHSEIN
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHSEIN + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhsein.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhsein.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhsein.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
  22 *                          LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
  23 *                          IFAILR, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          EIGSRC, INITV, SIDE
  27 *       INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       LOGICAL            SELECT( * )
  31 *       INTEGER            IFAILL( * ), IFAILR( * )
  32 *       DOUBLE PRECISION   RWORK( * )
  33 *       COMPLEX*16         H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
  34 *      $                   W( * ), WORK( * )
  35 *       ..
  36 *
  37 *
  38 *> \par Purpose:
  39 *  =============
  40 *>
  41 *> \verbatim
  42 *>
  43 *> ZHSEIN uses inverse iteration to find specified right and/or left
  44 *> eigenvectors of a complex upper Hessenberg matrix H.
  45 *>
  46 *> The right eigenvector x and the left eigenvector y of the matrix H
  47 *> corresponding to an eigenvalue w are defined by:
  48 *>
  49 *>              H * x = w * x,     y**h * H = w * y**h
  50 *>
  51 *> where y**h denotes the conjugate transpose of the vector y.
  52 *> \endverbatim
  53 *
  54 *  Arguments:
  55 *  ==========
  56 *
  57 *> \param[in] SIDE
  58 *> \verbatim
  59 *>          SIDE is CHARACTER*1
  60 *>          = 'R': compute right eigenvectors only;
  61 *>          = 'L': compute left eigenvectors only;
  62 *>          = 'B': compute both right and left eigenvectors.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] EIGSRC
  66 *> \verbatim
  67 *>          EIGSRC is CHARACTER*1
  68 *>          Specifies the source of eigenvalues supplied in W:
  69 *>          = 'Q': the eigenvalues were found using ZHSEQR; thus, if
  70 *>                 H has zero subdiagonal elements, and so is
  71 *>                 block-triangular, then the j-th eigenvalue can be
  72 *>                 assumed to be an eigenvalue of the block containing
  73 *>                 the j-th row/column.  This property allows ZHSEIN to
  74 *>                 perform inverse iteration on just one diagonal block.
  75 *>          = 'N': no assumptions are made on the correspondence
  76 *>                 between eigenvalues and diagonal blocks.  In this
  77 *>                 case, ZHSEIN must always perform inverse iteration
  78 *>                 using the whole matrix H.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] INITV
  82 *> \verbatim
  83 *>          INITV is CHARACTER*1
  84 *>          = 'N': no initial vectors are supplied;
  85 *>          = 'U': user-supplied initial vectors are stored in the arrays
  86 *>                 VL and/or VR.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] SELECT
  90 *> \verbatim
  91 *>          SELECT is LOGICAL array, dimension (N)
  92 *>          Specifies the eigenvectors to be computed. To select the
  93 *>          eigenvector corresponding to the eigenvalue W(j),
  94 *>          SELECT(j) must be set to .TRUE..
  95 *> \endverbatim
  96 *>
  97 *> \param[in] N
  98 *> \verbatim
  99 *>          N is INTEGER
 100 *>          The order of the matrix H.  N >= 0.
 101 *> \endverbatim
 102 *>
 103 *> \param[in] H
 104 *> \verbatim
 105 *>          H is COMPLEX*16 array, dimension (LDH,N)
 106 *>          The upper Hessenberg matrix H.
 107 *>          If a NaN is detected in H, the routine will return with INFO=-6.
 108 *> \endverbatim
 109 *>
 110 *> \param[in] LDH
 111 *> \verbatim
 112 *>          LDH is INTEGER
 113 *>          The leading dimension of the array H.  LDH >= max(1,N).
 114 *> \endverbatim
 115 *>
 116 *> \param[in,out] W
 117 *> \verbatim
 118 *>          W is COMPLEX*16 array, dimension (N)
 119 *>          On entry, the eigenvalues of H.
 120 *>          On exit, the real parts of W may have been altered since
 121 *>          close eigenvalues are perturbed slightly in searching for
 122 *>          independent eigenvectors.
 123 *> \endverbatim
 124 *>
 125 *> \param[in,out] VL
 126 *> \verbatim
 127 *>          VL is COMPLEX*16 array, dimension (LDVL,MM)
 128 *>          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
 129 *>          contain starting vectors for the inverse iteration for the
 130 *>          left eigenvectors; the starting vector for each eigenvector
 131 *>          must be in the same column in which the eigenvector will be
 132 *>          stored.
 133 *>          On exit, if SIDE = 'L' or 'B', the left eigenvectors
 134 *>          specified by SELECT will be stored consecutively in the
 135 *>          columns of VL, in the same order as their eigenvalues.
 136 *>          If SIDE = 'R', VL is not referenced.
 137 *> \endverbatim
 138 *>
 139 *> \param[in] LDVL
 140 *> \verbatim
 141 *>          LDVL is INTEGER
 142 *>          The leading dimension of the array VL.
 143 *>          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
 144 *> \endverbatim
 145 *>
 146 *> \param[in,out] VR
 147 *> \verbatim
 148 *>          VR is COMPLEX*16 array, dimension (LDVR,MM)
 149 *>          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
 150 *>          contain starting vectors for the inverse iteration for the
 151 *>          right eigenvectors; the starting vector for each eigenvector
 152 *>          must be in the same column in which the eigenvector will be
 153 *>          stored.
 154 *>          On exit, if SIDE = 'R' or 'B', the right eigenvectors
 155 *>          specified by SELECT will be stored consecutively in the
 156 *>          columns of VR, in the same order as their eigenvalues.
 157 *>          If SIDE = 'L', VR is not referenced.
 158 *> \endverbatim
 159 *>
 160 *> \param[in] LDVR
 161 *> \verbatim
 162 *>          LDVR is INTEGER
 163 *>          The leading dimension of the array VR.
 164 *>          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
 165 *> \endverbatim
 166 *>
 167 *> \param[in] MM
 168 *> \verbatim
 169 *>          MM is INTEGER
 170 *>          The number of columns in the arrays VL and/or VR. MM >= M.
 171 *> \endverbatim
 172 *>
 173 *> \param[out] M
 174 *> \verbatim
 175 *>          M is INTEGER
 176 *>          The number of columns in the arrays VL and/or VR required to
 177 *>          store the eigenvectors (= the number of .TRUE. elements in
 178 *>          SELECT).
 179 *> \endverbatim
 180 *>
 181 *> \param[out] WORK
 182 *> \verbatim
 183 *>          WORK is COMPLEX*16 array, dimension (N*N)
 184 *> \endverbatim
 185 *>
 186 *> \param[out] RWORK
 187 *> \verbatim
 188 *>          RWORK is DOUBLE PRECISION array, dimension (N)
 189 *> \endverbatim
 190 *>
 191 *> \param[out] IFAILL
 192 *> \verbatim
 193 *>          IFAILL is INTEGER array, dimension (MM)
 194 *>          If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
 195 *>          eigenvector in the i-th column of VL (corresponding to the
 196 *>          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
 197 *>          eigenvector converged satisfactorily.
 198 *>          If SIDE = 'R', IFAILL is not referenced.
 199 *> \endverbatim
 200 *>
 201 *> \param[out] IFAILR
 202 *> \verbatim
 203 *>          IFAILR is INTEGER array, dimension (MM)
 204 *>          If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
 205 *>          eigenvector in the i-th column of VR (corresponding to the
 206 *>          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
 207 *>          eigenvector converged satisfactorily.
 208 *>          If SIDE = 'L', IFAILR is not referenced.
 209 *> \endverbatim
 210 *>
 211 *> \param[out] INFO
 212 *> \verbatim
 213 *>          INFO is INTEGER
 214 *>          = 0:  successful exit
 215 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 216 *>          > 0:  if INFO = i, i is the number of eigenvectors which
 217 *>                failed to converge; see IFAILL and IFAILR for further
 218 *>                details.
 219 *> \endverbatim
 220 *
 221 *  Authors:
 222 *  ========
 223 *
 224 *> \author Univ. of Tennessee
 225 *> \author Univ. of California Berkeley
 226 *> \author Univ. of Colorado Denver
 227 *> \author NAG Ltd.
 228 *
 229 *> \date November 2013
 230 *
 231 *> \ingroup complex16OTHERcomputational
 232 *
 233 *> \par Further Details:
 234 *  =====================
 235 *>
 236 *> \verbatim
 237 *>
 238 *>  Each eigenvector is normalized so that the element of largest
 239 *>  magnitude has magnitude 1; here the magnitude of a complex number
 240 *>  (x,y) is taken to be |x|+|y|.
 241 *> \endverbatim
 242 *>
 243 *  =====================================================================
 244       SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
 245      $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
 246      $                   IFAILR, INFO )
 247 *
 248 *  -- LAPACK computational routine (version 3.5.0) --
 249 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 250 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 251 *     November 2013
 252 *
 253 *     .. Scalar Arguments ..
 254       CHARACTER          EIGSRC, INITV, SIDE
 255       INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
 256 *     ..
 257 *     .. Array Arguments ..
 258       LOGICAL            SELECT( * )
 259       INTEGER            IFAILL( * ), IFAILR( * )
 260       DOUBLE PRECISION   RWORK( * )
 261       COMPLEX*16         H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
 262      $                   W( * ), WORK( * )
 263 *     ..
 264 *
 265 *  =====================================================================
 266 *
 267 *     .. Parameters ..
 268       COMPLEX*16         ZERO
 269       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 270       DOUBLE PRECISION   RZERO
 271       PARAMETER          ( RZERO = 0.0D+0 )
 272 *     ..
 273 *     .. Local Scalars ..
 274       LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
 275       INTEGER            I, IINFO, K, KL, KLN, KR, KS, LDWORK
 276       DOUBLE PRECISION   EPS3, HNORM, SMLNUM, ULP, UNFL
 277       COMPLEX*16         CDUM, WK
 278 *     ..
 279 *     .. External Functions ..
 280       LOGICAL            LSAME, DISNAN
 281       DOUBLE PRECISION   DLAMCH, ZLANHS
 282       EXTERNAL           LSAME, DLAMCH, ZLANHS, DISNAN
 283 *     ..
 284 *     .. External Subroutines ..
 285       EXTERNAL           XERBLA, ZLAEIN
 286 *     ..
 287 *     .. Intrinsic Functions ..
 288       INTRINSIC          ABS, DBLE, DIMAG, MAX
 289 *     ..
 290 *     .. Statement Functions ..
 291       DOUBLE PRECISION   CABS1
 292 *     ..
 293 *     .. Statement Function definitions ..
 294       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
 295 *     ..
 296 *     .. Executable Statements ..
 297 *
 298 *     Decode and test the input parameters.
 299 *
 300       BOTHV = LSAME( SIDE, 'B' )
 301       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
 302       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
 303 *
 304       FROMQR = LSAME( EIGSRC, 'Q' )
 305 *
 306       NOINIT = LSAME( INITV, 'N' )
 307 *
 308 *     Set M to the number of columns required to store the selected
 309 *     eigenvectors.
 310 *
 311       M = 0
 312       DO 10 K = 1, N
 313          IF( SELECT( K ) )
 314      $      M = M + 1
 315    10 CONTINUE
 316 *
 317       INFO = 0
 318       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
 319          INFO = -1
 320       ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
 321          INFO = -2
 322       ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
 323          INFO = -3
 324       ELSE IF( N.LT.0 ) THEN
 325          INFO = -5
 326       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
 327          INFO = -7
 328       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
 329          INFO = -10
 330       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
 331          INFO = -12
 332       ELSE IF( MM.LT.M ) THEN
 333          INFO = -13
 334       END IF
 335       IF( INFO.NE.0 ) THEN
 336          CALL XERBLA( 'ZHSEIN', -INFO )
 337          RETURN
 338       END IF
 339 *
 340 *     Quick return if possible.
 341 *
 342       IF( N.EQ.0 )
 343      $   RETURN
 344 *
 345 *     Set machine-dependent constants.
 346 *
 347       UNFL = DLAMCH( 'Safe minimum' )
 348       ULP = DLAMCH( 'Precision' )
 349       SMLNUM = UNFL*( N / ULP )
 350 *
 351       LDWORK = N
 352 *
 353       KL = 1
 354       KLN = 0
 355       IF( FROMQR ) THEN
 356          KR = 0
 357       ELSE
 358          KR = N
 359       END IF
 360       KS = 1
 361 *
 362       DO 100 K = 1, N
 363          IF( SELECT( K ) ) THEN
 364 *
 365 *           Compute eigenvector(s) corresponding to W(K).
 366 *
 367             IF( FROMQR ) THEN
 368 *
 369 *              If affiliation of eigenvalues is known, check whether
 370 *              the matrix splits.
 371 *
 372 *              Determine KL and KR such that 1 <= KL <= K <= KR <= N
 373 *              and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
 374 *              KR = N).
 375 *
 376 *              Then inverse iteration can be performed with the
 377 *              submatrix H(KL:N,KL:N) for a left eigenvector, and with
 378 *              the submatrix H(1:KR,1:KR) for a right eigenvector.
 379 *
 380                DO 20 I = K, KL + 1, -1
 381                   IF( H( I, I-1 ).EQ.ZERO )
 382      $               GO TO 30
 383    20          CONTINUE
 384    30          CONTINUE
 385                KL = I
 386                IF( K.GT.KR ) THEN
 387                   DO 40 I = K, N - 1
 388                      IF( H( I+1, I ).EQ.ZERO )
 389      $                  GO TO 50
 390    40             CONTINUE
 391    50             CONTINUE
 392                   KR = I
 393                END IF
 394             END IF
 395 *
 396             IF( KL.NE.KLN ) THEN
 397                KLN = KL
 398 *
 399 *              Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
 400 *              has not ben computed before.
 401 *
 402                HNORM = ZLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
 403                IF( DISNAN( HNORM ) ) THEN
 404                   INFO = -6
 405                   RETURN
 406                ELSE IF( HNORM.GT.RZERO ) THEN
 407                   EPS3 = HNORM*ULP
 408                ELSE
 409                   EPS3 = SMLNUM
 410                END IF
 411             END IF
 412 *
 413 *           Perturb eigenvalue if it is close to any previous
 414 *           selected eigenvalues affiliated to the submatrix
 415 *           H(KL:KR,KL:KR). Close roots are modified by EPS3.
 416 *
 417             WK = W( K )
 418    60       CONTINUE
 419             DO 70 I = K - 1, KL, -1
 420                IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN
 421                   WK = WK + EPS3
 422                   GO TO 60
 423                END IF
 424    70       CONTINUE
 425             W( K ) = WK
 426 *
 427             IF( LEFTV ) THEN
 428 *
 429 *              Compute left eigenvector.
 430 *
 431                CALL ZLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
 432      $                      WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3,
 433      $                      SMLNUM, IINFO )
 434                IF( IINFO.GT.0 ) THEN
 435                   INFO = INFO + 1
 436                   IFAILL( KS ) = K
 437                ELSE
 438                   IFAILL( KS ) = 0
 439                END IF
 440                DO 80 I = 1, KL - 1
 441                   VL( I, KS ) = ZERO
 442    80          CONTINUE
 443             END IF
 444             IF( RIGHTV ) THEN
 445 *
 446 *              Compute right eigenvector.
 447 *
 448                CALL ZLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ),
 449      $                      WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO )
 450                IF( IINFO.GT.0 ) THEN
 451                   INFO = INFO + 1
 452                   IFAILR( KS ) = K
 453                ELSE
 454                   IFAILR( KS ) = 0
 455                END IF
 456                DO 90 I = KR + 1, N
 457                   VR( I, KS ) = ZERO
 458    90          CONTINUE
 459             END IF
 460             KS = KS + 1
 461          END IF
 462   100 CONTINUE
 463 *
 464       RETURN
 465 *
 466 *     End of ZHSEIN
 467 *
 468       END