1: *> \brief \b ZHSEIN
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZHSEIN + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhsein.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhsein.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhsein.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
22: * LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
23: * IFAILR, INFO )
24: *
25: * .. Scalar Arguments ..
26: * CHARACTER EIGSRC, INITV, SIDE
27: * INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
28: * ..
29: * .. Array Arguments ..
30: * LOGICAL SELECT( * )
31: * INTEGER IFAILL( * ), IFAILR( * )
32: * DOUBLE PRECISION RWORK( * )
33: * COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
34: * $ W( * ), WORK( * )
35: * ..
36: *
37: *
38: *> \par Purpose:
39: * =============
40: *>
41: *> \verbatim
42: *>
43: *> ZHSEIN uses inverse iteration to find specified right and/or left
44: *> eigenvectors of a complex upper Hessenberg matrix H.
45: *>
46: *> The right eigenvector x and the left eigenvector y of the matrix H
47: *> corresponding to an eigenvalue w are defined by:
48: *>
49: *> H * x = w * x, y**h * H = w * y**h
50: *>
51: *> where y**h denotes the conjugate transpose of the vector y.
52: *> \endverbatim
53: *
54: * Arguments:
55: * ==========
56: *
57: *> \param[in] SIDE
58: *> \verbatim
59: *> SIDE is CHARACTER*1
60: *> = 'R': compute right eigenvectors only;
61: *> = 'L': compute left eigenvectors only;
62: *> = 'B': compute both right and left eigenvectors.
63: *> \endverbatim
64: *>
65: *> \param[in] EIGSRC
66: *> \verbatim
67: *> EIGSRC is CHARACTER*1
68: *> Specifies the source of eigenvalues supplied in W:
69: *> = 'Q': the eigenvalues were found using ZHSEQR; thus, if
70: *> H has zero subdiagonal elements, and so is
71: *> block-triangular, then the j-th eigenvalue can be
72: *> assumed to be an eigenvalue of the block containing
73: *> the j-th row/column. This property allows ZHSEIN to
74: *> perform inverse iteration on just one diagonal block.
75: *> = 'N': no assumptions are made on the correspondence
76: *> between eigenvalues and diagonal blocks. In this
77: *> case, ZHSEIN must always perform inverse iteration
78: *> using the whole matrix H.
79: *> \endverbatim
80: *>
81: *> \param[in] INITV
82: *> \verbatim
83: *> INITV is CHARACTER*1
84: *> = 'N': no initial vectors are supplied;
85: *> = 'U': user-supplied initial vectors are stored in the arrays
86: *> VL and/or VR.
87: *> \endverbatim
88: *>
89: *> \param[in] SELECT
90: *> \verbatim
91: *> SELECT is LOGICAL array, dimension (N)
92: *> Specifies the eigenvectors to be computed. To select the
93: *> eigenvector corresponding to the eigenvalue W(j),
94: *> SELECT(j) must be set to .TRUE..
95: *> \endverbatim
96: *>
97: *> \param[in] N
98: *> \verbatim
99: *> N is INTEGER
100: *> The order of the matrix H. N >= 0.
101: *> \endverbatim
102: *>
103: *> \param[in] H
104: *> \verbatim
105: *> H is COMPLEX*16 array, dimension (LDH,N)
106: *> The upper Hessenberg matrix H.
107: *> If a NaN is detected in H, the routine will return with INFO=-6.
108: *> \endverbatim
109: *>
110: *> \param[in] LDH
111: *> \verbatim
112: *> LDH is INTEGER
113: *> The leading dimension of the array H. LDH >= max(1,N).
114: *> \endverbatim
115: *>
116: *> \param[in,out] W
117: *> \verbatim
118: *> W is COMPLEX*16 array, dimension (N)
119: *> On entry, the eigenvalues of H.
120: *> On exit, the real parts of W may have been altered since
121: *> close eigenvalues are perturbed slightly in searching for
122: *> independent eigenvectors.
123: *> \endverbatim
124: *>
125: *> \param[in,out] VL
126: *> \verbatim
127: *> VL is COMPLEX*16 array, dimension (LDVL,MM)
128: *> On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
129: *> contain starting vectors for the inverse iteration for the
130: *> left eigenvectors; the starting vector for each eigenvector
131: *> must be in the same column in which the eigenvector will be
132: *> stored.
133: *> On exit, if SIDE = 'L' or 'B', the left eigenvectors
134: *> specified by SELECT will be stored consecutively in the
135: *> columns of VL, in the same order as their eigenvalues.
136: *> If SIDE = 'R', VL is not referenced.
137: *> \endverbatim
138: *>
139: *> \param[in] LDVL
140: *> \verbatim
141: *> LDVL is INTEGER
142: *> The leading dimension of the array VL.
143: *> LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
144: *> \endverbatim
145: *>
146: *> \param[in,out] VR
147: *> \verbatim
148: *> VR is COMPLEX*16 array, dimension (LDVR,MM)
149: *> On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
150: *> contain starting vectors for the inverse iteration for the
151: *> right eigenvectors; the starting vector for each eigenvector
152: *> must be in the same column in which the eigenvector will be
153: *> stored.
154: *> On exit, if SIDE = 'R' or 'B', the right eigenvectors
155: *> specified by SELECT will be stored consecutively in the
156: *> columns of VR, in the same order as their eigenvalues.
157: *> If SIDE = 'L', VR is not referenced.
158: *> \endverbatim
159: *>
160: *> \param[in] LDVR
161: *> \verbatim
162: *> LDVR is INTEGER
163: *> The leading dimension of the array VR.
164: *> LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
165: *> \endverbatim
166: *>
167: *> \param[in] MM
168: *> \verbatim
169: *> MM is INTEGER
170: *> The number of columns in the arrays VL and/or VR. MM >= M.
171: *> \endverbatim
172: *>
173: *> \param[out] M
174: *> \verbatim
175: *> M is INTEGER
176: *> The number of columns in the arrays VL and/or VR required to
177: *> store the eigenvectors (= the number of .TRUE. elements in
178: *> SELECT).
179: *> \endverbatim
180: *>
181: *> \param[out] WORK
182: *> \verbatim
183: *> WORK is COMPLEX*16 array, dimension (N*N)
184: *> \endverbatim
185: *>
186: *> \param[out] RWORK
187: *> \verbatim
188: *> RWORK is DOUBLE PRECISION array, dimension (N)
189: *> \endverbatim
190: *>
191: *> \param[out] IFAILL
192: *> \verbatim
193: *> IFAILL is INTEGER array, dimension (MM)
194: *> If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
195: *> eigenvector in the i-th column of VL (corresponding to the
196: *> eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
197: *> eigenvector converged satisfactorily.
198: *> If SIDE = 'R', IFAILL is not referenced.
199: *> \endverbatim
200: *>
201: *> \param[out] IFAILR
202: *> \verbatim
203: *> IFAILR is INTEGER array, dimension (MM)
204: *> If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
205: *> eigenvector in the i-th column of VR (corresponding to the
206: *> eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
207: *> eigenvector converged satisfactorily.
208: *> If SIDE = 'L', IFAILR is not referenced.
209: *> \endverbatim
210: *>
211: *> \param[out] INFO
212: *> \verbatim
213: *> INFO is INTEGER
214: *> = 0: successful exit
215: *> < 0: if INFO = -i, the i-th argument had an illegal value
216: *> > 0: if INFO = i, i is the number of eigenvectors which
217: *> failed to converge; see IFAILL and IFAILR for further
218: *> details.
219: *> \endverbatim
220: *
221: * Authors:
222: * ========
223: *
224: *> \author Univ. of Tennessee
225: *> \author Univ. of California Berkeley
226: *> \author Univ. of Colorado Denver
227: *> \author NAG Ltd.
228: *
229: *> \date November 2013
230: *
231: *> \ingroup complex16OTHERcomputational
232: *
233: *> \par Further Details:
234: * =====================
235: *>
236: *> \verbatim
237: *>
238: *> Each eigenvector is normalized so that the element of largest
239: *> magnitude has magnitude 1; here the magnitude of a complex number
240: *> (x,y) is taken to be |x|+|y|.
241: *> \endverbatim
242: *>
243: * =====================================================================
244: SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
245: $ LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
246: $ IFAILR, INFO )
247: *
248: * -- LAPACK computational routine (version 3.5.0) --
249: * -- LAPACK is a software package provided by Univ. of Tennessee, --
250: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
251: * November 2013
252: *
253: * .. Scalar Arguments ..
254: CHARACTER EIGSRC, INITV, SIDE
255: INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
256: * ..
257: * .. Array Arguments ..
258: LOGICAL SELECT( * )
259: INTEGER IFAILL( * ), IFAILR( * )
260: DOUBLE PRECISION RWORK( * )
261: COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
262: $ W( * ), WORK( * )
263: * ..
264: *
265: * =====================================================================
266: *
267: * .. Parameters ..
268: COMPLEX*16 ZERO
269: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
270: DOUBLE PRECISION RZERO
271: PARAMETER ( RZERO = 0.0D+0 )
272: * ..
273: * .. Local Scalars ..
274: LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
275: INTEGER I, IINFO, K, KL, KLN, KR, KS, LDWORK
276: DOUBLE PRECISION EPS3, HNORM, SMLNUM, ULP, UNFL
277: COMPLEX*16 CDUM, WK
278: * ..
279: * .. External Functions ..
280: LOGICAL LSAME, DISNAN
281: DOUBLE PRECISION DLAMCH, ZLANHS
282: EXTERNAL LSAME, DLAMCH, ZLANHS, DISNAN
283: * ..
284: * .. External Subroutines ..
285: EXTERNAL XERBLA, ZLAEIN
286: * ..
287: * .. Intrinsic Functions ..
288: INTRINSIC ABS, DBLE, DIMAG, MAX
289: * ..
290: * .. Statement Functions ..
291: DOUBLE PRECISION CABS1
292: * ..
293: * .. Statement Function definitions ..
294: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
295: * ..
296: * .. Executable Statements ..
297: *
298: * Decode and test the input parameters.
299: *
300: BOTHV = LSAME( SIDE, 'B' )
301: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
302: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
303: *
304: FROMQR = LSAME( EIGSRC, 'Q' )
305: *
306: NOINIT = LSAME( INITV, 'N' )
307: *
308: * Set M to the number of columns required to store the selected
309: * eigenvectors.
310: *
311: M = 0
312: DO 10 K = 1, N
313: IF( SELECT( K ) )
314: $ M = M + 1
315: 10 CONTINUE
316: *
317: INFO = 0
318: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
319: INFO = -1
320: ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
321: INFO = -2
322: ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
323: INFO = -3
324: ELSE IF( N.LT.0 ) THEN
325: INFO = -5
326: ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
327: INFO = -7
328: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
329: INFO = -10
330: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
331: INFO = -12
332: ELSE IF( MM.LT.M ) THEN
333: INFO = -13
334: END IF
335: IF( INFO.NE.0 ) THEN
336: CALL XERBLA( 'ZHSEIN', -INFO )
337: RETURN
338: END IF
339: *
340: * Quick return if possible.
341: *
342: IF( N.EQ.0 )
343: $ RETURN
344: *
345: * Set machine-dependent constants.
346: *
347: UNFL = DLAMCH( 'Safe minimum' )
348: ULP = DLAMCH( 'Precision' )
349: SMLNUM = UNFL*( N / ULP )
350: *
351: LDWORK = N
352: *
353: KL = 1
354: KLN = 0
355: IF( FROMQR ) THEN
356: KR = 0
357: ELSE
358: KR = N
359: END IF
360: KS = 1
361: *
362: DO 100 K = 1, N
363: IF( SELECT( K ) ) THEN
364: *
365: * Compute eigenvector(s) corresponding to W(K).
366: *
367: IF( FROMQR ) THEN
368: *
369: * If affiliation of eigenvalues is known, check whether
370: * the matrix splits.
371: *
372: * Determine KL and KR such that 1 <= KL <= K <= KR <= N
373: * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
374: * KR = N).
375: *
376: * Then inverse iteration can be performed with the
377: * submatrix H(KL:N,KL:N) for a left eigenvector, and with
378: * the submatrix H(1:KR,1:KR) for a right eigenvector.
379: *
380: DO 20 I = K, KL + 1, -1
381: IF( H( I, I-1 ).EQ.ZERO )
382: $ GO TO 30
383: 20 CONTINUE
384: 30 CONTINUE
385: KL = I
386: IF( K.GT.KR ) THEN
387: DO 40 I = K, N - 1
388: IF( H( I+1, I ).EQ.ZERO )
389: $ GO TO 50
390: 40 CONTINUE
391: 50 CONTINUE
392: KR = I
393: END IF
394: END IF
395: *
396: IF( KL.NE.KLN ) THEN
397: KLN = KL
398: *
399: * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
400: * has not ben computed before.
401: *
402: HNORM = ZLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
403: IF( DISNAN( HNORM ) ) THEN
404: INFO = -6
405: RETURN
406: ELSE IF( HNORM.GT.RZERO ) THEN
407: EPS3 = HNORM*ULP
408: ELSE
409: EPS3 = SMLNUM
410: END IF
411: END IF
412: *
413: * Perturb eigenvalue if it is close to any previous
414: * selected eigenvalues affiliated to the submatrix
415: * H(KL:KR,KL:KR). Close roots are modified by EPS3.
416: *
417: WK = W( K )
418: 60 CONTINUE
419: DO 70 I = K - 1, KL, -1
420: IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN
421: WK = WK + EPS3
422: GO TO 60
423: END IF
424: 70 CONTINUE
425: W( K ) = WK
426: *
427: IF( LEFTV ) THEN
428: *
429: * Compute left eigenvector.
430: *
431: CALL ZLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
432: $ WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3,
433: $ SMLNUM, IINFO )
434: IF( IINFO.GT.0 ) THEN
435: INFO = INFO + 1
436: IFAILL( KS ) = K
437: ELSE
438: IFAILL( KS ) = 0
439: END IF
440: DO 80 I = 1, KL - 1
441: VL( I, KS ) = ZERO
442: 80 CONTINUE
443: END IF
444: IF( RIGHTV ) THEN
445: *
446: * Compute right eigenvector.
447: *
448: CALL ZLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ),
449: $ WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO )
450: IF( IINFO.GT.0 ) THEN
451: INFO = INFO + 1
452: IFAILR( KS ) = K
453: ELSE
454: IFAILR( KS ) = 0
455: END IF
456: DO 90 I = KR + 1, N
457: VR( I, KS ) = ZERO
458: 90 CONTINUE
459: END IF
460: KS = KS + 1
461: END IF
462: 100 CONTINUE
463: *
464: RETURN
465: *
466: * End of ZHSEIN
467: *
468: 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