001:       SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
002:      \$                   EPS3, SMLNUM, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       LOGICAL            NOINIT, RIGHTV
011:       INTEGER            INFO, LDB, LDH, N
012:       DOUBLE PRECISION   EPS3, SMLNUM
013:       COMPLEX*16         W
014: *     ..
015: *     .. Array Arguments ..
016:       DOUBLE PRECISION   RWORK( * )
017:       COMPLEX*16         B( LDB, * ), H( LDH, * ), V( * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  ZLAEIN uses inverse iteration to find a right or left eigenvector
024: *  corresponding to the eigenvalue W of a complex upper Hessenberg
025: *  matrix H.
026: *
027: *  Arguments
028: *  =========
029: *
030: *  RIGHTV   (input) LOGICAL
031: *          = .TRUE. : compute right eigenvector;
032: *          = .FALSE.: compute left eigenvector.
033: *
034: *  NOINIT   (input) LOGICAL
035: *          = .TRUE. : no initial vector supplied in V
036: *          = .FALSE.: initial vector supplied in V.
037: *
038: *  N       (input) INTEGER
039: *          The order of the matrix H.  N >= 0.
040: *
041: *  H       (input) COMPLEX*16 array, dimension (LDH,N)
042: *          The upper Hessenberg matrix H.
043: *
044: *  LDH     (input) INTEGER
045: *          The leading dimension of the array H.  LDH >= max(1,N).
046: *
047: *  W       (input) COMPLEX*16
048: *          The eigenvalue of H whose corresponding right or left
049: *          eigenvector is to be computed.
050: *
051: *  V       (input/output) COMPLEX*16 array, dimension (N)
052: *          On entry, if NOINIT = .FALSE., V must contain a starting
053: *          vector for inverse iteration; otherwise V need not be set.
054: *          On exit, V contains the computed eigenvector, normalized so
055: *          that the component of largest magnitude has magnitude 1; here
056: *          the magnitude of a complex number (x,y) is taken to be
057: *          |x| + |y|.
058: *
059: *  B       (workspace) COMPLEX*16 array, dimension (LDB,N)
060: *
061: *  LDB     (input) INTEGER
062: *          The leading dimension of the array B.  LDB >= max(1,N).
063: *
064: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
065: *
066: *  EPS3    (input) DOUBLE PRECISION
067: *          A small machine-dependent value which is used to perturb
068: *          close eigenvalues, and to replace zero pivots.
069: *
070: *  SMLNUM  (input) DOUBLE PRECISION
071: *          A machine-dependent value close to the underflow threshold.
072: *
073: *  INFO    (output) INTEGER
074: *          = 0:  successful exit
075: *          = 1:  inverse iteration did not converge; V is set to the
076: *                last iterate.
077: *
078: *  =====================================================================
079: *
080: *     .. Parameters ..
081:       DOUBLE PRECISION   ONE, TENTH
082:       PARAMETER          ( ONE = 1.0D+0, TENTH = 1.0D-1 )
083:       COMPLEX*16         ZERO
084:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
085: *     ..
086: *     .. Local Scalars ..
087:       CHARACTER          NORMIN, TRANS
088:       INTEGER            I, IERR, ITS, J
089:       DOUBLE PRECISION   GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM
090:       COMPLEX*16         CDUM, EI, EJ, TEMP, X
091: *     ..
092: *     .. External Functions ..
093:       INTEGER            IZAMAX
094:       DOUBLE PRECISION   DZASUM, DZNRM2
096:       EXTERNAL           IZAMAX, DZASUM, DZNRM2, ZLADIV
097: *     ..
098: *     .. External Subroutines ..
099:       EXTERNAL           ZDSCAL, ZLATRS
100: *     ..
101: *     .. Intrinsic Functions ..
102:       INTRINSIC          ABS, DBLE, DIMAG, MAX, SQRT
103: *     ..
104: *     .. Statement Functions ..
105:       DOUBLE PRECISION   CABS1
106: *     ..
107: *     .. Statement Function definitions ..
108:       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
109: *     ..
110: *     .. Executable Statements ..
111: *
112:       INFO = 0
113: *
114: *     GROWTO is the threshold used in the acceptance test for an
115: *     eigenvector.
116: *
117:       ROOTN = SQRT( DBLE( N ) )
118:       GROWTO = TENTH / ROOTN
119:       NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
120: *
121: *     Form B = H - W*I (except that the subdiagonal elements are not
122: *     stored).
123: *
124:       DO 20 J = 1, N
125:          DO 10 I = 1, J - 1
126:             B( I, J ) = H( I, J )
127:    10    CONTINUE
128:          B( J, J ) = H( J, J ) - W
129:    20 CONTINUE
130: *
131:       IF( NOINIT ) THEN
132: *
133: *        Initialize V.
134: *
135:          DO 30 I = 1, N
136:             V( I ) = EPS3
137:    30    CONTINUE
138:       ELSE
139: *
140: *        Scale supplied initial vector.
141: *
142:          VNORM = DZNRM2( N, V, 1 )
143:          CALL ZDSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 )
144:       END IF
145: *
146:       IF( RIGHTV ) THEN
147: *
148: *        LU decomposition with partial pivoting of B, replacing zero
149: *        pivots by EPS3.
150: *
151:          DO 60 I = 1, N - 1
152:             EI = H( I+1, I )
153:             IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN
154: *
155: *              Interchange rows and eliminate.
156: *
157:                X = ZLADIV( B( I, I ), EI )
158:                B( I, I ) = EI
159:                DO 40 J = I + 1, N
160:                   TEMP = B( I+1, J )
161:                   B( I+1, J ) = B( I, J ) - X*TEMP
162:                   B( I, J ) = TEMP
163:    40          CONTINUE
164:             ELSE
165: *
166: *              Eliminate without interchange.
167: *
168:                IF( B( I, I ).EQ.ZERO )
169:      \$            B( I, I ) = EPS3
170:                X = ZLADIV( EI, B( I, I ) )
171:                IF( X.NE.ZERO ) THEN
172:                   DO 50 J = I + 1, N
173:                      B( I+1, J ) = B( I+1, J ) - X*B( I, J )
174:    50             CONTINUE
175:                END IF
176:             END IF
177:    60    CONTINUE
178:          IF( B( N, N ).EQ.ZERO )
179:      \$      B( N, N ) = EPS3
180: *
181:          TRANS = 'N'
182: *
183:       ELSE
184: *
185: *        UL decomposition with partial pivoting of B, replacing zero
186: *        pivots by EPS3.
187: *
188:          DO 90 J = N, 2, -1
189:             EJ = H( J, J-1 )
190:             IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN
191: *
192: *              Interchange columns and eliminate.
193: *
194:                X = ZLADIV( B( J, J ), EJ )
195:                B( J, J ) = EJ
196:                DO 70 I = 1, J - 1
197:                   TEMP = B( I, J-1 )
198:                   B( I, J-1 ) = B( I, J ) - X*TEMP
199:                   B( I, J ) = TEMP
200:    70          CONTINUE
201:             ELSE
202: *
203: *              Eliminate without interchange.
204: *
205:                IF( B( J, J ).EQ.ZERO )
206:      \$            B( J, J ) = EPS3
207:                X = ZLADIV( EJ, B( J, J ) )
208:                IF( X.NE.ZERO ) THEN
209:                   DO 80 I = 1, J - 1
210:                      B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
211:    80             CONTINUE
212:                END IF
213:             END IF
214:    90    CONTINUE
215:          IF( B( 1, 1 ).EQ.ZERO )
216:      \$      B( 1, 1 ) = EPS3
217: *
218:          TRANS = 'C'
219: *
220:       END IF
221: *
222:       NORMIN = 'N'
223:       DO 110 ITS = 1, N
224: *
225: *        Solve U*x = scale*v for a right eigenvector
226: *          or U'*x = scale*v for a left eigenvector,
227: *        overwriting x on v.
228: *
229:          CALL ZLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V,
230:      \$                SCALE, RWORK, IERR )
231:          NORMIN = 'Y'
232: *
233: *        Test for sufficient growth in the norm of v.
234: *
235:          VNORM = DZASUM( N, V, 1 )
236:          IF( VNORM.GE.GROWTO*SCALE )
237:      \$      GO TO 120
238: *
239: *        Choose new orthogonal starting vector and try again.
240: *
241:          RTEMP = EPS3 / ( ROOTN+ONE )
242:          V( 1 ) = EPS3
243:          DO 100 I = 2, N
244:             V( I ) = RTEMP
245:   100    CONTINUE
246:          V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN
247:   110 CONTINUE
248: *
249: *     Failure to find eigenvector in N iterations.
250: *
251:       INFO = 1
252: *
253:   120 CONTINUE
254: *
255: *     Normalize eigenvector.
256: *
257:       I = IZAMAX( N, V, 1 )
258:       CALL ZDSCAL( N, ONE / CABS1( V( I ) ), V, 1 )
259: *
260:       RETURN
261: *
262: *     End of ZLAEIN
263: *
264:       END
265: