SRC/shsein.f

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