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