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