LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
shsein.f
Go to the documentation of this file.
1 *> \brief \b SHSEIN
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SHSEIN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/shsein.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/shsein.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/shsein.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI,
22 * VL, LDVL, VR, LDVR, MM, M, WORK, 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 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
33 * $ WI( * ), WORK( * ), WR( * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> SHSEIN uses inverse iteration to find specified right and/or left
43 *> eigenvectors of a real upper Hessenberg matrix H.
44 *>
45 *> The right eigenvector x and the left eigenvector y of the matrix H
46 *> corresponding to an eigenvalue w are defined by:
47 *>
48 *> H * x = w * x, y**h * H = w * y**h
49 *>
50 *> where y**h denotes the conjugate transpose of the vector y.
51 *> \endverbatim
52 *
53 * Arguments:
54 * ==========
55 *
56 *> \param[in] SIDE
57 *> \verbatim
58 *> SIDE is CHARACTER*1
59 *> = 'R': compute right eigenvectors only;
60 *> = 'L': compute left eigenvectors only;
61 *> = 'B': compute both right and left eigenvectors.
62 *> \endverbatim
63 *>
64 *> \param[in] EIGSRC
65 *> \verbatim
66 *> EIGSRC is CHARACTER*1
67 *> Specifies the source of eigenvalues supplied in (WR,WI):
68 *> = 'Q': the eigenvalues were found using SHSEQR; thus, if
69 *> H has zero subdiagonal elements, and so is
70 *> block-triangular, then the j-th eigenvalue can be
71 *> assumed to be an eigenvalue of the block containing
72 *> the j-th row/column. This property allows SHSEIN to
73 *> perform inverse iteration on just one diagonal block.
74 *> = 'N': no assumptions are made on the correspondence
75 *> between eigenvalues and diagonal blocks. In this
76 *> case, SHSEIN must always perform inverse iteration
77 *> using the whole matrix H.
78 *> \endverbatim
79 *>
80 *> \param[in] INITV
81 *> \verbatim
82 *> INITV is CHARACTER*1
83 *> = 'N': no initial vectors are supplied;
84 *> = 'U': user-supplied initial vectors are stored in the arrays
85 *> VL and/or VR.
86 *> \endverbatim
87 *>
88 *> \param[in,out] SELECT
89 *> \verbatim
90 *> SELECT is LOGICAL array, dimension (N)
91 *> Specifies the eigenvectors to be computed. To select the
92 *> real eigenvector corresponding to a real eigenvalue WR(j),
93 *> SELECT(j) must be set to .TRUE.. To select the complex
94 *> eigenvector corresponding to a complex eigenvalue
95 *> (WR(j),WI(j)), with complex conjugate (WR(j+1),WI(j+1)),
96 *> either SELECT(j) or SELECT(j+1) or both must be set to
97 *> .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is
98 *> .FALSE..
99 *> \endverbatim
100 *>
101 *> \param[in] N
102 *> \verbatim
103 *> N is INTEGER
104 *> The order of the matrix H. N >= 0.
105 *> \endverbatim
106 *>
107 *> \param[in] H
108 *> \verbatim
109 *> H is REAL array, dimension (LDH,N)
110 *> The upper Hessenberg matrix H.
111 *> \endverbatim
112 *>
113 *> \param[in] LDH
114 *> \verbatim
115 *> LDH is INTEGER
116 *> The leading dimension of the array H. LDH >= max(1,N).
117 *> \endverbatim
118 *>
119 *> \param[in,out] WR
120 *> \verbatim
121 *> WR is REAL array, dimension (N)
122 *> \endverbatim
123 *>
124 *> \param[in] WI
125 *> \verbatim
126 *> WI is REAL array, dimension (N)
127 *>
128 *> On entry, the real and imaginary parts of the eigenvalues of
129 *> H; a complex conjugate pair of eigenvalues must be stored in
130 *> consecutive elements of WR and WI.
131 *> On exit, WR may have been altered since close eigenvalues
132 *> are perturbed slightly in searching for independent
133 *> eigenvectors.
134 *> \endverbatim
135 *>
136 *> \param[in,out] VL
137 *> \verbatim
138 *> VL is REAL array, dimension (LDVL,MM)
139 *> On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
140 *> contain starting vectors for the inverse iteration for the
141 *> left eigenvectors; the starting vector for each eigenvector
142 *> must be in the same column(s) in which the eigenvector will
143 *> be stored.
144 *> On exit, if SIDE = 'L' or 'B', the left eigenvectors
145 *> specified by SELECT will be stored consecutively in the
146 *> columns of VL, in the same order as their eigenvalues. A
147 *> complex eigenvector corresponding to a complex eigenvalue is
148 *> stored in two consecutive columns, the first holding the real
149 *> part and the second the imaginary part.
150 *> If SIDE = 'R', VL is not referenced.
151 *> \endverbatim
152 *>
153 *> \param[in] LDVL
154 *> \verbatim
155 *> LDVL is INTEGER
156 *> The leading dimension of the array VL.
157 *> LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
158 *> \endverbatim
159 *>
160 *> \param[in,out] VR
161 *> \verbatim
162 *> VR is REAL array, dimension (LDVR,MM)
163 *> On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
164 *> contain starting vectors for the inverse iteration for the
165 *> right eigenvectors; the starting vector for each eigenvector
166 *> must be in the same column(s) in which the eigenvector will
167 *> be stored.
168 *> On exit, if SIDE = 'R' or 'B', the right eigenvectors
169 *> specified by SELECT will be stored consecutively in the
170 *> columns of VR, in the same order as their eigenvalues. A
171 *> complex eigenvector corresponding to a complex eigenvalue is
172 *> stored in two consecutive columns, the first holding the real
173 *> part and the second the imaginary part.
174 *> If SIDE = 'L', VR is not referenced.
175 *> \endverbatim
176 *>
177 *> \param[in] LDVR
178 *> \verbatim
179 *> LDVR is INTEGER
180 *> The leading dimension of the array VR.
181 *> LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
182 *> \endverbatim
183 *>
184 *> \param[in] MM
185 *> \verbatim
186 *> MM is INTEGER
187 *> The number of columns in the arrays VL and/or VR. MM >= M.
188 *> \endverbatim
189 *>
190 *> \param[out] M
191 *> \verbatim
192 *> M is INTEGER
193 *> The number of columns in the arrays VL and/or VR required to
194 *> store the eigenvectors; each selected real eigenvector
195 *> occupies one column and each selected complex eigenvector
196 *> occupies two columns.
197 *> \endverbatim
198 *>
199 *> \param[out] WORK
200 *> \verbatim
201 *> WORK is REAL array, dimension ((N+2)*N)
202 *> \endverbatim
203 *>
204 *> \param[out] IFAILL
205 *> \verbatim
206 *> IFAILL is INTEGER array, dimension (MM)
207 *> If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
208 *> eigenvector in the i-th column of VL (corresponding to the
209 *> eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
210 *> eigenvector converged satisfactorily. If the i-th and (i+1)th
211 *> columns of VL hold a complex eigenvector, then IFAILL(i) and
212 *> IFAILL(i+1) are set to the same value.
213 *> If SIDE = 'R', IFAILL is not referenced.
214 *> \endverbatim
215 *>
216 *> \param[out] IFAILR
217 *> \verbatim
218 *> IFAILR is INTEGER array, dimension (MM)
219 *> If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
220 *> eigenvector in the i-th column of VR (corresponding to the
221 *> eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
222 *> eigenvector converged satisfactorily. If the i-th and (i+1)th
223 *> columns of VR hold a complex eigenvector, then IFAILR(i) and
224 *> IFAILR(i+1) are set to the same value.
225 *> If SIDE = 'L', IFAILR is not referenced.
226 *> \endverbatim
227 *>
228 *> \param[out] INFO
229 *> \verbatim
230 *> INFO is INTEGER
231 *> = 0: successful exit
232 *> < 0: if INFO = -i, the i-th argument had an illegal value
233 *> > 0: if INFO = i, i is the number of eigenvectors which
234 *> failed to converge; see IFAILL and IFAILR for further
235 *> details.
236 *> \endverbatim
237 *
238 * Authors:
239 * ========
240 *
241 *> \author Univ. of Tennessee
242 *> \author Univ. of California Berkeley
243 *> \author Univ. of Colorado Denver
244 *> \author NAG Ltd.
245 *
246 *> \date November 2011
247 *
248 *> \ingroup realOTHERcomputational
249 *
250 *> \par Further Details:
251 * =====================
252 *>
253 *> \verbatim
254 *>
255 *> Each eigenvector is normalized so that the element of largest
256 *> magnitude has magnitude 1; here the magnitude of a complex number
257 *> (x,y) is taken to be |x|+|y|.
258 *> \endverbatim
259 *>
260 * =====================================================================
261  SUBROUTINE shsein( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI,
262  $ vl, ldvl, vr, ldvr, mm, m, work, ifaill,
263  $ ifailr, info )
264 *
265 * -- LAPACK computational routine (version 3.4.0) --
266 * -- LAPACK is a software package provided by Univ. of Tennessee, --
267 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268 * November 2011
269 *
270 * .. Scalar Arguments ..
271  CHARACTER eigsrc, initv, side
272  INTEGER info, ldh, ldvl, ldvr, m, mm, n
273 * ..
274 * .. Array Arguments ..
275  LOGICAL select( * )
276  INTEGER ifaill( * ), ifailr( * )
277  REAL h( ldh, * ), vl( ldvl, * ), vr( ldvr, * ),
278  $ wi( * ), work( * ), wr( * )
279 * ..
280 *
281 * =====================================================================
282 *
283 * .. Parameters ..
284  REAL zero, one
285  parameter( zero = 0.0e+0, one = 1.0e+0 )
286 * ..
287 * .. Local Scalars ..
288  LOGICAL bothv, fromqr, leftv, noinit, pair, rightv
289  INTEGER i, iinfo, k, kl, kln, kr, ksi, ksr, ldwork
290  REAL bignum, eps3, hnorm, smlnum, ulp, unfl, wki,
291  $ wkr
292 * ..
293 * .. External Functions ..
294  LOGICAL lsame
295  REAL slamch, slanhs
296  EXTERNAL lsame, slamch, slanhs
297 * ..
298 * .. External Subroutines ..
299  EXTERNAL slaein, xerbla
300 * ..
301 * .. Intrinsic Functions ..
302  INTRINSIC abs, max
303 * ..
304 * .. Executable Statements ..
305 *
306 * Decode and test the input parameters.
307 *
308  bothv = lsame( side, 'B' )
309  rightv = lsame( side, 'R' ) .OR. bothv
310  leftv = lsame( side, 'L' ) .OR. bothv
311 *
312  fromqr = lsame( eigsrc, 'Q' )
313 *
314  noinit = lsame( initv, 'N' )
315 *
316 * Set M to the number of columns required to store the selected
317 * eigenvectors, and standardize the array SELECT.
318 *
319  m = 0
320  pair = .false.
321  DO 10 k = 1, n
322  IF( pair ) THEN
323  pair = .false.
324  SELECT( k ) = .false.
325  ELSE
326  IF( wi( k ).EQ.zero ) THEN
327  IF( SELECT( k ) )
328  $ m = m + 1
329  ELSE
330  pair = .true.
331  IF( SELECT( k ) .OR. SELECT( k+1 ) ) THEN
332  SELECT( k ) = .true.
333  m = m + 2
334  END IF
335  END IF
336  END IF
337  10 continue
338 *
339  info = 0
340  IF( .NOT.rightv .AND. .NOT.leftv ) THEN
341  info = -1
342  ELSE IF( .NOT.fromqr .AND. .NOT.lsame( eigsrc, 'N' ) ) THEN
343  info = -2
344  ELSE IF( .NOT.noinit .AND. .NOT.lsame( initv, 'U' ) ) THEN
345  info = -3
346  ELSE IF( n.LT.0 ) THEN
347  info = -5
348  ELSE IF( ldh.LT.max( 1, n ) ) THEN
349  info = -7
350  ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
351  info = -11
352  ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
353  info = -13
354  ELSE IF( mm.LT.m ) THEN
355  info = -14
356  END IF
357  IF( info.NE.0 ) THEN
358  CALL xerbla( 'SHSEIN', -info )
359  return
360  END IF
361 *
362 * Quick return if possible.
363 *
364  IF( n.EQ.0 )
365  $ return
366 *
367 * Set machine-dependent constants.
368 *
369  unfl = slamch( 'Safe minimum' )
370  ulp = slamch( 'Precision' )
371  smlnum = unfl*( n / ulp )
372  bignum = ( one-ulp ) / smlnum
373 *
374  ldwork = n + 1
375 *
376  kl = 1
377  kln = 0
378  IF( fromqr ) THEN
379  kr = 0
380  ELSE
381  kr = n
382  END IF
383  ksr = 1
384 *
385  DO 120 k = 1, n
386  IF( SELECT( k ) ) THEN
387 *
388 * Compute eigenvector(s) corresponding to W(K).
389 *
390  IF( fromqr ) THEN
391 *
392 * If affiliation of eigenvalues is known, check whether
393 * the matrix splits.
394 *
395 * Determine KL and KR such that 1 <= KL <= K <= KR <= N
396 * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
397 * KR = N).
398 *
399 * Then inverse iteration can be performed with the
400 * submatrix H(KL:N,KL:N) for a left eigenvector, and with
401 * the submatrix H(1:KR,1:KR) for a right eigenvector.
402 *
403  DO 20 i = k, kl + 1, -1
404  IF( h( i, i-1 ).EQ.zero )
405  $ go to 30
406  20 continue
407  30 continue
408  kl = i
409  IF( k.GT.kr ) THEN
410  DO 40 i = k, n - 1
411  IF( h( i+1, i ).EQ.zero )
412  $ go to 50
413  40 continue
414  50 continue
415  kr = i
416  END IF
417  END IF
418 *
419  IF( kl.NE.kln ) THEN
420  kln = kl
421 *
422 * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
423 * has not ben computed before.
424 *
425  hnorm = slanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work )
426  IF( hnorm.GT.zero ) THEN
427  eps3 = hnorm*ulp
428  ELSE
429  eps3 = smlnum
430  END IF
431  END IF
432 *
433 * Perturb eigenvalue if it is close to any previous
434 * selected eigenvalues affiliated to the submatrix
435 * H(KL:KR,KL:KR). Close roots are modified by EPS3.
436 *
437  wkr = wr( k )
438  wki = wi( k )
439  60 continue
440  DO 70 i = k - 1, kl, -1
441  IF( SELECT( i ) .AND. abs( wr( i )-wkr )+
442  $ abs( wi( i )-wki ).LT.eps3 ) THEN
443  wkr = wkr + eps3
444  go to 60
445  END IF
446  70 continue
447  wr( k ) = wkr
448 *
449  pair = wki.NE.zero
450  IF( pair ) THEN
451  ksi = ksr + 1
452  ELSE
453  ksi = ksr
454  END IF
455  IF( leftv ) THEN
456 *
457 * Compute left eigenvector.
458 *
459  CALL slaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,
460  $ wkr, wki, vl( kl, ksr ), vl( kl, ksi ),
461  $ work, ldwork, work( n*n+n+1 ), eps3, smlnum,
462  $ bignum, iinfo )
463  IF( iinfo.GT.0 ) THEN
464  IF( pair ) THEN
465  info = info + 2
466  ELSE
467  info = info + 1
468  END IF
469  ifaill( ksr ) = k
470  ifaill( ksi ) = k
471  ELSE
472  ifaill( ksr ) = 0
473  ifaill( ksi ) = 0
474  END IF
475  DO 80 i = 1, kl - 1
476  vl( i, ksr ) = zero
477  80 continue
478  IF( pair ) THEN
479  DO 90 i = 1, kl - 1
480  vl( i, ksi ) = zero
481  90 continue
482  END IF
483  END IF
484  IF( rightv ) THEN
485 *
486 * Compute right eigenvector.
487 *
488  CALL slaein( .true., noinit, kr, h, ldh, wkr, wki,
489  $ vr( 1, ksr ), vr( 1, ksi ), work, ldwork,
490  $ work( n*n+n+1 ), eps3, smlnum, bignum,
491  $ iinfo )
492  IF( iinfo.GT.0 ) THEN
493  IF( pair ) THEN
494  info = info + 2
495  ELSE
496  info = info + 1
497  END IF
498  ifailr( ksr ) = k
499  ifailr( ksi ) = k
500  ELSE
501  ifailr( ksr ) = 0
502  ifailr( ksi ) = 0
503  END IF
504  DO 100 i = kr + 1, n
505  vr( i, ksr ) = zero
506  100 continue
507  IF( pair ) THEN
508  DO 110 i = kr + 1, n
509  vr( i, ksi ) = zero
510  110 continue
511  END IF
512  END IF
513 *
514  IF( pair ) THEN
515  ksr = ksr + 2
516  ELSE
517  ksr = ksr + 1
518  END IF
519  END IF
520  120 continue
521 *
522  return
523 *
524 * End of SHSEIN
525 *
526  END