LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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
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 *> \endverbatim
108 *>
109 *> \param[in] LDH
110 *> \verbatim
111 *> LDH is INTEGER
112 *> The leading dimension of the array H. LDH >= max(1,N).
113 *> \endverbatim
114 *>
115 *> \param[in,out] W
116 *> \verbatim
117 *> W is COMPLEX*16 array, dimension (N)
118 *> On entry, the eigenvalues of H.
119 *> On exit, the real parts of W may have been altered since
120 *> close eigenvalues are perturbed slightly in searching for
121 *> independent eigenvectors.
122 *> \endverbatim
123 *>
124 *> \param[in,out] VL
125 *> \verbatim
126 *> VL is COMPLEX*16 array, dimension (LDVL,MM)
127 *> On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
128 *> contain starting vectors for the inverse iteration for the
129 *> left eigenvectors; the starting vector for each eigenvector
130 *> must be in the same column in which the eigenvector will be
131 *> stored.
132 *> On exit, if SIDE = 'L' or 'B', the left eigenvectors
133 *> specified by SELECT will be stored consecutively in the
134 *> columns of VL, in the same order as their eigenvalues.
135 *> If SIDE = 'R', VL is not referenced.
136 *> \endverbatim
137 *>
138 *> \param[in] LDVL
139 *> \verbatim
140 *> LDVL is INTEGER
141 *> The leading dimension of the array VL.
142 *> LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
143 *> \endverbatim
144 *>
145 *> \param[in,out] VR
146 *> \verbatim
147 *> VR is COMPLEX*16 array, dimension (LDVR,MM)
148 *> On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
149 *> contain starting vectors for the inverse iteration for the
150 *> right eigenvectors; the starting vector for each eigenvector
151 *> must be in the same column in which the eigenvector will be
152 *> stored.
153 *> On exit, if SIDE = 'R' or 'B', the right eigenvectors
154 *> specified by SELECT will be stored consecutively in the
155 *> columns of VR, in the same order as their eigenvalues.
156 *> If SIDE = 'L', VR is not referenced.
157 *> \endverbatim
158 *>
159 *> \param[in] LDVR
160 *> \verbatim
161 *> LDVR is INTEGER
162 *> The leading dimension of the array VR.
163 *> LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
164 *> \endverbatim
165 *>
166 *> \param[in] MM
167 *> \verbatim
168 *> MM is INTEGER
169 *> The number of columns in the arrays VL and/or VR. MM >= M.
170 *> \endverbatim
171 *>
172 *> \param[out] M
173 *> \verbatim
174 *> M is INTEGER
175 *> The number of columns in the arrays VL and/or VR required to
176 *> store the eigenvectors (= the number of .TRUE. elements in
177 *> SELECT).
178 *> \endverbatim
179 *>
180 *> \param[out] WORK
181 *> \verbatim
182 *> WORK is COMPLEX*16 array, dimension (N*N)
183 *> \endverbatim
184 *>
185 *> \param[out] RWORK
186 *> \verbatim
187 *> RWORK is DOUBLE PRECISION array, dimension (N)
188 *> \endverbatim
189 *>
190 *> \param[out] IFAILL
191 *> \verbatim
192 *> IFAILL is INTEGER array, dimension (MM)
193 *> If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
194 *> eigenvector in the i-th column of VL (corresponding to the
195 *> eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
196 *> eigenvector converged satisfactorily.
197 *> If SIDE = 'R', IFAILL is not referenced.
198 *> \endverbatim
199 *>
200 *> \param[out] IFAILR
201 *> \verbatim
202 *> IFAILR is INTEGER array, dimension (MM)
203 *> If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
204 *> eigenvector in the i-th column of VR (corresponding to the
205 *> eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
206 *> eigenvector converged satisfactorily.
207 *> If SIDE = 'L', IFAILR is not referenced.
208 *> \endverbatim
209 *>
210 *> \param[out] INFO
211 *> \verbatim
212 *> INFO is INTEGER
213 *> = 0: successful exit
214 *> < 0: if INFO = -i, the i-th argument had an illegal value
215 *> > 0: if INFO = i, i is the number of eigenvectors which
216 *> failed to converge; see IFAILL and IFAILR for further
217 *> details.
218 *> \endverbatim
219 *
220 * Authors:
221 * ========
222 *
223 *> \author Univ. of Tennessee
224 *> \author Univ. of California Berkeley
225 *> \author Univ. of Colorado Denver
226 *> \author NAG Ltd.
227 *
228 *> \date November 2011
229 *
230 *> \ingroup complex16OTHERcomputational
231 *
232 *> \par Further Details:
233 * =====================
234 *>
235 *> \verbatim
236 *>
237 *> Each eigenvector is normalized so that the element of largest
238 *> magnitude has magnitude 1; here the magnitude of a complex number
239 *> (x,y) is taken to be |x|+|y|.
240 *> \endverbatim
241 *>
242 * =====================================================================
243  SUBROUTINE zhsein( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
244  $ ldvl, vr, ldvr, mm, m, work, rwork, ifaill,
245  $ ifailr, info )
246 *
247 * -- LAPACK computational routine (version 3.4.0) --
248 * -- LAPACK is a software package provided by Univ. of Tennessee, --
249 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
250 * November 2011
251 *
252 * .. Scalar Arguments ..
253  CHARACTER eigsrc, initv, side
254  INTEGER info, ldh, ldvl, ldvr, m, mm, n
255 * ..
256 * .. Array Arguments ..
257  LOGICAL select( * )
258  INTEGER ifaill( * ), ifailr( * )
259  DOUBLE PRECISION rwork( * )
260  COMPLEX*16 h( ldh, * ), vl( ldvl, * ), vr( ldvr, * ),
261  $ w( * ), work( * )
262 * ..
263 *
264 * =====================================================================
265 *
266 * .. Parameters ..
267  COMPLEX*16 zero
268  parameter( zero = ( 0.0d+0, 0.0d+0 ) )
269  DOUBLE PRECISION rzero
270  parameter( rzero = 0.0d+0 )
271 * ..
272 * .. Local Scalars ..
273  LOGICAL bothv, fromqr, leftv, noinit, rightv
274  INTEGER i, iinfo, k, kl, kln, kr, ks, ldwork
275  DOUBLE PRECISION eps3, hnorm, smlnum, ulp, unfl
276  COMPLEX*16 cdum, wk
277 * ..
278 * .. External Functions ..
279  LOGICAL lsame
280  DOUBLE PRECISION dlamch, zlanhs
281  EXTERNAL lsame, dlamch, zlanhs
282 * ..
283 * .. External Subroutines ..
284  EXTERNAL xerbla, zlaein
285 * ..
286 * .. Intrinsic Functions ..
287  INTRINSIC abs, dble, dimag, max
288 * ..
289 * .. Statement Functions ..
290  DOUBLE PRECISION cabs1
291 * ..
292 * .. Statement Function definitions ..
293  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
294 * ..
295 * .. Executable Statements ..
296 *
297 * Decode and test the input parameters.
298 *
299  bothv = lsame( side, 'B' )
300  rightv = lsame( side, 'R' ) .OR. bothv
301  leftv = lsame( side, 'L' ) .OR. bothv
302 *
303  fromqr = lsame( eigsrc, 'Q' )
304 *
305  noinit = lsame( initv, 'N' )
306 *
307 * Set M to the number of columns required to store the selected
308 * eigenvectors.
309 *
310  m = 0
311  DO 10 k = 1, n
312  IF( SELECT( k ) )
313  $ m = m + 1
314  10 continue
315 *
316  info = 0
317  IF( .NOT.rightv .AND. .NOT.leftv ) THEN
318  info = -1
319  ELSE IF( .NOT.fromqr .AND. .NOT.lsame( eigsrc, 'N' ) ) THEN
320  info = -2
321  ELSE IF( .NOT.noinit .AND. .NOT.lsame( initv, 'U' ) ) THEN
322  info = -3
323  ELSE IF( n.LT.0 ) THEN
324  info = -5
325  ELSE IF( ldh.LT.max( 1, n ) ) THEN
326  info = -7
327  ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
328  info = -10
329  ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
330  info = -12
331  ELSE IF( mm.LT.m ) THEN
332  info = -13
333  END IF
334  IF( info.NE.0 ) THEN
335  CALL xerbla( 'ZHSEIN', -info )
336  return
337  END IF
338 *
339 * Quick return if possible.
340 *
341  IF( n.EQ.0 )
342  $ return
343 *
344 * Set machine-dependent constants.
345 *
346  unfl = dlamch( 'Safe minimum' )
347  ulp = dlamch( 'Precision' )
348  smlnum = unfl*( n / ulp )
349 *
350  ldwork = n
351 *
352  kl = 1
353  kln = 0
354  IF( fromqr ) THEN
355  kr = 0
356  ELSE
357  kr = n
358  END IF
359  ks = 1
360 *
361  DO 100 k = 1, n
362  IF( SELECT( k ) ) THEN
363 *
364 * Compute eigenvector(s) corresponding to W(K).
365 *
366  IF( fromqr ) THEN
367 *
368 * If affiliation of eigenvalues is known, check whether
369 * the matrix splits.
370 *
371 * Determine KL and KR such that 1 <= KL <= K <= KR <= N
372 * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
373 * KR = N).
374 *
375 * Then inverse iteration can be performed with the
376 * submatrix H(KL:N,KL:N) for a left eigenvector, and with
377 * the submatrix H(1:KR,1:KR) for a right eigenvector.
378 *
379  DO 20 i = k, kl + 1, -1
380  IF( h( i, i-1 ).EQ.zero )
381  $ go to 30
382  20 continue
383  30 continue
384  kl = i
385  IF( k.GT.kr ) THEN
386  DO 40 i = k, n - 1
387  IF( h( i+1, i ).EQ.zero )
388  $ go to 50
389  40 continue
390  50 continue
391  kr = i
392  END IF
393  END IF
394 *
395  IF( kl.NE.kln ) THEN
396  kln = kl
397 *
398 * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
399 * has not ben computed before.
400 *
401  hnorm = zlanhs( 'I', kr-kl+1, h( kl, kl ), ldh, rwork )
402  IF( hnorm.GT.rzero ) THEN
403  eps3 = hnorm*ulp
404  ELSE
405  eps3 = smlnum
406  END IF
407  END IF
408 *
409 * Perturb eigenvalue if it is close to any previous
410 * selected eigenvalues affiliated to the submatrix
411 * H(KL:KR,KL:KR). Close roots are modified by EPS3.
412 *
413  wk = w( k )
414  60 continue
415  DO 70 i = k - 1, kl, -1
416  IF( SELECT( i ) .AND. cabs1( w( i )-wk ).LT.eps3 ) THEN
417  wk = wk + eps3
418  go to 60
419  END IF
420  70 continue
421  w( k ) = wk
422 *
423  IF( leftv ) THEN
424 *
425 * Compute left eigenvector.
426 *
427  CALL zlaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,
428  $ wk, vl( kl, ks ), work, ldwork, rwork, eps3,
429  $ smlnum, iinfo )
430  IF( iinfo.GT.0 ) THEN
431  info = info + 1
432  ifaill( ks ) = k
433  ELSE
434  ifaill( ks ) = 0
435  END IF
436  DO 80 i = 1, kl - 1
437  vl( i, ks ) = zero
438  80 continue
439  END IF
440  IF( rightv ) THEN
441 *
442 * Compute right eigenvector.
443 *
444  CALL zlaein( .true., noinit, kr, h, ldh, wk, vr( 1, ks ),
445  $ work, ldwork, rwork, eps3, smlnum, iinfo )
446  IF( iinfo.GT.0 ) THEN
447  info = info + 1
448  ifailr( ks ) = k
449  ELSE
450  ifailr( ks ) = 0
451  END IF
452  DO 90 i = kr + 1, n
453  vr( i, ks ) = zero
454  90 continue
455  END IF
456  ks = ks + 1
457  END IF
458  100 continue
459 *
460  return
461 *
462 * End of ZHSEIN
463 *
464  END