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