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