```      SUBROUTINE SLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
\$                   LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
*
*  -- LAPACK auxiliary routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
LOGICAL            NOINIT, RIGHTV
INTEGER            INFO, LDB, LDH, N
REAL               BIGNUM, EPS3, SMLNUM, WI, WR
*     ..
*     .. Array Arguments ..
REAL               B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
\$                   WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  SLAEIN uses inverse iteration to find a right or left eigenvector
*  corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
*  matrix H.
*
*  Arguments
*  =========
*
*  RIGHTV   (input) LOGICAL
*          = .TRUE. : compute right eigenvector;
*          = .FALSE.: compute left eigenvector.
*
*  NOINIT   (input) LOGICAL
*          = .TRUE. : no initial vector supplied in (VR,VI).
*          = .FALSE.: initial vector supplied in (VR,VI).
*
*  N       (input) INTEGER
*          The order of the matrix H.  N >= 0.
*
*  H       (input) REAL array, dimension (LDH,N)
*          The upper Hessenberg matrix H.
*
*  LDH     (input) INTEGER
*          The leading dimension of the array H.  LDH >= max(1,N).
*
*  WR      (input) REAL
*  WI      (input) REAL
*          The real and imaginary parts of the eigenvalue of H whose
*          corresponding right or left eigenvector is to be computed.
*
*  VR      (input/output) REAL array, dimension (N)
*  VI      (input/output) REAL array, dimension (N)
*          On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
*          a real starting vector for inverse iteration using the real
*          eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
*          must contain the real and imaginary parts of a complex
*          starting vector for inverse iteration using the complex
*          eigenvalue (WR,WI); otherwise VR and VI need not be set.
*          On exit, if WI = 0.0 (real eigenvalue), VR contains the
*          computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
*          VR and VI contain the real and imaginary parts of the
*          computed complex eigenvector. The eigenvector is normalized
*          so that the component of largest magnitude has magnitude 1;
*          here the magnitude of a complex number (x,y) is taken to be
*          |x| + |y|.
*          VI is not referenced if WI = 0.0.
*
*  B       (workspace) REAL array, dimension (LDB,N)
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B.  LDB >= N+1.
*
*  WORK   (workspace) REAL array, dimension (N)
*
*  EPS3    (input) REAL
*          A small machine-dependent value which is used to perturb
*          close eigenvalues, and to replace zero pivots.
*
*  SMLNUM  (input) REAL
*          A machine-dependent value close to the underflow threshold.
*
*  BIGNUM  (input) REAL
*          A machine-dependent value close to the overflow threshold.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          = 1:  inverse iteration did not converge; VR is set to the
*                last iterate, and so is VI if WI.ne.0.0.
*
*  =====================================================================
*
*     .. Parameters ..
REAL               ZERO, ONE, TENTH
PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TENTH = 1.0E-1 )
*     ..
*     .. Local Scalars ..
CHARACTER          NORMIN, TRANS
INTEGER            I, I1, I2, I3, IERR, ITS, J
REAL               ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
\$                   REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W,
\$                   W1, X, XI, XR, Y
*     ..
*     .. External Functions ..
INTEGER            ISAMAX
REAL               SASUM, SLAPY2, SNRM2
EXTERNAL           ISAMAX, SASUM, SLAPY2, SNRM2
*     ..
*     .. External Subroutines ..
*     ..
*     .. Intrinsic Functions ..
INTRINSIC          ABS, MAX, REAL, SQRT
*     ..
*     .. Executable Statements ..
*
INFO = 0
*
*     GROWTO is the threshold used in the acceptance test for an
*     eigenvector.
*
ROOTN = SQRT( REAL( N ) )
GROWTO = TENTH / ROOTN
NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
*
*     Form B = H - (WR,WI)*I (except that the subdiagonal elements and
*     the imaginary parts of the diagonal elements are not stored).
*
DO 20 J = 1, N
DO 10 I = 1, J - 1
B( I, J ) = H( I, J )
10    CONTINUE
B( J, J ) = H( J, J ) - WR
20 CONTINUE
*
IF( WI.EQ.ZERO ) THEN
*
*        Real eigenvalue.
*
IF( NOINIT ) THEN
*
*           Set initial vector.
*
DO 30 I = 1, N
VR( I ) = EPS3
30       CONTINUE
ELSE
*
*           Scale supplied initial vector.
*
VNORM = SNRM2( N, VR, 1 )
CALL SSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR,
\$                  1 )
END IF
*
IF( RIGHTV ) THEN
*
*           LU decomposition with partial pivoting of B, replacing zero
*           pivots by EPS3.
*
DO 60 I = 1, N - 1
EI = H( I+1, I )
IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN
*
*                 Interchange rows and eliminate.
*
X = B( I, I ) / EI
B( I, I ) = EI
DO 40 J = I + 1, N
TEMP = B( I+1, J )
B( I+1, J ) = B( I, J ) - X*TEMP
B( I, J ) = TEMP
40             CONTINUE
ELSE
*
*                 Eliminate without interchange.
*
IF( B( I, I ).EQ.ZERO )
\$               B( I, I ) = EPS3
X = EI / B( I, I )
IF( X.NE.ZERO ) THEN
DO 50 J = I + 1, N
B( I+1, J ) = B( I+1, J ) - X*B( I, J )
50                CONTINUE
END IF
END IF
60       CONTINUE
IF( B( N, N ).EQ.ZERO )
\$         B( N, N ) = EPS3
*
TRANS = 'N'
*
ELSE
*
*           UL decomposition with partial pivoting of B, replacing zero
*           pivots by EPS3.
*
DO 90 J = N, 2, -1
EJ = H( J, J-1 )
IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
*
*                 Interchange columns and eliminate.
*
X = B( J, J ) / EJ
B( J, J ) = EJ
DO 70 I = 1, J - 1
TEMP = B( I, J-1 )
B( I, J-1 ) = B( I, J ) - X*TEMP
B( I, J ) = TEMP
70             CONTINUE
ELSE
*
*                 Eliminate without interchange.
*
IF( B( J, J ).EQ.ZERO )
\$               B( J, J ) = EPS3
X = EJ / B( J, J )
IF( X.NE.ZERO ) THEN
DO 80 I = 1, J - 1
B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
80                CONTINUE
END IF
END IF
90       CONTINUE
IF( B( 1, 1 ).EQ.ZERO )
\$         B( 1, 1 ) = EPS3
*
TRANS = 'T'
*
END IF
*
NORMIN = 'N'
DO 110 ITS = 1, N
*
*           Solve U*x = scale*v for a right eigenvector
*             or U'*x = scale*v for a left eigenvector,
*           overwriting x on v.
*
CALL SLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB,
\$                   VR, SCALE, WORK, IERR )
NORMIN = 'Y'
*
*           Test for sufficient growth in the norm of v.
*
VNORM = SASUM( N, VR, 1 )
IF( VNORM.GE.GROWTO*SCALE )
\$         GO TO 120
*
*           Choose new orthogonal starting vector and try again.
*
TEMP = EPS3 / ( ROOTN+ONE )
VR( 1 ) = EPS3
DO 100 I = 2, N
VR( I ) = TEMP
100       CONTINUE
VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
110    CONTINUE
*
*        Failure to find eigenvector in N iterations.
*
INFO = 1
*
120    CONTINUE
*
*        Normalize eigenvector.
*
I = ISAMAX( N, VR, 1 )
CALL SSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
ELSE
*
*        Complex eigenvalue.
*
IF( NOINIT ) THEN
*
*           Set initial vector.
*
DO 130 I = 1, N
VR( I ) = EPS3
VI( I ) = ZERO
130       CONTINUE
ELSE
*
*           Scale supplied initial vector.
*
NORM = SLAPY2( SNRM2( N, VR, 1 ), SNRM2( N, VI, 1 ) )
REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
CALL SSCAL( N, REC, VR, 1 )
CALL SSCAL( N, REC, VI, 1 )
END IF
*
IF( RIGHTV ) THEN
*
*           LU decomposition with partial pivoting of B, replacing zero
*           pivots by EPS3.
*
*           The imaginary part of the (i,j)-th element of U is stored in
*           B(j+1,i).
*
B( 2, 1 ) = -WI
DO 140 I = 2, N
B( I+1, 1 ) = ZERO
140       CONTINUE
*
DO 170 I = 1, N - 1
ABSBII = SLAPY2( B( I, I ), B( I+1, I ) )
EI = H( I+1, I )
IF( ABSBII.LT.ABS( EI ) ) THEN
*
*                 Interchange rows and eliminate.
*
XR = B( I, I ) / EI
XI = B( I+1, I ) / EI
B( I, I ) = EI
B( I+1, I ) = ZERO
DO 150 J = I + 1, N
TEMP = B( I+1, J )
B( I+1, J ) = B( I, J ) - XR*TEMP
B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP
B( I, J ) = TEMP
B( J+1, I ) = ZERO
150             CONTINUE
B( I+2, I ) = -WI
B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI
B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI
ELSE
*
*                 Eliminate without interchanging rows.
*
IF( ABSBII.EQ.ZERO ) THEN
B( I, I ) = EPS3
B( I+1, I ) = ZERO
ABSBII = EPS3
END IF
EI = ( EI / ABSBII ) / ABSBII
XR = B( I, I )*EI
XI = -B( I+1, I )*EI
DO 160 J = I + 1, N
B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) +
\$                             XI*B( J+1, I )
B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J )
160             CONTINUE
B( I+2, I+1 ) = B( I+2, I+1 ) - WI
END IF
*
*              Compute 1-norm of offdiagonal elements of i-th row.
*
WORK( I ) = SASUM( N-I, B( I, I+1 ), LDB ) +
\$                     SASUM( N-I, B( I+2, I ), 1 )
170       CONTINUE
IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO )
\$         B( N, N ) = EPS3
WORK( N ) = ZERO
*
I1 = N
I2 = 1
I3 = -1
ELSE
*
*           UL decomposition with partial pivoting of conjg(B),
*           replacing zero pivots by EPS3.
*
*           The imaginary part of the (i,j)-th element of U is stored in
*           B(j+1,i).
*
B( N+1, N ) = WI
DO 180 J = 1, N - 1
B( N+1, J ) = ZERO
180       CONTINUE
*
DO 210 J = N, 2, -1
EJ = H( J, J-1 )
ABSBJJ = SLAPY2( B( J, J ), B( J+1, J ) )
IF( ABSBJJ.LT.ABS( EJ ) ) THEN
*
*                 Interchange columns and eliminate
*
XR = B( J, J ) / EJ
XI = B( J+1, J ) / EJ
B( J, J ) = EJ
B( J+1, J ) = ZERO
DO 190 I = 1, J - 1
TEMP = B( I, J-1 )
B( I, J-1 ) = B( I, J ) - XR*TEMP
B( J, I ) = B( J+1, I ) - XI*TEMP
B( I, J ) = TEMP
B( J+1, I ) = ZERO
190             CONTINUE
B( J+1, J-1 ) = WI
B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI
B( J, J-1 ) = B( J, J-1 ) - XR*WI
ELSE
*
*                 Eliminate without interchange.
*
IF( ABSBJJ.EQ.ZERO ) THEN
B( J, J ) = EPS3
B( J+1, J ) = ZERO
ABSBJJ = EPS3
END IF
EJ = ( EJ / ABSBJJ ) / ABSBJJ
XR = B( J, J )*EJ
XI = -B( J+1, J )*EJ
DO 200 I = 1, J - 1
B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) +
\$                             XI*B( J+1, I )
B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J )
200             CONTINUE
B( J, J-1 ) = B( J, J-1 ) + WI
END IF
*
*              Compute 1-norm of offdiagonal elements of j-th column.
*
WORK( J ) = SASUM( J-1, B( 1, J ), 1 ) +
\$                     SASUM( J-1, B( J+1, 1 ), LDB )
210       CONTINUE
IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
\$         B( 1, 1 ) = EPS3
WORK( 1 ) = ZERO
*
I1 = 1
I2 = N
I3 = 1
END IF
*
DO 270 ITS = 1, N
SCALE = ONE
VMAX = ONE
VCRIT = BIGNUM
*
*           Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
*             or U'*(xr,xi) = scale*(vr,vi) for a left eigenvector,
*           overwriting (xr,xi) on (vr,vi).
*
DO 250 I = I1, I2, I3
*
IF( WORK( I ).GT.VCRIT ) THEN
REC = ONE / VMAX
CALL SSCAL( N, REC, VR, 1 )
CALL SSCAL( N, REC, VI, 1 )
SCALE = SCALE*REC
VMAX = ONE
VCRIT = BIGNUM
END IF
*
XR = VR( I )
XI = VI( I )
IF( RIGHTV ) THEN
DO 220 J = I + 1, N
XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J )
XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J )
220             CONTINUE
ELSE
DO 230 J = 1, I - 1
XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J )
XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J )
230             CONTINUE
END IF
*
W = ABS( B( I, I ) ) + ABS( B( I+1, I ) )
IF( W.GT.SMLNUM ) THEN
IF( W.LT.ONE ) THEN
W1 = ABS( XR ) + ABS( XI )
IF( W1.GT.W*BIGNUM ) THEN
REC = ONE / W1
CALL SSCAL( N, REC, VR, 1 )
CALL SSCAL( N, REC, VI, 1 )
XR = VR( I )
XI = VI( I )
SCALE = SCALE*REC
VMAX = VMAX*REC
END IF
END IF
*
*                 Divide by diagonal element of B.
*
CALL SLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ),
\$                         VI( I ) )
VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX )
VCRIT = BIGNUM / VMAX
ELSE
DO 240 J = 1, N
VR( J ) = ZERO
VI( J ) = ZERO
240             CONTINUE
VR( I ) = ONE
VI( I ) = ONE
SCALE = ZERO
VMAX = ONE
VCRIT = BIGNUM
END IF
250       CONTINUE
*
*           Test for sufficient growth in the norm of (VR,VI).
*
VNORM = SASUM( N, VR, 1 ) + SASUM( N, VI, 1 )
IF( VNORM.GE.GROWTO*SCALE )
\$         GO TO 280
*
*           Choose a new orthogonal starting vector and try again.
*
Y = EPS3 / ( ROOTN+ONE )
VR( 1 ) = EPS3
VI( 1 ) = ZERO
*
DO 260 I = 2, N
VR( I ) = Y
VI( I ) = ZERO
260       CONTINUE
VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
270    CONTINUE
*
*        Failure to find eigenvector in N iterations
*
INFO = 1
*
280    CONTINUE
*
*        Normalize eigenvector.
*
VNORM = ZERO
DO 290 I = 1, N
VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) )
290    CONTINUE
CALL SSCAL( N, ONE / VNORM, VR, 1 )
CALL SSCAL( N, ONE / VNORM, VI, 1 )
*
END IF
*
RETURN
*
*     End of SLAEIN
*
END

```