LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
zhsein.f
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 *> \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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
zlaein
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
zhsein
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