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