1: *> \brief \b ZLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZLAEIN + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaein.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaein.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaein.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
22: * EPS3, SMLNUM, INFO )
23: *
24: * .. Scalar Arguments ..
25: * LOGICAL NOINIT, RIGHTV
26: * INTEGER INFO, LDB, LDH, N
27: * DOUBLE PRECISION EPS3, SMLNUM
28: * COMPLEX*16 W
29: * ..
30: * .. Array Arguments ..
31: * DOUBLE PRECISION RWORK( * )
32: * COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * )
33: * ..
34: *
35: *
36: *> \par Purpose:
37: * =============
38: *>
39: *> \verbatim
40: *>
41: *> ZLAEIN uses inverse iteration to find a right or left eigenvector
42: *> corresponding to the eigenvalue W of a complex upper Hessenberg
43: *> matrix H.
44: *> \endverbatim
45: *
46: * Arguments:
47: * ==========
48: *
49: *> \param[in] RIGHTV
50: *> \verbatim
51: *> RIGHTV is LOGICAL
52: *> = .TRUE. : compute right eigenvector;
53: *> = .FALSE.: compute left eigenvector.
54: *> \endverbatim
55: *>
56: *> \param[in] NOINIT
57: *> \verbatim
58: *> NOINIT is LOGICAL
59: *> = .TRUE. : no initial vector supplied in V
60: *> = .FALSE.: initial vector supplied in V.
61: *> \endverbatim
62: *>
63: *> \param[in] N
64: *> \verbatim
65: *> N is INTEGER
66: *> The order of the matrix H. N >= 0.
67: *> \endverbatim
68: *>
69: *> \param[in] H
70: *> \verbatim
71: *> H is COMPLEX*16 array, dimension (LDH,N)
72: *> The upper Hessenberg matrix H.
73: *> \endverbatim
74: *>
75: *> \param[in] LDH
76: *> \verbatim
77: *> LDH is INTEGER
78: *> The leading dimension of the array H. LDH >= max(1,N).
79: *> \endverbatim
80: *>
81: *> \param[in] W
82: *> \verbatim
83: *> W is COMPLEX*16
84: *> The eigenvalue of H whose corresponding right or left
85: *> eigenvector is to be computed.
86: *> \endverbatim
87: *>
88: *> \param[in,out] V
89: *> \verbatim
90: *> V is COMPLEX*16 array, dimension (N)
91: *> On entry, if NOINIT = .FALSE., V must contain a starting
92: *> vector for inverse iteration; otherwise V need not be set.
93: *> On exit, V contains the computed eigenvector, normalized so
94: *> that the component of largest magnitude has magnitude 1; here
95: *> the magnitude of a complex number (x,y) is taken to be
96: *> |x| + |y|.
97: *> \endverbatim
98: *>
99: *> \param[out] B
100: *> \verbatim
101: *> B is COMPLEX*16 array, dimension (LDB,N)
102: *> \endverbatim
103: *>
104: *> \param[in] LDB
105: *> \verbatim
106: *> LDB is INTEGER
107: *> The leading dimension of the array B. LDB >= max(1,N).
108: *> \endverbatim
109: *>
110: *> \param[out] RWORK
111: *> \verbatim
112: *> RWORK is DOUBLE PRECISION array, dimension (N)
113: *> \endverbatim
114: *>
115: *> \param[in] EPS3
116: *> \verbatim
117: *> EPS3 is DOUBLE PRECISION
118: *> A small machine-dependent value which is used to perturb
119: *> close eigenvalues, and to replace zero pivots.
120: *> \endverbatim
121: *>
122: *> \param[in] SMLNUM
123: *> \verbatim
124: *> SMLNUM is DOUBLE PRECISION
125: *> A machine-dependent value close to the underflow threshold.
126: *> \endverbatim
127: *>
128: *> \param[out] INFO
129: *> \verbatim
130: *> INFO is INTEGER
131: *> = 0: successful exit
132: *> = 1: inverse iteration did not converge; V is set to the
133: *> last iterate.
134: *> \endverbatim
135: *
136: * Authors:
137: * ========
138: *
139: *> \author Univ. of Tennessee
140: *> \author Univ. of California Berkeley
141: *> \author Univ. of Colorado Denver
142: *> \author NAG Ltd.
143: *
144: *> \ingroup complex16OTHERauxiliary
145: *
146: * =====================================================================
147: SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
148: $ EPS3, SMLNUM, INFO )
149: *
150: * -- LAPACK auxiliary routine --
151: * -- LAPACK is a software package provided by Univ. of Tennessee, --
152: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
153: *
154: * .. Scalar Arguments ..
155: LOGICAL NOINIT, RIGHTV
156: INTEGER INFO, LDB, LDH, N
157: DOUBLE PRECISION EPS3, SMLNUM
158: COMPLEX*16 W
159: * ..
160: * .. Array Arguments ..
161: DOUBLE PRECISION RWORK( * )
162: COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * )
163: * ..
164: *
165: * =====================================================================
166: *
167: * .. Parameters ..
168: DOUBLE PRECISION ONE, TENTH
169: PARAMETER ( ONE = 1.0D+0, TENTH = 1.0D-1 )
170: COMPLEX*16 ZERO
171: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
172: * ..
173: * .. Local Scalars ..
174: CHARACTER NORMIN, TRANS
175: INTEGER I, IERR, ITS, J
176: DOUBLE PRECISION GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM
177: COMPLEX*16 CDUM, EI, EJ, TEMP, X
178: * ..
179: * .. External Functions ..
180: INTEGER IZAMAX
181: DOUBLE PRECISION DZASUM, DZNRM2
182: COMPLEX*16 ZLADIV
183: EXTERNAL IZAMAX, DZASUM, DZNRM2, ZLADIV
184: * ..
185: * .. External Subroutines ..
186: EXTERNAL ZDSCAL, ZLATRS
187: * ..
188: * .. Intrinsic Functions ..
189: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
190: * ..
191: * .. Statement Functions ..
192: DOUBLE PRECISION CABS1
193: * ..
194: * .. Statement Function definitions ..
195: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
196: * ..
197: * .. Executable Statements ..
198: *
199: INFO = 0
200: *
201: * GROWTO is the threshold used in the acceptance test for an
202: * eigenvector.
203: *
204: ROOTN = SQRT( DBLE( N ) )
205: GROWTO = TENTH / ROOTN
206: NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
207: *
208: * Form B = H - W*I (except that the subdiagonal elements are not
209: * stored).
210: *
211: DO 20 J = 1, N
212: DO 10 I = 1, J - 1
213: B( I, J ) = H( I, J )
214: 10 CONTINUE
215: B( J, J ) = H( J, J ) - W
216: 20 CONTINUE
217: *
218: IF( NOINIT ) THEN
219: *
220: * Initialize V.
221: *
222: DO 30 I = 1, N
223: V( I ) = EPS3
224: 30 CONTINUE
225: ELSE
226: *
227: * Scale supplied initial vector.
228: *
229: VNORM = DZNRM2( N, V, 1 )
230: CALL ZDSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 )
231: END IF
232: *
233: IF( RIGHTV ) THEN
234: *
235: * LU decomposition with partial pivoting of B, replacing zero
236: * pivots by EPS3.
237: *
238: DO 60 I = 1, N - 1
239: EI = H( I+1, I )
240: IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN
241: *
242: * Interchange rows and eliminate.
243: *
244: X = ZLADIV( B( I, I ), EI )
245: B( I, I ) = EI
246: DO 40 J = I + 1, N
247: TEMP = B( I+1, J )
248: B( I+1, J ) = B( I, J ) - X*TEMP
249: B( I, J ) = TEMP
250: 40 CONTINUE
251: ELSE
252: *
253: * Eliminate without interchange.
254: *
255: IF( B( I, I ).EQ.ZERO )
256: $ B( I, I ) = EPS3
257: X = ZLADIV( EI, B( I, I ) )
258: IF( X.NE.ZERO ) THEN
259: DO 50 J = I + 1, N
260: B( I+1, J ) = B( I+1, J ) - X*B( I, J )
261: 50 CONTINUE
262: END IF
263: END IF
264: 60 CONTINUE
265: IF( B( N, N ).EQ.ZERO )
266: $ B( N, N ) = EPS3
267: *
268: TRANS = 'N'
269: *
270: ELSE
271: *
272: * UL decomposition with partial pivoting of B, replacing zero
273: * pivots by EPS3.
274: *
275: DO 90 J = N, 2, -1
276: EJ = H( J, J-1 )
277: IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN
278: *
279: * Interchange columns and eliminate.
280: *
281: X = ZLADIV( B( J, J ), EJ )
282: B( J, J ) = EJ
283: DO 70 I = 1, J - 1
284: TEMP = B( I, J-1 )
285: B( I, J-1 ) = B( I, J ) - X*TEMP
286: B( I, J ) = TEMP
287: 70 CONTINUE
288: ELSE
289: *
290: * Eliminate without interchange.
291: *
292: IF( B( J, J ).EQ.ZERO )
293: $ B( J, J ) = EPS3
294: X = ZLADIV( EJ, B( J, J ) )
295: IF( X.NE.ZERO ) THEN
296: DO 80 I = 1, J - 1
297: B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
298: 80 CONTINUE
299: END IF
300: END IF
301: 90 CONTINUE
302: IF( B( 1, 1 ).EQ.ZERO )
303: $ B( 1, 1 ) = EPS3
304: *
305: TRANS = 'C'
306: *
307: END IF
308: *
309: NORMIN = 'N'
310: DO 110 ITS = 1, N
311: *
312: * Solve U*x = scale*v for a right eigenvector
313: * or U**H *x = scale*v for a left eigenvector,
314: * overwriting x on v.
315: *
316: CALL ZLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V,
317: $ SCALE, RWORK, IERR )
318: NORMIN = 'Y'
319: *
320: * Test for sufficient growth in the norm of v.
321: *
322: VNORM = DZASUM( N, V, 1 )
323: IF( VNORM.GE.GROWTO*SCALE )
324: $ GO TO 120
325: *
326: * Choose new orthogonal starting vector and try again.
327: *
328: RTEMP = EPS3 / ( ROOTN+ONE )
329: V( 1 ) = EPS3
330: DO 100 I = 2, N
331: V( I ) = RTEMP
332: 100 CONTINUE
333: V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN
334: 110 CONTINUE
335: *
336: * Failure to find eigenvector in N iterations.
337: *
338: INFO = 1
339: *
340: 120 CONTINUE
341: *
342: * Normalize eigenvector.
343: *
344: I = IZAMAX( N, V, 1 )
345: CALL ZDSCAL( N, ONE / CABS1( V( I ) ), V, 1 )
346: *
347: RETURN
348: *
349: * End of ZLAEIN
350: *
351: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaein.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack