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.2.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * June 2010
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'*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'*(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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaein.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack