LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
chsein.f
Go to the documentation of this file.
1 *> \brief \b CHSEIN
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CHSEIN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chsein.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chsein.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chsein.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHSEIN( 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 * REAL RWORK( * )
33 * COMPLEX H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
34 * $ W( * ), WORK( * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> CHSEIN 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 CHSEQR; 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 CHSEIN 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, CHSEIN 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 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 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 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 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 array, dimension (N*N)
183 *> \endverbatim
184 *>
185 *> \param[out] RWORK
186 *> \verbatim
187 *> RWORK is REAL 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 complexOTHERcomputational
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 chsein( 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  REAL rwork( * )
260  COMPLEX h( ldh, * ), vl( ldvl, * ), vr( ldvr, * ),
261  $ w( * ), work( * )
262 * ..
263 *
264 * =====================================================================
265 *
266 * .. Parameters ..
267  COMPLEX zero
268  parameter( zero = ( 0.0e+0, 0.0e+0 ) )
269  REAL rzero
270  parameter( rzero = 0.0e+0 )
271 * ..
272 * .. Local Scalars ..
273  LOGICAL bothv, fromqr, leftv, noinit, rightv
274  INTEGER i, iinfo, k, kl, kln, kr, ks, ldwork
275  REAL eps3, hnorm, smlnum, ulp, unfl
276  COMPLEX cdum, wk
277 * ..
278 * .. External Functions ..
279  LOGICAL lsame
280  REAL clanhs, slamch
281  EXTERNAL lsame, clanhs, slamch
282 * ..
283 * .. External Subroutines ..
284  EXTERNAL claein, xerbla
285 * ..
286 * .. Intrinsic Functions ..
287  INTRINSIC abs, aimag, max, real
288 * ..
289 * .. Statement Functions ..
290  REAL cabs1
291 * ..
292 * .. Statement Function definitions ..
293  cabs1( cdum ) = abs( REAL( CDUM ) ) + abs( aimag( 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( 'CHSEIN', -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 = slamch( 'Safe minimum' )
347  ulp = slamch( '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 = clanhs( '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 claein( .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 claein( .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 CHSEIN
463 *
464  END