LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
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
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 *> If a NaN is detected in H, the routine will return with INFO=-6.
112 *> \endverbatim
113 *>
114 *> \param[in] LDH
115 *> \verbatim
116 *> LDH is INTEGER
117 *> The leading dimension of the array H. LDH >= max(1,N).
118 *> \endverbatim
119 *>
120 *> \param[in,out] WR
121 *> \verbatim
122 *> WR is REAL array, dimension (N)
123 *> \endverbatim
124 *>
125 *> \param[in] WI
126 *> \verbatim
127 *> WI is REAL array, dimension (N)
128 *>
129 *> On entry, the real and imaginary parts of the eigenvalues of
130 *> H; a complex conjugate pair of eigenvalues must be stored in
131 *> consecutive elements of WR and WI.
132 *> On exit, WR may have been altered since close eigenvalues
133 *> are perturbed slightly in searching for independent
134 *> eigenvectors.
135 *> \endverbatim
136 *>
137 *> \param[in,out] VL
138 *> \verbatim
139 *> VL is REAL array, dimension (LDVL,MM)
140 *> On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
141 *> contain starting vectors for the inverse iteration for the
142 *> left eigenvectors; the starting vector for each eigenvector
143 *> must be in the same column(s) in which the eigenvector will
144 *> be stored.
145 *> On exit, if SIDE = 'L' or 'B', the left eigenvectors
146 *> specified by SELECT will be stored consecutively in the
147 *> columns of VL, in the same order as their eigenvalues. A
148 *> complex eigenvector corresponding to a complex eigenvalue is
149 *> stored in two consecutive columns, the first holding the real
150 *> part and the second the imaginary part.
151 *> If SIDE = 'R', VL is not referenced.
152 *> \endverbatim
153 *>
154 *> \param[in] LDVL
155 *> \verbatim
156 *> LDVL is INTEGER
157 *> The leading dimension of the array VL.
158 *> LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
159 *> \endverbatim
160 *>
161 *> \param[in,out] VR
162 *> \verbatim
163 *> VR is REAL array, dimension (LDVR,MM)
164 *> On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
165 *> contain starting vectors for the inverse iteration for the
166 *> right eigenvectors; the starting vector for each eigenvector
167 *> must be in the same column(s) in which the eigenvector will
168 *> be stored.
169 *> On exit, if SIDE = 'R' or 'B', the right eigenvectors
170 *> specified by SELECT will be stored consecutively in the
171 *> columns of VR, in the same order as their eigenvalues. A
172 *> complex eigenvector corresponding to a complex eigenvalue is
173 *> stored in two consecutive columns, the first holding the real
174 *> part and the second the imaginary part.
175 *> If SIDE = 'L', VR is not referenced.
176 *> \endverbatim
177 *>
178 *> \param[in] LDVR
179 *> \verbatim
180 *> LDVR is INTEGER
181 *> The leading dimension of the array VR.
182 *> LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
183 *> \endverbatim
184 *>
185 *> \param[in] MM
186 *> \verbatim
187 *> MM is INTEGER
188 *> The number of columns in the arrays VL and/or VR. MM >= M.
189 *> \endverbatim
190 *>
191 *> \param[out] M
192 *> \verbatim
193 *> M is INTEGER
194 *> The number of columns in the arrays VL and/or VR required to
195 *> store the eigenvectors; each selected real eigenvector
196 *> occupies one column and each selected complex eigenvector
197 *> occupies two columns.
198 *> \endverbatim
199 *>
200 *> \param[out] WORK
201 *> \verbatim
202 *> WORK is REAL array, dimension ((N+2)*N)
203 *> \endverbatim
204 *>
205 *> \param[out] IFAILL
206 *> \verbatim
207 *> IFAILL is INTEGER array, dimension (MM)
208 *> If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
209 *> eigenvector in the i-th column of VL (corresponding to the
210 *> eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
211 *> eigenvector converged satisfactorily. If the i-th and (i+1)th
212 *> columns of VL hold a complex eigenvector, then IFAILL(i) and
213 *> IFAILL(i+1) are set to the same value.
214 *> If SIDE = 'R', IFAILL is not referenced.
215 *> \endverbatim
216 *>
217 *> \param[out] IFAILR
218 *> \verbatim
219 *> IFAILR is INTEGER array, dimension (MM)
220 *> If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
221 *> eigenvector in the i-th column of VR (corresponding to the
222 *> eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
223 *> eigenvector converged satisfactorily. If the i-th and (i+1)th
224 *> columns of VR hold a complex eigenvector, then IFAILR(i) and
225 *> IFAILR(i+1) are set to the same value.
226 *> If SIDE = 'L', IFAILR is not referenced.
227 *> \endverbatim
228 *>
229 *> \param[out] INFO
230 *> \verbatim
231 *> INFO is INTEGER
232 *> = 0: successful exit
233 *> < 0: if INFO = -i, the i-th argument had an illegal value
234 *> > 0: if INFO = i, i is the number of eigenvectors which
235 *> failed to converge; see IFAILL and IFAILR for further
236 *> details.
237 *> \endverbatim
238 *
239 * Authors:
240 * ========
241 *
242 *> \author Univ. of Tennessee
243 *> \author Univ. of California Berkeley
244 *> \author Univ. of Colorado Denver
245 *> \author NAG Ltd.
246 *
247 *> \date November 2013
248 *
249 *> \ingroup realOTHERcomputational
250 *
251 *> \par Further Details:
252 * =====================
253 *>
254 *> \verbatim
255 *>
256 *> Each eigenvector is normalized so that the element of largest
257 *> magnitude has magnitude 1; here the magnitude of a complex number
258 *> (x,y) is taken to be |x|+|y|.
259 *> \endverbatim
260 *>
261 * =====================================================================
262  SUBROUTINE shsein( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI,
263  \$ vl, ldvl, vr, ldvr, mm, m, work, ifaill,
264  \$ ifailr, info )
265 *
266 * -- LAPACK computational routine (version 3.5.0) --
267 * -- LAPACK is a software package provided by Univ. of Tennessee, --
268 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269 * November 2013
270 *
271 * .. Scalar Arguments ..
272  CHARACTER EIGSRC, INITV, SIDE
273  INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
274 * ..
275 * .. Array Arguments ..
276  LOGICAL SELECT( * )
277  INTEGER IFAILL( * ), IFAILR( * )
278  REAL H( ldh, * ), VL( ldvl, * ), VR( ldvr, * ),
279  \$ wi( * ), work( * ), wr( * )
280 * ..
281 *
282 * =====================================================================
283 *
284 * .. Parameters ..
285  REAL ZERO, ONE
286  parameter ( zero = 0.0e+0, one = 1.0e+0 )
287 * ..
288 * .. Local Scalars ..
289  LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, PAIR, RIGHTV
290  INTEGER I, IINFO, K, KL, KLN, KR, KSI, KSR, LDWORK
291  REAL BIGNUM, EPS3, HNORM, SMLNUM, ULP, UNFL, WKI,
292  \$ wkr
293 * ..
294 * .. External Functions ..
295  LOGICAL LSAME, SISNAN
296  REAL SLAMCH, SLANHS
297  EXTERNAL lsame, slamch, slanhs, sisnan
298 * ..
299 * .. External Subroutines ..
300  EXTERNAL slaein, xerbla
301 * ..
302 * .. Intrinsic Functions ..
303  INTRINSIC abs, max
304 * ..
305 * .. Executable Statements ..
306 *
307 * Decode and test the input parameters.
308 *
309  bothv = lsame( side, 'B' )
310  rightv = lsame( side, 'R' ) .OR. bothv
311  leftv = lsame( side, 'L' ) .OR. bothv
312 *
313  fromqr = lsame( eigsrc, 'Q' )
314 *
315  noinit = lsame( initv, 'N' )
316 *
317 * Set M to the number of columns required to store the selected
318 * eigenvectors, and standardize the array SELECT.
319 *
320  m = 0
321  pair = .false.
322  DO 10 k = 1, n
323  IF( pair ) THEN
324  pair = .false.
325  SELECT( k ) = .false.
326  ELSE
327  IF( wi( k ).EQ.zero ) THEN
328  IF( SELECT( k ) )
329  \$ m = m + 1
330  ELSE
331  pair = .true.
332  IF( SELECT( k ) .OR. SELECT( k+1 ) ) THEN
333  SELECT( k ) = .true.
334  m = m + 2
335  END IF
336  END IF
337  END IF
338  10 CONTINUE
339 *
340  info = 0
341  IF( .NOT.rightv .AND. .NOT.leftv ) THEN
342  info = -1
343  ELSE IF( .NOT.fromqr .AND. .NOT.lsame( eigsrc, 'N' ) ) THEN
344  info = -2
345  ELSE IF( .NOT.noinit .AND. .NOT.lsame( initv, 'U' ) ) THEN
346  info = -3
347  ELSE IF( n.LT.0 ) THEN
348  info = -5
349  ELSE IF( ldh.LT.max( 1, n ) ) THEN
350  info = -7
351  ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
352  info = -11
353  ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
354  info = -13
355  ELSE IF( mm.LT.m ) THEN
356  info = -14
357  END IF
358  IF( info.NE.0 ) THEN
359  CALL xerbla( 'SHSEIN', -info )
360  RETURN
361  END IF
362 *
363 * Quick return if possible.
364 *
365  IF( n.EQ.0 )
366  \$ RETURN
367 *
368 * Set machine-dependent constants.
369 *
370  unfl = slamch( 'Safe minimum' )
371  ulp = slamch( 'Precision' )
372  smlnum = unfl*( n / ulp )
373  bignum = ( one-ulp ) / smlnum
374 *
375  ldwork = n + 1
376 *
377  kl = 1
378  kln = 0
379  IF( fromqr ) THEN
380  kr = 0
381  ELSE
382  kr = n
383  END IF
384  ksr = 1
385 *
386  DO 120 k = 1, n
387  IF( SELECT( k ) ) THEN
388 *
389 * Compute eigenvector(s) corresponding to W(K).
390 *
391  IF( fromqr ) THEN
392 *
393 * If affiliation of eigenvalues is known, check whether
394 * the matrix splits.
395 *
396 * Determine KL and KR such that 1 <= KL <= K <= KR <= N
397 * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
398 * KR = N).
399 *
400 * Then inverse iteration can be performed with the
401 * submatrix H(KL:N,KL:N) for a left eigenvector, and with
402 * the submatrix H(1:KR,1:KR) for a right eigenvector.
403 *
404  DO 20 i = k, kl + 1, -1
405  IF( h( i, i-1 ).EQ.zero )
406  \$ GO TO 30
407  20 CONTINUE
408  30 CONTINUE
409  kl = i
410  IF( k.GT.kr ) THEN
411  DO 40 i = k, n - 1
412  IF( h( i+1, i ).EQ.zero )
413  \$ GO TO 50
414  40 CONTINUE
415  50 CONTINUE
416  kr = i
417  END IF
418  END IF
419 *
420  IF( kl.NE.kln ) THEN
421  kln = kl
422 *
423 * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
424 * has not ben computed before.
425 *
426  hnorm = slanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work )
427  IF( sisnan( hnorm ) ) THEN
428  info = -6
429  RETURN
430  ELSE IF( hnorm.GT.zero ) THEN
431  eps3 = hnorm*ulp
432  ELSE
433  eps3 = smlnum
434  END IF
435  END IF
436 *
437 * Perturb eigenvalue if it is close to any previous
438 * selected eigenvalues affiliated to the submatrix
439 * H(KL:KR,KL:KR). Close roots are modified by EPS3.
440 *
441  wkr = wr( k )
442  wki = wi( k )
443  60 CONTINUE
444  DO 70 i = k - 1, kl, -1
445  IF( SELECT( i ) .AND. abs( wr( i )-wkr )+
446  \$ abs( wi( i )-wki ).LT.eps3 ) THEN
447  wkr = wkr + eps3
448  GO TO 60
449  END IF
450  70 CONTINUE
451  wr( k ) = wkr
452 *
453  pair = wki.NE.zero
454  IF( pair ) THEN
455  ksi = ksr + 1
456  ELSE
457  ksi = ksr
458  END IF
459  IF( leftv ) THEN
460 *
461 * Compute left eigenvector.
462 *
463  CALL slaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,
464  \$ wkr, wki, vl( kl, ksr ), vl( kl, ksi ),
465  \$ work, ldwork, work( n*n+n+1 ), eps3, smlnum,
466  \$ bignum, iinfo )
467  IF( iinfo.GT.0 ) THEN
468  IF( pair ) THEN
469  info = info + 2
470  ELSE
471  info = info + 1
472  END IF
473  ifaill( ksr ) = k
474  ifaill( ksi ) = k
475  ELSE
476  ifaill( ksr ) = 0
477  ifaill( ksi ) = 0
478  END IF
479  DO 80 i = 1, kl - 1
480  vl( i, ksr ) = zero
481  80 CONTINUE
482  IF( pair ) THEN
483  DO 90 i = 1, kl - 1
484  vl( i, ksi ) = zero
485  90 CONTINUE
486  END IF
487  END IF
488  IF( rightv ) THEN
489 *
490 * Compute right eigenvector.
491 *
492  CALL slaein( .true., noinit, kr, h, ldh, wkr, wki,
493  \$ vr( 1, ksr ), vr( 1, ksi ), work, ldwork,
494  \$ work( n*n+n+1 ), eps3, smlnum, bignum,
495  \$ iinfo )
496  IF( iinfo.GT.0 ) THEN
497  IF( pair ) THEN
498  info = info + 2
499  ELSE
500  info = info + 1
501  END IF
502  ifailr( ksr ) = k
503  ifailr( ksi ) = k
504  ELSE
505  ifailr( ksr ) = 0
506  ifailr( ksi ) = 0
507  END IF
508  DO 100 i = kr + 1, n
509  vr( i, ksr ) = zero
510  100 CONTINUE
511  IF( pair ) THEN
512  DO 110 i = kr + 1, n
513  vr( i, ksi ) = zero
514  110 CONTINUE
515  END IF
516  END IF
517 *
518  IF( pair ) THEN
519  ksr = ksr + 2
520  ELSE
521  ksr = ksr + 1
522  END IF
523  END IF
524  120 CONTINUE
525 *
526  RETURN
527 *
528 * End of SHSEIN
529 *
530  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slaein(RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO)
SLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iterat...
Definition: slaein.f:174
subroutine shsein(SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI, VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL, IFAILR, INFO)
SHSEIN
Definition: shsein.f:265