LAPACK 3.3.1 Linear Algebra PACKage

# shsein.f

Go to the documentation of this file.
```00001       SUBROUTINE SHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI,
00002      \$                   VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL,
00003      \$                   IFAILR, INFO )
00004 *
00005 *  -- LAPACK routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          EIGSRC, INITV, SIDE
00012       INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            SELECT( * )
00016       INTEGER            IFAILL( * ), IFAILR( * )
00017       REAL               H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
00018      \$                   WI( * ), WORK( * ), WR( * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  SHSEIN uses inverse iteration to find specified right and/or left
00025 *  eigenvectors of a real upper Hessenberg matrix H.
00026 *
00027 *  The right eigenvector x and the left eigenvector y of the matrix H
00028 *  corresponding to an eigenvalue w are defined by:
00029 *
00030 *               H * x = w * x,     y**h * H = w * y**h
00031 *
00032 *  where y**h denotes the conjugate transpose of the vector y.
00033 *
00034 *  Arguments
00035 *  =========
00036 *
00037 *  SIDE    (input) CHARACTER*1
00038 *          = 'R': compute right eigenvectors only;
00039 *          = 'L': compute left eigenvectors only;
00040 *          = 'B': compute both right and left eigenvectors.
00041 *
00042 *  EIGSRC  (input) CHARACTER*1
00043 *          Specifies the source of eigenvalues supplied in (WR,WI):
00044 *          = 'Q': the eigenvalues were found using SHSEQR; thus, if
00045 *                 H has zero subdiagonal elements, and so is
00046 *                 block-triangular, then the j-th eigenvalue can be
00047 *                 assumed to be an eigenvalue of the block containing
00048 *                 the j-th row/column.  This property allows SHSEIN to
00049 *                 perform inverse iteration on just one diagonal block.
00050 *          = 'N': no assumptions are made on the correspondence
00051 *                 between eigenvalues and diagonal blocks.  In this
00052 *                 case, SHSEIN must always perform inverse iteration
00053 *                 using the whole matrix H.
00054 *
00055 *  INITV   (input) CHARACTER*1
00056 *          = 'N': no initial vectors are supplied;
00057 *          = 'U': user-supplied initial vectors are stored in the arrays
00058 *                 VL and/or VR.
00059 *
00060 *  SELECT  (input/output) LOGICAL array, dimension (N)
00061 *          Specifies the eigenvectors to be computed. To select the
00062 *          real eigenvector corresponding to a real eigenvalue WR(j),
00063 *          SELECT(j) must be set to .TRUE.. To select the complex
00064 *          eigenvector corresponding to a complex eigenvalue
00065 *          (WR(j),WI(j)), with complex conjugate (WR(j+1),WI(j+1)),
00066 *          either SELECT(j) or SELECT(j+1) or both must be set to
00067 *          .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is
00068 *          .FALSE..
00069 *
00070 *  N       (input) INTEGER
00071 *          The order of the matrix H.  N >= 0.
00072 *
00073 *  H       (input) REAL array, dimension (LDH,N)
00074 *          The upper Hessenberg matrix H.
00075 *
00076 *  LDH     (input) INTEGER
00077 *          The leading dimension of the array H.  LDH >= max(1,N).
00078 *
00079 *  WR      (input/output) REAL array, dimension (N)
00080 *  WI      (input) REAL array, dimension (N)
00081 *          On entry, the real and imaginary parts of the eigenvalues of
00082 *          H; a complex conjugate pair of eigenvalues must be stored in
00083 *          consecutive elements of WR and WI.
00084 *          On exit, WR may have been altered since close eigenvalues
00085 *          are perturbed slightly in searching for independent
00086 *          eigenvectors.
00087 *
00088 *  VL      (input/output) REAL array, dimension (LDVL,MM)
00089 *          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
00090 *          contain starting vectors for the inverse iteration for the
00091 *          left eigenvectors; the starting vector for each eigenvector
00092 *          must be in the same column(s) in which the eigenvector will
00093 *          be stored.
00094 *          On exit, if SIDE = 'L' or 'B', the left eigenvectors
00095 *          specified by SELECT will be stored consecutively in the
00096 *          columns of VL, in the same order as their eigenvalues. A
00097 *          complex eigenvector corresponding to a complex eigenvalue is
00098 *          stored in two consecutive columns, the first holding the real
00099 *          part and the second the imaginary part.
00100 *          If SIDE = 'R', VL is not referenced.
00101 *
00102 *  LDVL    (input) INTEGER
00103 *          The leading dimension of the array VL.
00104 *          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
00105 *
00106 *  VR      (input/output) REAL array, dimension (LDVR,MM)
00107 *          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
00108 *          contain starting vectors for the inverse iteration for the
00109 *          right eigenvectors; the starting vector for each eigenvector
00110 *          must be in the same column(s) in which the eigenvector will
00111 *          be stored.
00112 *          On exit, if SIDE = 'R' or 'B', the right eigenvectors
00113 *          specified by SELECT will be stored consecutively in the
00114 *          columns of VR, in the same order as their eigenvalues. A
00115 *          complex eigenvector corresponding to a complex eigenvalue is
00116 *          stored in two consecutive columns, the first holding the real
00117 *          part and the second the imaginary part.
00118 *          If SIDE = 'L', VR is not referenced.
00119 *
00120 *  LDVR    (input) INTEGER
00121 *          The leading dimension of the array VR.
00122 *          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
00123 *
00124 *  MM      (input) INTEGER
00125 *          The number of columns in the arrays VL and/or VR. MM >= M.
00126 *
00127 *  M       (output) INTEGER
00128 *          The number of columns in the arrays VL and/or VR required to
00129 *          store the eigenvectors; each selected real eigenvector
00130 *          occupies one column and each selected complex eigenvector
00131 *          occupies two columns.
00132 *
00133 *  WORK    (workspace) REAL array, dimension ((N+2)*N)
00134 *
00135 *  IFAILL  (output) INTEGER array, dimension (MM)
00136 *          If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
00137 *          eigenvector in the i-th column of VL (corresponding to the
00138 *          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
00139 *          eigenvector converged satisfactorily. If the i-th and (i+1)th
00140 *          columns of VL hold a complex eigenvector, then IFAILL(i) and
00141 *          IFAILL(i+1) are set to the same value.
00142 *          If SIDE = 'R', IFAILL is not referenced.
00143 *
00144 *  IFAILR  (output) INTEGER array, dimension (MM)
00145 *          If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
00146 *          eigenvector in the i-th column of VR (corresponding to the
00147 *          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
00148 *          eigenvector converged satisfactorily. If the i-th and (i+1)th
00149 *          columns of VR hold a complex eigenvector, then IFAILR(i) and
00150 *          IFAILR(i+1) are set to the same value.
00151 *          If SIDE = 'L', IFAILR is not referenced.
00152 *
00153 *  INFO    (output) INTEGER
00154 *          = 0:  successful exit
00155 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00156 *          > 0:  if INFO = i, i is the number of eigenvectors which
00157 *                failed to converge; see IFAILL and IFAILR for further
00158 *                details.
00159 *
00160 *  Further Details
00161 *  ===============
00162 *
00163 *  Each eigenvector is normalized so that the element of largest
00164 *  magnitude has magnitude 1; here the magnitude of a complex number
00165 *  (x,y) is taken to be |x|+|y|.
00166 *
00167 *  =====================================================================
00168 *
00169 *     .. Parameters ..
00170       REAL               ZERO, ONE
00171       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00172 *     ..
00173 *     .. Local Scalars ..
00174       LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, PAIR, RIGHTV
00175       INTEGER            I, IINFO, K, KL, KLN, KR, KSI, KSR, LDWORK
00176       REAL               BIGNUM, EPS3, HNORM, SMLNUM, ULP, UNFL, WKI,
00177      \$                   WKR
00178 *     ..
00179 *     .. External Functions ..
00180       LOGICAL            LSAME
00181       REAL               SLAMCH, SLANHS
00182       EXTERNAL           LSAME, SLAMCH, SLANHS
00183 *     ..
00184 *     .. External Subroutines ..
00185       EXTERNAL           SLAEIN, XERBLA
00186 *     ..
00187 *     .. Intrinsic Functions ..
00188       INTRINSIC          ABS, MAX
00189 *     ..
00190 *     .. Executable Statements ..
00191 *
00192 *     Decode and test the input parameters.
00193 *
00194       BOTHV = LSAME( SIDE, 'B' )
00195       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
00196       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
00197 *
00198       FROMQR = LSAME( EIGSRC, 'Q' )
00199 *
00200       NOINIT = LSAME( INITV, 'N' )
00201 *
00202 *     Set M to the number of columns required to store the selected
00203 *     eigenvectors, and standardize the array SELECT.
00204 *
00205       M = 0
00206       PAIR = .FALSE.
00207       DO 10 K = 1, N
00208          IF( PAIR ) THEN
00209             PAIR = .FALSE.
00210             SELECT( K ) = .FALSE.
00211          ELSE
00212             IF( WI( K ).EQ.ZERO ) THEN
00213                IF( SELECT( K ) )
00214      \$            M = M + 1
00215             ELSE
00216                PAIR = .TRUE.
00217                IF( SELECT( K ) .OR. SELECT( K+1 ) ) THEN
00218                   SELECT( K ) = .TRUE.
00219                   M = M + 2
00220                END IF
00221             END IF
00222          END IF
00223    10 CONTINUE
00224 *
00225       INFO = 0
00226       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
00227          INFO = -1
00228       ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
00229          INFO = -2
00230       ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
00231          INFO = -3
00232       ELSE IF( N.LT.0 ) THEN
00233          INFO = -5
00234       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
00235          INFO = -7
00236       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
00237          INFO = -11
00238       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
00239          INFO = -13
00240       ELSE IF( MM.LT.M ) THEN
00241          INFO = -14
00242       END IF
00243       IF( INFO.NE.0 ) THEN
00244          CALL XERBLA( 'SHSEIN', -INFO )
00245          RETURN
00246       END IF
00247 *
00248 *     Quick return if possible.
00249 *
00250       IF( N.EQ.0 )
00251      \$   RETURN
00252 *
00253 *     Set machine-dependent constants.
00254 *
00255       UNFL = SLAMCH( 'Safe minimum' )
00256       ULP = SLAMCH( 'Precision' )
00257       SMLNUM = UNFL*( N / ULP )
00258       BIGNUM = ( ONE-ULP ) / SMLNUM
00259 *
00260       LDWORK = N + 1
00261 *
00262       KL = 1
00263       KLN = 0
00264       IF( FROMQR ) THEN
00265          KR = 0
00266       ELSE
00267          KR = N
00268       END IF
00269       KSR = 1
00270 *
00271       DO 120 K = 1, N
00272          IF( SELECT( K ) ) THEN
00273 *
00274 *           Compute eigenvector(s) corresponding to W(K).
00275 *
00276             IF( FROMQR ) THEN
00277 *
00278 *              If affiliation of eigenvalues is known, check whether
00279 *              the matrix splits.
00280 *
00281 *              Determine KL and KR such that 1 <= KL <= K <= KR <= N
00282 *              and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
00283 *              KR = N).
00284 *
00285 *              Then inverse iteration can be performed with the
00286 *              submatrix H(KL:N,KL:N) for a left eigenvector, and with
00287 *              the submatrix H(1:KR,1:KR) for a right eigenvector.
00288 *
00289                DO 20 I = K, KL + 1, -1
00290                   IF( H( I, I-1 ).EQ.ZERO )
00291      \$               GO TO 30
00292    20          CONTINUE
00293    30          CONTINUE
00294                KL = I
00295                IF( K.GT.KR ) THEN
00296                   DO 40 I = K, N - 1
00297                      IF( H( I+1, I ).EQ.ZERO )
00298      \$                  GO TO 50
00299    40             CONTINUE
00300    50             CONTINUE
00301                   KR = I
00302                END IF
00303             END IF
00304 *
00305             IF( KL.NE.KLN ) THEN
00306                KLN = KL
00307 *
00308 *              Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
00309 *              has not ben computed before.
00310 *
00311                HNORM = SLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, WORK )
00312                IF( HNORM.GT.ZERO ) THEN
00313                   EPS3 = HNORM*ULP
00314                ELSE
00315                   EPS3 = SMLNUM
00316                END IF
00317             END IF
00318 *
00319 *           Perturb eigenvalue if it is close to any previous
00320 *           selected eigenvalues affiliated to the submatrix
00321 *           H(KL:KR,KL:KR). Close roots are modified by EPS3.
00322 *
00323             WKR = WR( K )
00324             WKI = WI( K )
00325    60       CONTINUE
00326             DO 70 I = K - 1, KL, -1
00327                IF( SELECT( I ) .AND. ABS( WR( I )-WKR )+
00328      \$             ABS( WI( I )-WKI ).LT.EPS3 ) THEN
00329                   WKR = WKR + EPS3
00330                   GO TO 60
00331                END IF
00332    70       CONTINUE
00333             WR( K ) = WKR
00334 *
00335             PAIR = WKI.NE.ZERO
00336             IF( PAIR ) THEN
00337                KSI = KSR + 1
00338             ELSE
00339                KSI = KSR
00340             END IF
00341             IF( LEFTV ) THEN
00342 *
00343 *              Compute left eigenvector.
00344 *
00345                CALL SLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
00346      \$                      WKR, WKI, VL( KL, KSR ), VL( KL, KSI ),
00347      \$                      WORK, LDWORK, WORK( N*N+N+1 ), EPS3, SMLNUM,
00348      \$                      BIGNUM, IINFO )
00349                IF( IINFO.GT.0 ) THEN
00350                   IF( PAIR ) THEN
00351                      INFO = INFO + 2
00352                   ELSE
00353                      INFO = INFO + 1
00354                   END IF
00355                   IFAILL( KSR ) = K
00356                   IFAILL( KSI ) = K
00357                ELSE
00358                   IFAILL( KSR ) = 0
00359                   IFAILL( KSI ) = 0
00360                END IF
00361                DO 80 I = 1, KL - 1
00362                   VL( I, KSR ) = ZERO
00363    80          CONTINUE
00364                IF( PAIR ) THEN
00365                   DO 90 I = 1, KL - 1
00366                      VL( I, KSI ) = ZERO
00367    90             CONTINUE
00368                END IF
00369             END IF
00370             IF( RIGHTV ) THEN
00371 *
00372 *              Compute right eigenvector.
00373 *
00374                CALL SLAEIN( .TRUE., NOINIT, KR, H, LDH, WKR, WKI,
00375      \$                      VR( 1, KSR ), VR( 1, KSI ), WORK, LDWORK,
00376      \$                      WORK( N*N+N+1 ), EPS3, SMLNUM, BIGNUM,
00377      \$                      IINFO )
00378                IF( IINFO.GT.0 ) THEN
00379                   IF( PAIR ) THEN
00380                      INFO = INFO + 2
00381                   ELSE
00382                      INFO = INFO + 1
00383                   END IF
00384                   IFAILR( KSR ) = K
00385                   IFAILR( KSI ) = K
00386                ELSE
00387                   IFAILR( KSR ) = 0
00388                   IFAILR( KSI ) = 0
00389                END IF
00390                DO 100 I = KR + 1, N
00391                   VR( I, KSR ) = ZERO
00392   100          CONTINUE
00393                IF( PAIR ) THEN
00394                   DO 110 I = KR + 1, N
00395                      VR( I, KSI ) = ZERO
00396   110             CONTINUE
00397                END IF
00398             END IF
00399 *
00400             IF( PAIR ) THEN
00401                KSR = KSR + 2
00402             ELSE
00403                KSR = KSR + 1
00404             END IF
00405          END IF
00406   120 CONTINUE
00407 *
00408       RETURN
00409 *
00410 *     End of SHSEIN
00411 *
00412       END
```