1 SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
2 $ LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
3 $ IFAILR, INFO )
4 *
5 * -- LAPACK routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER EIGSRC, INITV, SIDE
12 INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
13 * ..
14 * .. Array Arguments ..
15 LOGICAL SELECT( * )
16 INTEGER IFAILL( * ), IFAILR( * )
17 DOUBLE PRECISION RWORK( * )
18 COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
19 $ W( * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * ZHSEIN uses inverse iteration to find specified right and/or left
26 * eigenvectors of a complex upper Hessenberg matrix H.
27 *
28 * The right eigenvector x and the left eigenvector y of the matrix H
29 * corresponding to an eigenvalue w are defined by:
30 *
31 * H * x = w * x, y**h * H = w * y**h
32 *
33 * where y**h denotes the conjugate transpose of the vector y.
34 *
35 * Arguments
36 * =========
37 *
38 * SIDE (input) CHARACTER*1
39 * = 'R': compute right eigenvectors only;
40 * = 'L': compute left eigenvectors only;
41 * = 'B': compute both right and left eigenvectors.
42 *
43 * EIGSRC (input) CHARACTER*1
44 * Specifies the source of eigenvalues supplied in W:
45 * = 'Q': the eigenvalues were found using ZHSEQR; thus, if
46 * H has zero subdiagonal elements, and so is
47 * block-triangular, then the j-th eigenvalue can be
48 * assumed to be an eigenvalue of the block containing
49 * the j-th row/column. This property allows ZHSEIN to
50 * perform inverse iteration on just one diagonal block.
51 * = 'N': no assumptions are made on the correspondence
52 * between eigenvalues and diagonal blocks. In this
53 * case, ZHSEIN must always perform inverse iteration
54 * using the whole matrix H.
55 *
56 * INITV (input) CHARACTER*1
57 * = 'N': no initial vectors are supplied;
58 * = 'U': user-supplied initial vectors are stored in the arrays
59 * VL and/or VR.
60 *
61 * SELECT (input) LOGICAL array, dimension (N)
62 * Specifies the eigenvectors to be computed. To select the
63 * eigenvector corresponding to the eigenvalue W(j),
64 * SELECT(j) must be set to .TRUE..
65 *
66 * N (input) INTEGER
67 * The order of the matrix H. N >= 0.
68 *
69 * H (input) COMPLEX*16 array, dimension (LDH,N)
70 * The upper Hessenberg matrix H.
71 *
72 * LDH (input) INTEGER
73 * The leading dimension of the array H. LDH >= max(1,N).
74 *
75 * W (input/output) COMPLEX*16 array, dimension (N)
76 * On entry, the eigenvalues of H.
77 * On exit, the real parts of W may have been altered since
78 * close eigenvalues are perturbed slightly in searching for
79 * independent eigenvectors.
80 *
81 * VL (input/output) COMPLEX*16 array, dimension (LDVL,MM)
82 * On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
83 * contain starting vectors for the inverse iteration for the
84 * left eigenvectors; the starting vector for each eigenvector
85 * must be in the same column in which the eigenvector will be
86 * stored.
87 * On exit, if SIDE = 'L' or 'B', the left eigenvectors
88 * specified by SELECT will be stored consecutively in the
89 * columns of VL, in the same order as their eigenvalues.
90 * If SIDE = 'R', VL is not referenced.
91 *
92 * LDVL (input) INTEGER
93 * The leading dimension of the array VL.
94 * LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
95 *
96 * VR (input/output) COMPLEX*16 array, dimension (LDVR,MM)
97 * On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
98 * contain starting vectors for the inverse iteration for the
99 * right eigenvectors; the starting vector for each eigenvector
100 * must be in the same column in which the eigenvector will be
101 * stored.
102 * On exit, if SIDE = 'R' or 'B', the right eigenvectors
103 * specified by SELECT will be stored consecutively in the
104 * columns of VR, in the same order as their eigenvalues.
105 * If SIDE = 'L', VR is not referenced.
106 *
107 * LDVR (input) INTEGER
108 * The leading dimension of the array VR.
109 * LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
110 *
111 * MM (input) INTEGER
112 * The number of columns in the arrays VL and/or VR. MM >= M.
113 *
114 * M (output) INTEGER
115 * The number of columns in the arrays VL and/or VR required to
116 * store the eigenvectors (= the number of .TRUE. elements in
117 * SELECT).
118 *
119 * WORK (workspace) COMPLEX*16 array, dimension (N*N)
120 *
121 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
122 *
123 * IFAILL (output) INTEGER array, dimension (MM)
124 * If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
125 * eigenvector in the i-th column of VL (corresponding to the
126 * eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
127 * eigenvector converged satisfactorily.
128 * If SIDE = 'R', IFAILL is not referenced.
129 *
130 * IFAILR (output) INTEGER array, dimension (MM)
131 * If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
132 * eigenvector in the i-th column of VR (corresponding to the
133 * eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
134 * eigenvector converged satisfactorily.
135 * If SIDE = 'L', IFAILR is not referenced.
136 *
137 * INFO (output) INTEGER
138 * = 0: successful exit
139 * < 0: if INFO = -i, the i-th argument had an illegal value
140 * > 0: if INFO = i, i is the number of eigenvectors which
141 * failed to converge; see IFAILL and IFAILR for further
142 * details.
143 *
144 * Further Details
145 * ===============
146 *
147 * Each eigenvector is normalized so that the element of largest
148 * magnitude has magnitude 1; here the magnitude of a complex number
149 * (x,y) is taken to be |x|+|y|.
150 *
151 * =====================================================================
152 *
153 * .. Parameters ..
154 COMPLEX*16 ZERO
155 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
156 DOUBLE PRECISION RZERO
157 PARAMETER ( RZERO = 0.0D+0 )
158 * ..
159 * .. Local Scalars ..
160 LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
161 INTEGER I, IINFO, K, KL, KLN, KR, KS, LDWORK
162 DOUBLE PRECISION EPS3, HNORM, SMLNUM, ULP, UNFL
163 COMPLEX*16 CDUM, WK
164 * ..
165 * .. External Functions ..
166 LOGICAL LSAME
167 DOUBLE PRECISION DLAMCH, ZLANHS
168 EXTERNAL LSAME, DLAMCH, ZLANHS
169 * ..
170 * .. External Subroutines ..
171 EXTERNAL XERBLA, ZLAEIN
172 * ..
173 * .. Intrinsic Functions ..
174 INTRINSIC ABS, DBLE, DIMAG, MAX
175 * ..
176 * .. Statement Functions ..
177 DOUBLE PRECISION CABS1
178 * ..
179 * .. Statement Function definitions ..
180 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
181 * ..
182 * .. Executable Statements ..
183 *
184 * Decode and test the input parameters.
185 *
186 BOTHV = LSAME( SIDE, 'B' )
187 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
188 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
189 *
190 FROMQR = LSAME( EIGSRC, 'Q' )
191 *
192 NOINIT = LSAME( INITV, 'N' )
193 *
194 * Set M to the number of columns required to store the selected
195 * eigenvectors.
196 *
197 M = 0
198 DO 10 K = 1, N
199 IF( SELECT( K ) )
200 $ M = M + 1
201 10 CONTINUE
202 *
203 INFO = 0
204 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
205 INFO = -1
206 ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
207 INFO = -2
208 ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
209 INFO = -3
210 ELSE IF( N.LT.0 ) THEN
211 INFO = -5
212 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
213 INFO = -7
214 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
215 INFO = -10
216 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
217 INFO = -12
218 ELSE IF( MM.LT.M ) THEN
219 INFO = -13
220 END IF
221 IF( INFO.NE.0 ) THEN
222 CALL XERBLA( 'ZHSEIN', -INFO )
223 RETURN
224 END IF
225 *
226 * Quick return if possible.
227 *
228 IF( N.EQ.0 )
229 $ RETURN
230 *
231 * Set machine-dependent constants.
232 *
233 UNFL = DLAMCH( 'Safe minimum' )
234 ULP = DLAMCH( 'Precision' )
235 SMLNUM = UNFL*( N / ULP )
236 *
237 LDWORK = N
238 *
239 KL = 1
240 KLN = 0
241 IF( FROMQR ) THEN
242 KR = 0
243 ELSE
244 KR = N
245 END IF
246 KS = 1
247 *
248 DO 100 K = 1, N
249 IF( SELECT( K ) ) THEN
250 *
251 * Compute eigenvector(s) corresponding to W(K).
252 *
253 IF( FROMQR ) THEN
254 *
255 * If affiliation of eigenvalues is known, check whether
256 * the matrix splits.
257 *
258 * Determine KL and KR such that 1 <= KL <= K <= KR <= N
259 * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
260 * KR = N).
261 *
262 * Then inverse iteration can be performed with the
263 * submatrix H(KL:N,KL:N) for a left eigenvector, and with
264 * the submatrix H(1:KR,1:KR) for a right eigenvector.
265 *
266 DO 20 I = K, KL + 1, -1
267 IF( H( I, I-1 ).EQ.ZERO )
268 $ GO TO 30
269 20 CONTINUE
270 30 CONTINUE
271 KL = I
272 IF( K.GT.KR ) THEN
273 DO 40 I = K, N - 1
274 IF( H( I+1, I ).EQ.ZERO )
275 $ GO TO 50
276 40 CONTINUE
277 50 CONTINUE
278 KR = I
279 END IF
280 END IF
281 *
282 IF( KL.NE.KLN ) THEN
283 KLN = KL
284 *
285 * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
286 * has not ben computed before.
287 *
288 HNORM = ZLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
289 IF( HNORM.GT.RZERO ) THEN
290 EPS3 = HNORM*ULP
291 ELSE
292 EPS3 = SMLNUM
293 END IF
294 END IF
295 *
296 * Perturb eigenvalue if it is close to any previous
297 * selected eigenvalues affiliated to the submatrix
298 * H(KL:KR,KL:KR). Close roots are modified by EPS3.
299 *
300 WK = W( K )
301 60 CONTINUE
302 DO 70 I = K - 1, KL, -1
303 IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN
304 WK = WK + EPS3
305 GO TO 60
306 END IF
307 70 CONTINUE
308 W( K ) = WK
309 *
310 IF( LEFTV ) THEN
311 *
312 * Compute left eigenvector.
313 *
314 CALL ZLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
315 $ WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3,
316 $ SMLNUM, IINFO )
317 IF( IINFO.GT.0 ) THEN
318 INFO = INFO + 1
319 IFAILL( KS ) = K
320 ELSE
321 IFAILL( KS ) = 0
322 END IF
323 DO 80 I = 1, KL - 1
324 VL( I, KS ) = ZERO
325 80 CONTINUE
326 END IF
327 IF( RIGHTV ) THEN
328 *
329 * Compute right eigenvector.
330 *
331 CALL ZLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ),
332 $ WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO )
333 IF( IINFO.GT.0 ) THEN
334 INFO = INFO + 1
335 IFAILR( KS ) = K
336 ELSE
337 IFAILR( KS ) = 0
338 END IF
339 DO 90 I = KR + 1, N
340 VR( I, KS ) = ZERO
341 90 CONTINUE
342 END IF
343 KS = KS + 1
344 END IF
345 100 CONTINUE
346 *
347 RETURN
348 *
349 * End of ZHSEIN
350 *
351 END
2 $ LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
3 $ IFAILR, INFO )
4 *
5 * -- LAPACK routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER EIGSRC, INITV, SIDE
12 INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
13 * ..
14 * .. Array Arguments ..
15 LOGICAL SELECT( * )
16 INTEGER IFAILL( * ), IFAILR( * )
17 DOUBLE PRECISION RWORK( * )
18 COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
19 $ W( * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * ZHSEIN uses inverse iteration to find specified right and/or left
26 * eigenvectors of a complex upper Hessenberg matrix H.
27 *
28 * The right eigenvector x and the left eigenvector y of the matrix H
29 * corresponding to an eigenvalue w are defined by:
30 *
31 * H * x = w * x, y**h * H = w * y**h
32 *
33 * where y**h denotes the conjugate transpose of the vector y.
34 *
35 * Arguments
36 * =========
37 *
38 * SIDE (input) CHARACTER*1
39 * = 'R': compute right eigenvectors only;
40 * = 'L': compute left eigenvectors only;
41 * = 'B': compute both right and left eigenvectors.
42 *
43 * EIGSRC (input) CHARACTER*1
44 * Specifies the source of eigenvalues supplied in W:
45 * = 'Q': the eigenvalues were found using ZHSEQR; thus, if
46 * H has zero subdiagonal elements, and so is
47 * block-triangular, then the j-th eigenvalue can be
48 * assumed to be an eigenvalue of the block containing
49 * the j-th row/column. This property allows ZHSEIN to
50 * perform inverse iteration on just one diagonal block.
51 * = 'N': no assumptions are made on the correspondence
52 * between eigenvalues and diagonal blocks. In this
53 * case, ZHSEIN must always perform inverse iteration
54 * using the whole matrix H.
55 *
56 * INITV (input) CHARACTER*1
57 * = 'N': no initial vectors are supplied;
58 * = 'U': user-supplied initial vectors are stored in the arrays
59 * VL and/or VR.
60 *
61 * SELECT (input) LOGICAL array, dimension (N)
62 * Specifies the eigenvectors to be computed. To select the
63 * eigenvector corresponding to the eigenvalue W(j),
64 * SELECT(j) must be set to .TRUE..
65 *
66 * N (input) INTEGER
67 * The order of the matrix H. N >= 0.
68 *
69 * H (input) COMPLEX*16 array, dimension (LDH,N)
70 * The upper Hessenberg matrix H.
71 *
72 * LDH (input) INTEGER
73 * The leading dimension of the array H. LDH >= max(1,N).
74 *
75 * W (input/output) COMPLEX*16 array, dimension (N)
76 * On entry, the eigenvalues of H.
77 * On exit, the real parts of W may have been altered since
78 * close eigenvalues are perturbed slightly in searching for
79 * independent eigenvectors.
80 *
81 * VL (input/output) COMPLEX*16 array, dimension (LDVL,MM)
82 * On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
83 * contain starting vectors for the inverse iteration for the
84 * left eigenvectors; the starting vector for each eigenvector
85 * must be in the same column in which the eigenvector will be
86 * stored.
87 * On exit, if SIDE = 'L' or 'B', the left eigenvectors
88 * specified by SELECT will be stored consecutively in the
89 * columns of VL, in the same order as their eigenvalues.
90 * If SIDE = 'R', VL is not referenced.
91 *
92 * LDVL (input) INTEGER
93 * The leading dimension of the array VL.
94 * LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
95 *
96 * VR (input/output) COMPLEX*16 array, dimension (LDVR,MM)
97 * On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
98 * contain starting vectors for the inverse iteration for the
99 * right eigenvectors; the starting vector for each eigenvector
100 * must be in the same column in which the eigenvector will be
101 * stored.
102 * On exit, if SIDE = 'R' or 'B', the right eigenvectors
103 * specified by SELECT will be stored consecutively in the
104 * columns of VR, in the same order as their eigenvalues.
105 * If SIDE = 'L', VR is not referenced.
106 *
107 * LDVR (input) INTEGER
108 * The leading dimension of the array VR.
109 * LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
110 *
111 * MM (input) INTEGER
112 * The number of columns in the arrays VL and/or VR. MM >= M.
113 *
114 * M (output) INTEGER
115 * The number of columns in the arrays VL and/or VR required to
116 * store the eigenvectors (= the number of .TRUE. elements in
117 * SELECT).
118 *
119 * WORK (workspace) COMPLEX*16 array, dimension (N*N)
120 *
121 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
122 *
123 * IFAILL (output) INTEGER array, dimension (MM)
124 * If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
125 * eigenvector in the i-th column of VL (corresponding to the
126 * eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
127 * eigenvector converged satisfactorily.
128 * If SIDE = 'R', IFAILL is not referenced.
129 *
130 * IFAILR (output) INTEGER array, dimension (MM)
131 * If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
132 * eigenvector in the i-th column of VR (corresponding to the
133 * eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
134 * eigenvector converged satisfactorily.
135 * If SIDE = 'L', IFAILR is not referenced.
136 *
137 * INFO (output) INTEGER
138 * = 0: successful exit
139 * < 0: if INFO = -i, the i-th argument had an illegal value
140 * > 0: if INFO = i, i is the number of eigenvectors which
141 * failed to converge; see IFAILL and IFAILR for further
142 * details.
143 *
144 * Further Details
145 * ===============
146 *
147 * Each eigenvector is normalized so that the element of largest
148 * magnitude has magnitude 1; here the magnitude of a complex number
149 * (x,y) is taken to be |x|+|y|.
150 *
151 * =====================================================================
152 *
153 * .. Parameters ..
154 COMPLEX*16 ZERO
155 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
156 DOUBLE PRECISION RZERO
157 PARAMETER ( RZERO = 0.0D+0 )
158 * ..
159 * .. Local Scalars ..
160 LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
161 INTEGER I, IINFO, K, KL, KLN, KR, KS, LDWORK
162 DOUBLE PRECISION EPS3, HNORM, SMLNUM, ULP, UNFL
163 COMPLEX*16 CDUM, WK
164 * ..
165 * .. External Functions ..
166 LOGICAL LSAME
167 DOUBLE PRECISION DLAMCH, ZLANHS
168 EXTERNAL LSAME, DLAMCH, ZLANHS
169 * ..
170 * .. External Subroutines ..
171 EXTERNAL XERBLA, ZLAEIN
172 * ..
173 * .. Intrinsic Functions ..
174 INTRINSIC ABS, DBLE, DIMAG, MAX
175 * ..
176 * .. Statement Functions ..
177 DOUBLE PRECISION CABS1
178 * ..
179 * .. Statement Function definitions ..
180 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
181 * ..
182 * .. Executable Statements ..
183 *
184 * Decode and test the input parameters.
185 *
186 BOTHV = LSAME( SIDE, 'B' )
187 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
188 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
189 *
190 FROMQR = LSAME( EIGSRC, 'Q' )
191 *
192 NOINIT = LSAME( INITV, 'N' )
193 *
194 * Set M to the number of columns required to store the selected
195 * eigenvectors.
196 *
197 M = 0
198 DO 10 K = 1, N
199 IF( SELECT( K ) )
200 $ M = M + 1
201 10 CONTINUE
202 *
203 INFO = 0
204 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
205 INFO = -1
206 ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
207 INFO = -2
208 ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
209 INFO = -3
210 ELSE IF( N.LT.0 ) THEN
211 INFO = -5
212 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
213 INFO = -7
214 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
215 INFO = -10
216 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
217 INFO = -12
218 ELSE IF( MM.LT.M ) THEN
219 INFO = -13
220 END IF
221 IF( INFO.NE.0 ) THEN
222 CALL XERBLA( 'ZHSEIN', -INFO )
223 RETURN
224 END IF
225 *
226 * Quick return if possible.
227 *
228 IF( N.EQ.0 )
229 $ RETURN
230 *
231 * Set machine-dependent constants.
232 *
233 UNFL = DLAMCH( 'Safe minimum' )
234 ULP = DLAMCH( 'Precision' )
235 SMLNUM = UNFL*( N / ULP )
236 *
237 LDWORK = N
238 *
239 KL = 1
240 KLN = 0
241 IF( FROMQR ) THEN
242 KR = 0
243 ELSE
244 KR = N
245 END IF
246 KS = 1
247 *
248 DO 100 K = 1, N
249 IF( SELECT( K ) ) THEN
250 *
251 * Compute eigenvector(s) corresponding to W(K).
252 *
253 IF( FROMQR ) THEN
254 *
255 * If affiliation of eigenvalues is known, check whether
256 * the matrix splits.
257 *
258 * Determine KL and KR such that 1 <= KL <= K <= KR <= N
259 * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
260 * KR = N).
261 *
262 * Then inverse iteration can be performed with the
263 * submatrix H(KL:N,KL:N) for a left eigenvector, and with
264 * the submatrix H(1:KR,1:KR) for a right eigenvector.
265 *
266 DO 20 I = K, KL + 1, -1
267 IF( H( I, I-1 ).EQ.ZERO )
268 $ GO TO 30
269 20 CONTINUE
270 30 CONTINUE
271 KL = I
272 IF( K.GT.KR ) THEN
273 DO 40 I = K, N - 1
274 IF( H( I+1, I ).EQ.ZERO )
275 $ GO TO 50
276 40 CONTINUE
277 50 CONTINUE
278 KR = I
279 END IF
280 END IF
281 *
282 IF( KL.NE.KLN ) THEN
283 KLN = KL
284 *
285 * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
286 * has not ben computed before.
287 *
288 HNORM = ZLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
289 IF( HNORM.GT.RZERO ) THEN
290 EPS3 = HNORM*ULP
291 ELSE
292 EPS3 = SMLNUM
293 END IF
294 END IF
295 *
296 * Perturb eigenvalue if it is close to any previous
297 * selected eigenvalues affiliated to the submatrix
298 * H(KL:KR,KL:KR). Close roots are modified by EPS3.
299 *
300 WK = W( K )
301 60 CONTINUE
302 DO 70 I = K - 1, KL, -1
303 IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN
304 WK = WK + EPS3
305 GO TO 60
306 END IF
307 70 CONTINUE
308 W( K ) = WK
309 *
310 IF( LEFTV ) THEN
311 *
312 * Compute left eigenvector.
313 *
314 CALL ZLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
315 $ WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3,
316 $ SMLNUM, IINFO )
317 IF( IINFO.GT.0 ) THEN
318 INFO = INFO + 1
319 IFAILL( KS ) = K
320 ELSE
321 IFAILL( KS ) = 0
322 END IF
323 DO 80 I = 1, KL - 1
324 VL( I, KS ) = ZERO
325 80 CONTINUE
326 END IF
327 IF( RIGHTV ) THEN
328 *
329 * Compute right eigenvector.
330 *
331 CALL ZLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ),
332 $ WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO )
333 IF( IINFO.GT.0 ) THEN
334 INFO = INFO + 1
335 IFAILR( KS ) = K
336 ELSE
337 IFAILR( KS ) = 0
338 END IF
339 DO 90 I = KR + 1, N
340 VR( I, KS ) = ZERO
341 90 CONTINUE
342 END IF
343 KS = KS + 1
344 END IF
345 100 CONTINUE
346 *
347 RETURN
348 *
349 * End of ZHSEIN
350 *
351 END