LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zhsein.f
Go to the documentation of this file.
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
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*> \ingroup hsein
230*
231*> \par Further Details:
232* =====================
233*>
234*> \verbatim
235*>
236*> Each eigenvector is normalized so that the element of largest
237*> magnitude has magnitude 1; here the magnitude of a complex number
238*> (x,y) is taken to be |x|+|y|.
239*> \endverbatim
240*>
241* =====================================================================
242 SUBROUTINE zhsein( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
243 \$ LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
244 \$ IFAILR, INFO )
245*
246* -- LAPACK computational routine --
247* -- LAPACK is a software package provided by Univ. of Tennessee, --
248* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
249*
250* .. Scalar Arguments ..
251 CHARACTER EIGSRC, INITV, SIDE
252 INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
253* ..
254* .. Array Arguments ..
255 LOGICAL SELECT( * )
256 INTEGER IFAILL( * ), IFAILR( * )
257 DOUBLE PRECISION RWORK( * )
258 COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
259 \$ w( * ), work( * )
260* ..
261*
262* =====================================================================
263*
264* .. Parameters ..
265 COMPLEX*16 ZERO
266 PARAMETER ( ZERO = ( 0.0d+0, 0.0d+0 ) )
267 DOUBLE PRECISION RZERO
268 parameter( rzero = 0.0d+0 )
269* ..
270* .. Local Scalars ..
271 LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
272 INTEGER I, IINFO, K, KL, KLN, KR, KS, LDWORK
273 DOUBLE PRECISION EPS3, HNORM, SMLNUM, ULP, UNFL
274 COMPLEX*16 CDUM, WK
275* ..
276* .. External Functions ..
277 LOGICAL LSAME, DISNAN
278 DOUBLE PRECISION DLAMCH, ZLANHS
279 EXTERNAL lsame, dlamch, zlanhs, disnan
280* ..
281* .. External Subroutines ..
282 EXTERNAL xerbla, zlaein
283* ..
284* .. Intrinsic Functions ..
285 INTRINSIC abs, dble, dimag, max
286* ..
287* .. Statement Functions ..
288 DOUBLE PRECISION CABS1
289* ..
290* .. Statement Function definitions ..
291 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
292* ..
293* .. Executable Statements ..
294*
295* Decode and test the input parameters.
296*
297 bothv = lsame( side, 'B' )
298 rightv = lsame( side, 'R' ) .OR. bothv
299 leftv = lsame( side, 'L' ) .OR. bothv
300*
301 fromqr = lsame( eigsrc, 'Q' )
302*
303 noinit = lsame( initv, 'N' )
304*
305* Set M to the number of columns required to store the selected
306* eigenvectors.
307*
308 m = 0
309 DO 10 k = 1, n
310 IF( SELECT( k ) )
311 \$ m = m + 1
312 10 CONTINUE
313*
314 info = 0
315 IF( .NOT.rightv .AND. .NOT.leftv ) THEN
316 info = -1
317 ELSE IF( .NOT.fromqr .AND. .NOT.lsame( eigsrc, 'N' ) ) THEN
318 info = -2
319 ELSE IF( .NOT.noinit .AND. .NOT.lsame( initv, 'U' ) ) THEN
320 info = -3
321 ELSE IF( n.LT.0 ) THEN
322 info = -5
323 ELSE IF( ldh.LT.max( 1, n ) ) THEN
324 info = -7
325 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
326 info = -10
327 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
328 info = -12
329 ELSE IF( mm.LT.m ) THEN
330 info = -13
331 END IF
332 IF( info.NE.0 ) THEN
333 CALL xerbla( 'ZHSEIN', -info )
334 RETURN
335 END IF
336*
337* Quick return if possible.
338*
339 IF( n.EQ.0 )
340 \$ RETURN
341*
342* Set machine-dependent constants.
343*
344 unfl = dlamch( 'Safe minimum' )
345 ulp = dlamch( 'Precision' )
346 smlnum = unfl*( n / ulp )
347*
348 ldwork = n
349*
350 kl = 1
351 kln = 0
352 IF( fromqr ) THEN
353 kr = 0
354 ELSE
355 kr = n
356 END IF
357 ks = 1
358*
359 DO 100 k = 1, n
360 IF( SELECT( k ) ) THEN
361*
362* Compute eigenvector(s) corresponding to W(K).
363*
364 IF( fromqr ) THEN
365*
366* If affiliation of eigenvalues is known, check whether
367* the matrix splits.
368*
369* Determine KL and KR such that 1 <= KL <= K <= KR <= N
370* and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
371* KR = N).
372*
373* Then inverse iteration can be performed with the
374* submatrix H(KL:N,KL:N) for a left eigenvector, and with
375* the submatrix H(1:KR,1:KR) for a right eigenvector.
376*
377 DO 20 i = k, kl + 1, -1
378 IF( h( i, i-1 ).EQ.zero )
379 \$ GO TO 30
380 20 CONTINUE
381 30 CONTINUE
382 kl = i
383 IF( k.GT.kr ) THEN
384 DO 40 i = k, n - 1
385 IF( h( i+1, i ).EQ.zero )
386 \$ GO TO 50
387 40 CONTINUE
388 50 CONTINUE
389 kr = i
390 END IF
391 END IF
392*
393 IF( kl.NE.kln ) THEN
394 kln = kl
395*
396* Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
397* has not ben computed before.
398*
399 hnorm = zlanhs( 'I', kr-kl+1, h( kl, kl ), ldh, rwork )
400 IF( disnan( hnorm ) ) THEN
401 info = -6
402 RETURN
403 ELSE IF( hnorm.GT.rzero ) THEN
404 eps3 = hnorm*ulp
405 ELSE
406 eps3 = smlnum
407 END IF
408 END IF
409*
410* Perturb eigenvalue if it is close to any previous
411* selected eigenvalues affiliated to the submatrix
412* H(KL:KR,KL:KR). Close roots are modified by EPS3.
413*
414 wk = w( k )
415 60 CONTINUE
416 DO 70 i = k - 1, kl, -1
417 IF( SELECT( i ) .AND. cabs1( w( i )-wk ).LT.eps3 ) THEN
418 wk = wk + eps3
419 GO TO 60
420 END IF
421 70 CONTINUE
422 w( k ) = wk
423*
424 IF( leftv ) THEN
425*
426* Compute left eigenvector.
427*
428 CALL zlaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,
429 \$ wk, vl( kl, ks ), work, ldwork, rwork, eps3,
430 \$ smlnum, iinfo )
431 IF( iinfo.GT.0 ) THEN
432 info = info + 1
433 ifaill( ks ) = k
434 ELSE
435 ifaill( ks ) = 0
436 END IF
437 DO 80 i = 1, kl - 1
438 vl( i, ks ) = zero
439 80 CONTINUE
440 END IF
441 IF( rightv ) THEN
442*
443* Compute right eigenvector.
444*
445 CALL zlaein( .true., noinit, kr, h, ldh, wk, vr( 1, ks ),
446 \$ work, ldwork, rwork, eps3, smlnum, iinfo )
447 IF( iinfo.GT.0 ) THEN
448 info = info + 1
449 ifailr( ks ) = k
450 ELSE
451 ifailr( ks ) = 0
452 END IF
453 DO 90 i = kr + 1, n
454 vr( i, ks ) = zero
455 90 CONTINUE
456 END IF
457 ks = ks + 1
458 END IF
459 100 CONTINUE
460*
461 RETURN
462*
463* End of ZHSEIN
464*
465 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zhsein(side, eigsrc, initv, select, n, h, ldh, w, vl, ldvl, vr, ldvr, mm, m, work, rwork, ifaill, ifailr, info)
ZHSEIN
Definition zhsein.f:245
subroutine zlaein(rightv, noinit, n, h, ldh, w, v, b, ldb, rwork, eps3, smlnum, info)
ZLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iterat...
Definition zlaein.f:149