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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhsein.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack