SRC/chsein.f

   1 *> \brief \b CHSEIN
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CHSEIN + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chsein.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chsein.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chsein.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CHSEIN( 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 *       REAL               RWORK( * )
  33 *       COMPLEX            H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
  34 *      $                   W( * ), WORK( * )
  35 *       ..
  36 *
  37 *
  38 *> \par Purpose:
  39 *  =============
  40 *>
  41 *> \verbatim
  42 *>
  43 *> CHSEIN 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 CHSEQR; 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 CHSEIN 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, CHSEIN 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 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 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 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 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 array, dimension (N*N)
 184 *> \endverbatim
 185 *>
 186 *> \param[out] RWORK
 187 *> \verbatim
 188 *>          RWORK is REAL 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 complexOTHERcomputational
 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 CHSEIN( 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       REAL               RWORK( * )
 261       COMPLEX            H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
 262      $                   W( * ), WORK( * )
 263 *     ..
 264 *
 265 *  =====================================================================
 266 *
 267 *     .. Parameters ..
 268       COMPLEX            ZERO
 269       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) )
 270       REAL               RZERO
 271       PARAMETER          ( RZERO = 0.0E+0 )
 272 *     ..
 273 *     .. Local Scalars ..
 274       LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
 275       INTEGER            I, IINFO, K, KL, KLN, KR, KS, LDWORK
 276       REAL               EPS3, HNORM, SMLNUM, ULP, UNFL
 277       COMPLEX            CDUM, WK
 278 *     ..
 279 *     .. External Functions ..
 280       LOGICAL            LSAME, SISNAN
 281       REAL               CLANHS, SLAMCH
 282       EXTERNAL           LSAME, CLANHS, SLAMCH, SISNAN
 283 *     ..
 284 *     .. External Subroutines ..
 285       EXTERNAL           CLAEIN, XERBLA
 286 *     ..
 287 *     .. Intrinsic Functions ..
 288       INTRINSIC          ABS, AIMAG, MAX, REAL
 289 *     ..
 290 *     .. Statement Functions ..
 291       REAL               CABS1
 292 *     ..
 293 *     .. Statement Function definitions ..
 294       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( 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( 'CHSEIN', -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 = SLAMCH( 'Safe minimum' )
 348       ULP = SLAMCH( '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 = CLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
 403                IF( SISNAN( 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 CLAEIN( .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 CLAEIN( .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 CHSEIN
 467 *
 468       END