1 SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
2 $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 LOGICAL NOINIT, RIGHTV
11 INTEGER INFO, LDB, LDH, N
12 DOUBLE PRECISION BIGNUM, EPS3, SMLNUM, WI, WR
13 * ..
14 * .. Array Arguments ..
15 DOUBLE PRECISION B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
16 $ WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DLAEIN uses inverse iteration to find a right or left eigenvector
23 * corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
24 * matrix H.
25 *
26 * Arguments
27 * =========
28 *
29 * RIGHTV (input) LOGICAL
30 * = .TRUE. : compute right eigenvector;
31 * = .FALSE.: compute left eigenvector.
32 *
33 * NOINIT (input) LOGICAL
34 * = .TRUE. : no initial vector supplied in (VR,VI).
35 * = .FALSE.: initial vector supplied in (VR,VI).
36 *
37 * N (input) INTEGER
38 * The order of the matrix H. N >= 0.
39 *
40 * H (input) DOUBLE PRECISION array, dimension (LDH,N)
41 * The upper Hessenberg matrix H.
42 *
43 * LDH (input) INTEGER
44 * The leading dimension of the array H. LDH >= max(1,N).
45 *
46 * WR (input) DOUBLE PRECISION
47 * WI (input) DOUBLE PRECISION
48 * The real and imaginary parts of the eigenvalue of H whose
49 * corresponding right or left eigenvector is to be computed.
50 *
51 * VR (input/output) DOUBLE PRECISION array, dimension (N)
52 * VI (input/output) DOUBLE PRECISION array, dimension (N)
53 * On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
54 * a real starting vector for inverse iteration using the real
55 * eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
56 * must contain the real and imaginary parts of a complex
57 * starting vector for inverse iteration using the complex
58 * eigenvalue (WR,WI); otherwise VR and VI need not be set.
59 * On exit, if WI = 0.0 (real eigenvalue), VR contains the
60 * computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
61 * VR and VI contain the real and imaginary parts of the
62 * computed complex eigenvector. The eigenvector is normalized
63 * so that the component of largest magnitude has magnitude 1;
64 * here the magnitude of a complex number (x,y) is taken to be
65 * |x| + |y|.
66 * VI is not referenced if WI = 0.0.
67 *
68 * B (workspace) DOUBLE PRECISION array, dimension (LDB,N)
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of the array B. LDB >= N+1.
72 *
73 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
74 *
75 * EPS3 (input) DOUBLE PRECISION
76 * A small machine-dependent value which is used to perturb
77 * close eigenvalues, and to replace zero pivots.
78 *
79 * SMLNUM (input) DOUBLE PRECISION
80 * A machine-dependent value close to the underflow threshold.
81 *
82 * BIGNUM (input) DOUBLE PRECISION
83 * A machine-dependent value close to the overflow threshold.
84 *
85 * INFO (output) INTEGER
86 * = 0: successful exit
87 * = 1: inverse iteration did not converge; VR is set to the
88 * last iterate, and so is VI if WI.ne.0.0.
89 *
90 * =====================================================================
91 *
92 * .. Parameters ..
93 DOUBLE PRECISION ZERO, ONE, TENTH
94 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 )
95 * ..
96 * .. Local Scalars ..
97 CHARACTER NORMIN, TRANS
98 INTEGER I, I1, I2, I3, IERR, ITS, J
99 DOUBLE PRECISION ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
100 $ REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W,
101 $ W1, X, XI, XR, Y
102 * ..
103 * .. External Functions ..
104 INTEGER IDAMAX
105 DOUBLE PRECISION DASUM, DLAPY2, DNRM2
106 EXTERNAL IDAMAX, DASUM, DLAPY2, DNRM2
107 * ..
108 * .. External Subroutines ..
109 EXTERNAL DLADIV, DLATRS, DSCAL
110 * ..
111 * .. Intrinsic Functions ..
112 INTRINSIC ABS, DBLE, MAX, SQRT
113 * ..
114 * .. Executable Statements ..
115 *
116 INFO = 0
117 *
118 * GROWTO is the threshold used in the acceptance test for an
119 * eigenvector.
120 *
121 ROOTN = SQRT( DBLE( N ) )
122 GROWTO = TENTH / ROOTN
123 NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
124 *
125 * Form B = H - (WR,WI)*I (except that the subdiagonal elements and
126 * the imaginary parts of the diagonal elements are not stored).
127 *
128 DO 20 J = 1, N
129 DO 10 I = 1, J - 1
130 B( I, J ) = H( I, J )
131 10 CONTINUE
132 B( J, J ) = H( J, J ) - WR
133 20 CONTINUE
134 *
135 IF( WI.EQ.ZERO ) THEN
136 *
137 * Real eigenvalue.
138 *
139 IF( NOINIT ) THEN
140 *
141 * Set initial vector.
142 *
143 DO 30 I = 1, N
144 VR( I ) = EPS3
145 30 CONTINUE
146 ELSE
147 *
148 * Scale supplied initial vector.
149 *
150 VNORM = DNRM2( N, VR, 1 )
151 CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR,
152 $ 1 )
153 END IF
154 *
155 IF( RIGHTV ) THEN
156 *
157 * LU decomposition with partial pivoting of B, replacing zero
158 * pivots by EPS3.
159 *
160 DO 60 I = 1, N - 1
161 EI = H( I+1, I )
162 IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN
163 *
164 * Interchange rows and eliminate.
165 *
166 X = B( I, I ) / EI
167 B( I, I ) = EI
168 DO 40 J = I + 1, N
169 TEMP = B( I+1, J )
170 B( I+1, J ) = B( I, J ) - X*TEMP
171 B( I, J ) = TEMP
172 40 CONTINUE
173 ELSE
174 *
175 * Eliminate without interchange.
176 *
177 IF( B( I, I ).EQ.ZERO )
178 $ B( I, I ) = EPS3
179 X = EI / B( I, I )
180 IF( X.NE.ZERO ) THEN
181 DO 50 J = I + 1, N
182 B( I+1, J ) = B( I+1, J ) - X*B( I, J )
183 50 CONTINUE
184 END IF
185 END IF
186 60 CONTINUE
187 IF( B( N, N ).EQ.ZERO )
188 $ B( N, N ) = EPS3
189 *
190 TRANS = 'N'
191 *
192 ELSE
193 *
194 * UL decomposition with partial pivoting of B, replacing zero
195 * pivots by EPS3.
196 *
197 DO 90 J = N, 2, -1
198 EJ = H( J, J-1 )
199 IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
200 *
201 * Interchange columns and eliminate.
202 *
203 X = B( J, J ) / EJ
204 B( J, J ) = EJ
205 DO 70 I = 1, J - 1
206 TEMP = B( I, J-1 )
207 B( I, J-1 ) = B( I, J ) - X*TEMP
208 B( I, J ) = TEMP
209 70 CONTINUE
210 ELSE
211 *
212 * Eliminate without interchange.
213 *
214 IF( B( J, J ).EQ.ZERO )
215 $ B( J, J ) = EPS3
216 X = EJ / B( J, J )
217 IF( X.NE.ZERO ) THEN
218 DO 80 I = 1, J - 1
219 B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
220 80 CONTINUE
221 END IF
222 END IF
223 90 CONTINUE
224 IF( B( 1, 1 ).EQ.ZERO )
225 $ B( 1, 1 ) = EPS3
226 *
227 TRANS = 'T'
228 *
229 END IF
230 *
231 NORMIN = 'N'
232 DO 110 ITS = 1, N
233 *
234 * Solve U*x = scale*v for a right eigenvector
235 * or U**T*x = scale*v for a left eigenvector,
236 * overwriting x on v.
237 *
238 CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB,
239 $ VR, SCALE, WORK, IERR )
240 NORMIN = 'Y'
241 *
242 * Test for sufficient growth in the norm of v.
243 *
244 VNORM = DASUM( N, VR, 1 )
245 IF( VNORM.GE.GROWTO*SCALE )
246 $ GO TO 120
247 *
248 * Choose new orthogonal starting vector and try again.
249 *
250 TEMP = EPS3 / ( ROOTN+ONE )
251 VR( 1 ) = EPS3
252 DO 100 I = 2, N
253 VR( I ) = TEMP
254 100 CONTINUE
255 VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
256 110 CONTINUE
257 *
258 * Failure to find eigenvector in N iterations.
259 *
260 INFO = 1
261 *
262 120 CONTINUE
263 *
264 * Normalize eigenvector.
265 *
266 I = IDAMAX( N, VR, 1 )
267 CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
268 ELSE
269 *
270 * Complex eigenvalue.
271 *
272 IF( NOINIT ) THEN
273 *
274 * Set initial vector.
275 *
276 DO 130 I = 1, N
277 VR( I ) = EPS3
278 VI( I ) = ZERO
279 130 CONTINUE
280 ELSE
281 *
282 * Scale supplied initial vector.
283 *
284 NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) )
285 REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
286 CALL DSCAL( N, REC, VR, 1 )
287 CALL DSCAL( N, REC, VI, 1 )
288 END IF
289 *
290 IF( RIGHTV ) THEN
291 *
292 * LU decomposition with partial pivoting of B, replacing zero
293 * pivots by EPS3.
294 *
295 * The imaginary part of the (i,j)-th element of U is stored in
296 * B(j+1,i).
297 *
298 B( 2, 1 ) = -WI
299 DO 140 I = 2, N
300 B( I+1, 1 ) = ZERO
301 140 CONTINUE
302 *
303 DO 170 I = 1, N - 1
304 ABSBII = DLAPY2( B( I, I ), B( I+1, I ) )
305 EI = H( I+1, I )
306 IF( ABSBII.LT.ABS( EI ) ) THEN
307 *
308 * Interchange rows and eliminate.
309 *
310 XR = B( I, I ) / EI
311 XI = B( I+1, I ) / EI
312 B( I, I ) = EI
313 B( I+1, I ) = ZERO
314 DO 150 J = I + 1, N
315 TEMP = B( I+1, J )
316 B( I+1, J ) = B( I, J ) - XR*TEMP
317 B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP
318 B( I, J ) = TEMP
319 B( J+1, I ) = ZERO
320 150 CONTINUE
321 B( I+2, I ) = -WI
322 B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI
323 B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI
324 ELSE
325 *
326 * Eliminate without interchanging rows.
327 *
328 IF( ABSBII.EQ.ZERO ) THEN
329 B( I, I ) = EPS3
330 B( I+1, I ) = ZERO
331 ABSBII = EPS3
332 END IF
333 EI = ( EI / ABSBII ) / ABSBII
334 XR = B( I, I )*EI
335 XI = -B( I+1, I )*EI
336 DO 160 J = I + 1, N
337 B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) +
338 $ XI*B( J+1, I )
339 B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J )
340 160 CONTINUE
341 B( I+2, I+1 ) = B( I+2, I+1 ) - WI
342 END IF
343 *
344 * Compute 1-norm of offdiagonal elements of i-th row.
345 *
346 WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) +
347 $ DASUM( N-I, B( I+2, I ), 1 )
348 170 CONTINUE
349 IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO )
350 $ B( N, N ) = EPS3
351 WORK( N ) = ZERO
352 *
353 I1 = N
354 I2 = 1
355 I3 = -1
356 ELSE
357 *
358 * UL decomposition with partial pivoting of conjg(B),
359 * replacing zero pivots by EPS3.
360 *
361 * The imaginary part of the (i,j)-th element of U is stored in
362 * B(j+1,i).
363 *
364 B( N+1, N ) = WI
365 DO 180 J = 1, N - 1
366 B( N+1, J ) = ZERO
367 180 CONTINUE
368 *
369 DO 210 J = N, 2, -1
370 EJ = H( J, J-1 )
371 ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) )
372 IF( ABSBJJ.LT.ABS( EJ ) ) THEN
373 *
374 * Interchange columns and eliminate
375 *
376 XR = B( J, J ) / EJ
377 XI = B( J+1, J ) / EJ
378 B( J, J ) = EJ
379 B( J+1, J ) = ZERO
380 DO 190 I = 1, J - 1
381 TEMP = B( I, J-1 )
382 B( I, J-1 ) = B( I, J ) - XR*TEMP
383 B( J, I ) = B( J+1, I ) - XI*TEMP
384 B( I, J ) = TEMP
385 B( J+1, I ) = ZERO
386 190 CONTINUE
387 B( J+1, J-1 ) = WI
388 B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI
389 B( J, J-1 ) = B( J, J-1 ) - XR*WI
390 ELSE
391 *
392 * Eliminate without interchange.
393 *
394 IF( ABSBJJ.EQ.ZERO ) THEN
395 B( J, J ) = EPS3
396 B( J+1, J ) = ZERO
397 ABSBJJ = EPS3
398 END IF
399 EJ = ( EJ / ABSBJJ ) / ABSBJJ
400 XR = B( J, J )*EJ
401 XI = -B( J+1, J )*EJ
402 DO 200 I = 1, J - 1
403 B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) +
404 $ XI*B( J+1, I )
405 B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J )
406 200 CONTINUE
407 B( J, J-1 ) = B( J, J-1 ) + WI
408 END IF
409 *
410 * Compute 1-norm of offdiagonal elements of j-th column.
411 *
412 WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) +
413 $ DASUM( J-1, B( J+1, 1 ), LDB )
414 210 CONTINUE
415 IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
416 $ B( 1, 1 ) = EPS3
417 WORK( 1 ) = ZERO
418 *
419 I1 = 1
420 I2 = N
421 I3 = 1
422 END IF
423 *
424 DO 270 ITS = 1, N
425 SCALE = ONE
426 VMAX = ONE
427 VCRIT = BIGNUM
428 *
429 * Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
430 * or U**T*(xr,xi) = scale*(vr,vi) for a left eigenvector,
431 * overwriting (xr,xi) on (vr,vi).
432 *
433 DO 250 I = I1, I2, I3
434 *
435 IF( WORK( I ).GT.VCRIT ) THEN
436 REC = ONE / VMAX
437 CALL DSCAL( N, REC, VR, 1 )
438 CALL DSCAL( N, REC, VI, 1 )
439 SCALE = SCALE*REC
440 VMAX = ONE
441 VCRIT = BIGNUM
442 END IF
443 *
444 XR = VR( I )
445 XI = VI( I )
446 IF( RIGHTV ) THEN
447 DO 220 J = I + 1, N
448 XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J )
449 XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J )
450 220 CONTINUE
451 ELSE
452 DO 230 J = 1, I - 1
453 XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J )
454 XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J )
455 230 CONTINUE
456 END IF
457 *
458 W = ABS( B( I, I ) ) + ABS( B( I+1, I ) )
459 IF( W.GT.SMLNUM ) THEN
460 IF( W.LT.ONE ) THEN
461 W1 = ABS( XR ) + ABS( XI )
462 IF( W1.GT.W*BIGNUM ) THEN
463 REC = ONE / W1
464 CALL DSCAL( N, REC, VR, 1 )
465 CALL DSCAL( N, REC, VI, 1 )
466 XR = VR( I )
467 XI = VI( I )
468 SCALE = SCALE*REC
469 VMAX = VMAX*REC
470 END IF
471 END IF
472 *
473 * Divide by diagonal element of B.
474 *
475 CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ),
476 $ VI( I ) )
477 VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX )
478 VCRIT = BIGNUM / VMAX
479 ELSE
480 DO 240 J = 1, N
481 VR( J ) = ZERO
482 VI( J ) = ZERO
483 240 CONTINUE
484 VR( I ) = ONE
485 VI( I ) = ONE
486 SCALE = ZERO
487 VMAX = ONE
488 VCRIT = BIGNUM
489 END IF
490 250 CONTINUE
491 *
492 * Test for sufficient growth in the norm of (VR,VI).
493 *
494 VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 )
495 IF( VNORM.GE.GROWTO*SCALE )
496 $ GO TO 280
497 *
498 * Choose a new orthogonal starting vector and try again.
499 *
500 Y = EPS3 / ( ROOTN+ONE )
501 VR( 1 ) = EPS3
502 VI( 1 ) = ZERO
503 *
504 DO 260 I = 2, N
505 VR( I ) = Y
506 VI( I ) = ZERO
507 260 CONTINUE
508 VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
509 270 CONTINUE
510 *
511 * Failure to find eigenvector in N iterations
512 *
513 INFO = 1
514 *
515 280 CONTINUE
516 *
517 * Normalize eigenvector.
518 *
519 VNORM = ZERO
520 DO 290 I = 1, N
521 VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) )
522 290 CONTINUE
523 CALL DSCAL( N, ONE / VNORM, VR, 1 )
524 CALL DSCAL( N, ONE / VNORM, VI, 1 )
525 *
526 END IF
527 *
528 RETURN
529 *
530 * End of DLAEIN
531 *
532 END
2 $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 LOGICAL NOINIT, RIGHTV
11 INTEGER INFO, LDB, LDH, N
12 DOUBLE PRECISION BIGNUM, EPS3, SMLNUM, WI, WR
13 * ..
14 * .. Array Arguments ..
15 DOUBLE PRECISION B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
16 $ WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DLAEIN uses inverse iteration to find a right or left eigenvector
23 * corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
24 * matrix H.
25 *
26 * Arguments
27 * =========
28 *
29 * RIGHTV (input) LOGICAL
30 * = .TRUE. : compute right eigenvector;
31 * = .FALSE.: compute left eigenvector.
32 *
33 * NOINIT (input) LOGICAL
34 * = .TRUE. : no initial vector supplied in (VR,VI).
35 * = .FALSE.: initial vector supplied in (VR,VI).
36 *
37 * N (input) INTEGER
38 * The order of the matrix H. N >= 0.
39 *
40 * H (input) DOUBLE PRECISION array, dimension (LDH,N)
41 * The upper Hessenberg matrix H.
42 *
43 * LDH (input) INTEGER
44 * The leading dimension of the array H. LDH >= max(1,N).
45 *
46 * WR (input) DOUBLE PRECISION
47 * WI (input) DOUBLE PRECISION
48 * The real and imaginary parts of the eigenvalue of H whose
49 * corresponding right or left eigenvector is to be computed.
50 *
51 * VR (input/output) DOUBLE PRECISION array, dimension (N)
52 * VI (input/output) DOUBLE PRECISION array, dimension (N)
53 * On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
54 * a real starting vector for inverse iteration using the real
55 * eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
56 * must contain the real and imaginary parts of a complex
57 * starting vector for inverse iteration using the complex
58 * eigenvalue (WR,WI); otherwise VR and VI need not be set.
59 * On exit, if WI = 0.0 (real eigenvalue), VR contains the
60 * computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
61 * VR and VI contain the real and imaginary parts of the
62 * computed complex eigenvector. The eigenvector is normalized
63 * so that the component of largest magnitude has magnitude 1;
64 * here the magnitude of a complex number (x,y) is taken to be
65 * |x| + |y|.
66 * VI is not referenced if WI = 0.0.
67 *
68 * B (workspace) DOUBLE PRECISION array, dimension (LDB,N)
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of the array B. LDB >= N+1.
72 *
73 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
74 *
75 * EPS3 (input) DOUBLE PRECISION
76 * A small machine-dependent value which is used to perturb
77 * close eigenvalues, and to replace zero pivots.
78 *
79 * SMLNUM (input) DOUBLE PRECISION
80 * A machine-dependent value close to the underflow threshold.
81 *
82 * BIGNUM (input) DOUBLE PRECISION
83 * A machine-dependent value close to the overflow threshold.
84 *
85 * INFO (output) INTEGER
86 * = 0: successful exit
87 * = 1: inverse iteration did not converge; VR is set to the
88 * last iterate, and so is VI if WI.ne.0.0.
89 *
90 * =====================================================================
91 *
92 * .. Parameters ..
93 DOUBLE PRECISION ZERO, ONE, TENTH
94 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 )
95 * ..
96 * .. Local Scalars ..
97 CHARACTER NORMIN, TRANS
98 INTEGER I, I1, I2, I3, IERR, ITS, J
99 DOUBLE PRECISION ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
100 $ REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W,
101 $ W1, X, XI, XR, Y
102 * ..
103 * .. External Functions ..
104 INTEGER IDAMAX
105 DOUBLE PRECISION DASUM, DLAPY2, DNRM2
106 EXTERNAL IDAMAX, DASUM, DLAPY2, DNRM2
107 * ..
108 * .. External Subroutines ..
109 EXTERNAL DLADIV, DLATRS, DSCAL
110 * ..
111 * .. Intrinsic Functions ..
112 INTRINSIC ABS, DBLE, MAX, SQRT
113 * ..
114 * .. Executable Statements ..
115 *
116 INFO = 0
117 *
118 * GROWTO is the threshold used in the acceptance test for an
119 * eigenvector.
120 *
121 ROOTN = SQRT( DBLE( N ) )
122 GROWTO = TENTH / ROOTN
123 NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
124 *
125 * Form B = H - (WR,WI)*I (except that the subdiagonal elements and
126 * the imaginary parts of the diagonal elements are not stored).
127 *
128 DO 20 J = 1, N
129 DO 10 I = 1, J - 1
130 B( I, J ) = H( I, J )
131 10 CONTINUE
132 B( J, J ) = H( J, J ) - WR
133 20 CONTINUE
134 *
135 IF( WI.EQ.ZERO ) THEN
136 *
137 * Real eigenvalue.
138 *
139 IF( NOINIT ) THEN
140 *
141 * Set initial vector.
142 *
143 DO 30 I = 1, N
144 VR( I ) = EPS3
145 30 CONTINUE
146 ELSE
147 *
148 * Scale supplied initial vector.
149 *
150 VNORM = DNRM2( N, VR, 1 )
151 CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR,
152 $ 1 )
153 END IF
154 *
155 IF( RIGHTV ) THEN
156 *
157 * LU decomposition with partial pivoting of B, replacing zero
158 * pivots by EPS3.
159 *
160 DO 60 I = 1, N - 1
161 EI = H( I+1, I )
162 IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN
163 *
164 * Interchange rows and eliminate.
165 *
166 X = B( I, I ) / EI
167 B( I, I ) = EI
168 DO 40 J = I + 1, N
169 TEMP = B( I+1, J )
170 B( I+1, J ) = B( I, J ) - X*TEMP
171 B( I, J ) = TEMP
172 40 CONTINUE
173 ELSE
174 *
175 * Eliminate without interchange.
176 *
177 IF( B( I, I ).EQ.ZERO )
178 $ B( I, I ) = EPS3
179 X = EI / B( I, I )
180 IF( X.NE.ZERO ) THEN
181 DO 50 J = I + 1, N
182 B( I+1, J ) = B( I+1, J ) - X*B( I, J )
183 50 CONTINUE
184 END IF
185 END IF
186 60 CONTINUE
187 IF( B( N, N ).EQ.ZERO )
188 $ B( N, N ) = EPS3
189 *
190 TRANS = 'N'
191 *
192 ELSE
193 *
194 * UL decomposition with partial pivoting of B, replacing zero
195 * pivots by EPS3.
196 *
197 DO 90 J = N, 2, -1
198 EJ = H( J, J-1 )
199 IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
200 *
201 * Interchange columns and eliminate.
202 *
203 X = B( J, J ) / EJ
204 B( J, J ) = EJ
205 DO 70 I = 1, J - 1
206 TEMP = B( I, J-1 )
207 B( I, J-1 ) = B( I, J ) - X*TEMP
208 B( I, J ) = TEMP
209 70 CONTINUE
210 ELSE
211 *
212 * Eliminate without interchange.
213 *
214 IF( B( J, J ).EQ.ZERO )
215 $ B( J, J ) = EPS3
216 X = EJ / B( J, J )
217 IF( X.NE.ZERO ) THEN
218 DO 80 I = 1, J - 1
219 B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
220 80 CONTINUE
221 END IF
222 END IF
223 90 CONTINUE
224 IF( B( 1, 1 ).EQ.ZERO )
225 $ B( 1, 1 ) = EPS3
226 *
227 TRANS = 'T'
228 *
229 END IF
230 *
231 NORMIN = 'N'
232 DO 110 ITS = 1, N
233 *
234 * Solve U*x = scale*v for a right eigenvector
235 * or U**T*x = scale*v for a left eigenvector,
236 * overwriting x on v.
237 *
238 CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB,
239 $ VR, SCALE, WORK, IERR )
240 NORMIN = 'Y'
241 *
242 * Test for sufficient growth in the norm of v.
243 *
244 VNORM = DASUM( N, VR, 1 )
245 IF( VNORM.GE.GROWTO*SCALE )
246 $ GO TO 120
247 *
248 * Choose new orthogonal starting vector and try again.
249 *
250 TEMP = EPS3 / ( ROOTN+ONE )
251 VR( 1 ) = EPS3
252 DO 100 I = 2, N
253 VR( I ) = TEMP
254 100 CONTINUE
255 VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
256 110 CONTINUE
257 *
258 * Failure to find eigenvector in N iterations.
259 *
260 INFO = 1
261 *
262 120 CONTINUE
263 *
264 * Normalize eigenvector.
265 *
266 I = IDAMAX( N, VR, 1 )
267 CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
268 ELSE
269 *
270 * Complex eigenvalue.
271 *
272 IF( NOINIT ) THEN
273 *
274 * Set initial vector.
275 *
276 DO 130 I = 1, N
277 VR( I ) = EPS3
278 VI( I ) = ZERO
279 130 CONTINUE
280 ELSE
281 *
282 * Scale supplied initial vector.
283 *
284 NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) )
285 REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
286 CALL DSCAL( N, REC, VR, 1 )
287 CALL DSCAL( N, REC, VI, 1 )
288 END IF
289 *
290 IF( RIGHTV ) THEN
291 *
292 * LU decomposition with partial pivoting of B, replacing zero
293 * pivots by EPS3.
294 *
295 * The imaginary part of the (i,j)-th element of U is stored in
296 * B(j+1,i).
297 *
298 B( 2, 1 ) = -WI
299 DO 140 I = 2, N
300 B( I+1, 1 ) = ZERO
301 140 CONTINUE
302 *
303 DO 170 I = 1, N - 1
304 ABSBII = DLAPY2( B( I, I ), B( I+1, I ) )
305 EI = H( I+1, I )
306 IF( ABSBII.LT.ABS( EI ) ) THEN
307 *
308 * Interchange rows and eliminate.
309 *
310 XR = B( I, I ) / EI
311 XI = B( I+1, I ) / EI
312 B( I, I ) = EI
313 B( I+1, I ) = ZERO
314 DO 150 J = I + 1, N
315 TEMP = B( I+1, J )
316 B( I+1, J ) = B( I, J ) - XR*TEMP
317 B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP
318 B( I, J ) = TEMP
319 B( J+1, I ) = ZERO
320 150 CONTINUE
321 B( I+2, I ) = -WI
322 B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI
323 B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI
324 ELSE
325 *
326 * Eliminate without interchanging rows.
327 *
328 IF( ABSBII.EQ.ZERO ) THEN
329 B( I, I ) = EPS3
330 B( I+1, I ) = ZERO
331 ABSBII = EPS3
332 END IF
333 EI = ( EI / ABSBII ) / ABSBII
334 XR = B( I, I )*EI
335 XI = -B( I+1, I )*EI
336 DO 160 J = I + 1, N
337 B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) +
338 $ XI*B( J+1, I )
339 B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J )
340 160 CONTINUE
341 B( I+2, I+1 ) = B( I+2, I+1 ) - WI
342 END IF
343 *
344 * Compute 1-norm of offdiagonal elements of i-th row.
345 *
346 WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) +
347 $ DASUM( N-I, B( I+2, I ), 1 )
348 170 CONTINUE
349 IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO )
350 $ B( N, N ) = EPS3
351 WORK( N ) = ZERO
352 *
353 I1 = N
354 I2 = 1
355 I3 = -1
356 ELSE
357 *
358 * UL decomposition with partial pivoting of conjg(B),
359 * replacing zero pivots by EPS3.
360 *
361 * The imaginary part of the (i,j)-th element of U is stored in
362 * B(j+1,i).
363 *
364 B( N+1, N ) = WI
365 DO 180 J = 1, N - 1
366 B( N+1, J ) = ZERO
367 180 CONTINUE
368 *
369 DO 210 J = N, 2, -1
370 EJ = H( J, J-1 )
371 ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) )
372 IF( ABSBJJ.LT.ABS( EJ ) ) THEN
373 *
374 * Interchange columns and eliminate
375 *
376 XR = B( J, J ) / EJ
377 XI = B( J+1, J ) / EJ
378 B( J, J ) = EJ
379 B( J+1, J ) = ZERO
380 DO 190 I = 1, J - 1
381 TEMP = B( I, J-1 )
382 B( I, J-1 ) = B( I, J ) - XR*TEMP
383 B( J, I ) = B( J+1, I ) - XI*TEMP
384 B( I, J ) = TEMP
385 B( J+1, I ) = ZERO
386 190 CONTINUE
387 B( J+1, J-1 ) = WI
388 B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI
389 B( J, J-1 ) = B( J, J-1 ) - XR*WI
390 ELSE
391 *
392 * Eliminate without interchange.
393 *
394 IF( ABSBJJ.EQ.ZERO ) THEN
395 B( J, J ) = EPS3
396 B( J+1, J ) = ZERO
397 ABSBJJ = EPS3
398 END IF
399 EJ = ( EJ / ABSBJJ ) / ABSBJJ
400 XR = B( J, J )*EJ
401 XI = -B( J+1, J )*EJ
402 DO 200 I = 1, J - 1
403 B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) +
404 $ XI*B( J+1, I )
405 B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J )
406 200 CONTINUE
407 B( J, J-1 ) = B( J, J-1 ) + WI
408 END IF
409 *
410 * Compute 1-norm of offdiagonal elements of j-th column.
411 *
412 WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) +
413 $ DASUM( J-1, B( J+1, 1 ), LDB )
414 210 CONTINUE
415 IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
416 $ B( 1, 1 ) = EPS3
417 WORK( 1 ) = ZERO
418 *
419 I1 = 1
420 I2 = N
421 I3 = 1
422 END IF
423 *
424 DO 270 ITS = 1, N
425 SCALE = ONE
426 VMAX = ONE
427 VCRIT = BIGNUM
428 *
429 * Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
430 * or U**T*(xr,xi) = scale*(vr,vi) for a left eigenvector,
431 * overwriting (xr,xi) on (vr,vi).
432 *
433 DO 250 I = I1, I2, I3
434 *
435 IF( WORK( I ).GT.VCRIT ) THEN
436 REC = ONE / VMAX
437 CALL DSCAL( N, REC, VR, 1 )
438 CALL DSCAL( N, REC, VI, 1 )
439 SCALE = SCALE*REC
440 VMAX = ONE
441 VCRIT = BIGNUM
442 END IF
443 *
444 XR = VR( I )
445 XI = VI( I )
446 IF( RIGHTV ) THEN
447 DO 220 J = I + 1, N
448 XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J )
449 XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J )
450 220 CONTINUE
451 ELSE
452 DO 230 J = 1, I - 1
453 XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J )
454 XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J )
455 230 CONTINUE
456 END IF
457 *
458 W = ABS( B( I, I ) ) + ABS( B( I+1, I ) )
459 IF( W.GT.SMLNUM ) THEN
460 IF( W.LT.ONE ) THEN
461 W1 = ABS( XR ) + ABS( XI )
462 IF( W1.GT.W*BIGNUM ) THEN
463 REC = ONE / W1
464 CALL DSCAL( N, REC, VR, 1 )
465 CALL DSCAL( N, REC, VI, 1 )
466 XR = VR( I )
467 XI = VI( I )
468 SCALE = SCALE*REC
469 VMAX = VMAX*REC
470 END IF
471 END IF
472 *
473 * Divide by diagonal element of B.
474 *
475 CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ),
476 $ VI( I ) )
477 VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX )
478 VCRIT = BIGNUM / VMAX
479 ELSE
480 DO 240 J = 1, N
481 VR( J ) = ZERO
482 VI( J ) = ZERO
483 240 CONTINUE
484 VR( I ) = ONE
485 VI( I ) = ONE
486 SCALE = ZERO
487 VMAX = ONE
488 VCRIT = BIGNUM
489 END IF
490 250 CONTINUE
491 *
492 * Test for sufficient growth in the norm of (VR,VI).
493 *
494 VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 )
495 IF( VNORM.GE.GROWTO*SCALE )
496 $ GO TO 280
497 *
498 * Choose a new orthogonal starting vector and try again.
499 *
500 Y = EPS3 / ( ROOTN+ONE )
501 VR( 1 ) = EPS3
502 VI( 1 ) = ZERO
503 *
504 DO 260 I = 2, N
505 VR( I ) = Y
506 VI( I ) = ZERO
507 260 CONTINUE
508 VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
509 270 CONTINUE
510 *
511 * Failure to find eigenvector in N iterations
512 *
513 INFO = 1
514 *
515 280 CONTINUE
516 *
517 * Normalize eigenvector.
518 *
519 VNORM = ZERO
520 DO 290 I = 1, N
521 VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) )
522 290 CONTINUE
523 CALL DSCAL( N, ONE / VNORM, VR, 1 )
524 CALL DSCAL( N, ONE / VNORM, VI, 1 )
525 *
526 END IF
527 *
528 RETURN
529 *
530 * End of DLAEIN
531 *
532 END