LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> 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*> \ingroup hsein
248*
249*> \par Further Details:
250* =====================
251*>
252*> \verbatim
253*>
254*> Each eigenvector is normalized so that the element of largest
255*> magnitude has magnitude 1; here the magnitude of a complex number
256*> (x,y) is taken to be |x|+|y|.
257*> \endverbatim
258*>
259* =====================================================================
260 SUBROUTINE shsein( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI,
261 $ VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL,
262 $ IFAILR, INFO )
263*
264* -- LAPACK computational routine --
265* -- LAPACK is a software package provided by Univ. of Tennessee, --
266* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
267*
268* .. Scalar Arguments ..
269 CHARACTER EIGSRC, INITV, SIDE
270 INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
271* ..
272* .. Array Arguments ..
273 LOGICAL SELECT( * )
274 INTEGER IFAILL( * ), IFAILR( * )
275 REAL H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
276 $ wi( * ), work( * ), wr( * )
277* ..
278*
279* =====================================================================
280*
281* .. Parameters ..
282 REAL ZERO, ONE
283 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
284* ..
285* .. Local Scalars ..
286 LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, PAIR, RIGHTV
287 INTEGER I, IINFO, K, KL, KLN, KR, KSI, KSR, LDWORK
288 REAL BIGNUM, EPS3, HNORM, SMLNUM, ULP, UNFL, WKI,
289 $ wkr
290* ..
291* .. External Functions ..
292 LOGICAL LSAME, SISNAN
293 REAL SLAMCH, SLANHS
294 EXTERNAL lsame, slamch, slanhs, sisnan
295* ..
296* .. External Subroutines ..
297 EXTERNAL slaein, xerbla
298* ..
299* .. Intrinsic Functions ..
300 INTRINSIC abs, max
301* ..
302* .. Executable Statements ..
303*
304* Decode and test the input parameters.
305*
306 bothv = lsame( side, 'B' )
307 rightv = lsame( side, 'R' ) .OR. bothv
308 leftv = lsame( side, 'L' ) .OR. bothv
309*
310 fromqr = lsame( eigsrc, 'Q' )
311*
312 noinit = lsame( initv, 'N' )
313*
314* Set M to the number of columns required to store the selected
315* eigenvectors, and standardize the array SELECT.
316*
317 m = 0
318 pair = .false.
319 DO 10 k = 1, n
320 IF( pair ) THEN
321 pair = .false.
322 SELECT( k ) = .false.
323 ELSE
324 IF( wi( k ).EQ.zero ) THEN
325 IF( SELECT( k ) )
326 $ m = m + 1
327 ELSE
328 pair = .true.
329 IF( SELECT( k ) .OR. SELECT( k+1 ) ) THEN
330 SELECT( k ) = .true.
331 m = m + 2
332 END IF
333 END IF
334 END IF
335 10 CONTINUE
336*
337 info = 0
338 IF( .NOT.rightv .AND. .NOT.leftv ) THEN
339 info = -1
340 ELSE IF( .NOT.fromqr .AND. .NOT.lsame( eigsrc, 'N' ) ) THEN
341 info = -2
342 ELSE IF( .NOT.noinit .AND. .NOT.lsame( initv, 'U' ) ) THEN
343 info = -3
344 ELSE IF( n.LT.0 ) THEN
345 info = -5
346 ELSE IF( ldh.LT.max( 1, n ) ) THEN
347 info = -7
348 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
349 info = -11
350 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
351 info = -13
352 ELSE IF( mm.LT.m ) THEN
353 info = -14
354 END IF
355 IF( info.NE.0 ) THEN
356 CALL xerbla( 'SHSEIN', -info )
357 RETURN
358 END IF
359*
360* Quick return if possible.
361*
362 IF( n.EQ.0 )
363 $ RETURN
364*
365* Set machine-dependent constants.
366*
367 unfl = slamch( 'Safe minimum' )
368 ulp = slamch( 'Precision' )
369 smlnum = unfl*( n / ulp )
370 bignum = ( one-ulp ) / smlnum
371*
372 ldwork = n + 1
373*
374 kl = 1
375 kln = 0
376 IF( fromqr ) THEN
377 kr = 0
378 ELSE
379 kr = n
380 END IF
381 ksr = 1
382*
383 DO 120 k = 1, n
384 IF( SELECT( k ) ) THEN
385*
386* Compute eigenvector(s) corresponding to W(K).
387*
388 IF( fromqr ) THEN
389*
390* If affiliation of eigenvalues is known, check whether
391* the matrix splits.
392*
393* Determine KL and KR such that 1 <= KL <= K <= KR <= N
394* and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
395* KR = N).
396*
397* Then inverse iteration can be performed with the
398* submatrix H(KL:N,KL:N) for a left eigenvector, and with
399* the submatrix H(1:KR,1:KR) for a right eigenvector.
400*
401 DO 20 i = k, kl + 1, -1
402 IF( h( i, i-1 ).EQ.zero )
403 $ GO TO 30
404 20 CONTINUE
405 30 CONTINUE
406 kl = i
407 IF( k.GT.kr ) THEN
408 DO 40 i = k, n - 1
409 IF( h( i+1, i ).EQ.zero )
410 $ GO TO 50
411 40 CONTINUE
412 50 CONTINUE
413 kr = i
414 END IF
415 END IF
416*
417 IF( kl.NE.kln ) THEN
418 kln = kl
419*
420* Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
421* has not ben computed before.
422*
423 hnorm = slanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work )
424 IF( sisnan( hnorm ) ) THEN
425 info = -6
426 RETURN
427 ELSE IF( hnorm.GT.zero ) THEN
428 eps3 = hnorm*ulp
429 ELSE
430 eps3 = smlnum
431 END IF
432 END IF
433*
434* Perturb eigenvalue if it is close to any previous
435* selected eigenvalues affiliated to the submatrix
436* H(KL:KR,KL:KR). Close roots are modified by EPS3.
437*
438 wkr = wr( k )
439 wki = wi( k )
440 60 CONTINUE
441 DO 70 i = k - 1, kl, -1
442 IF( SELECT( i ) .AND. abs( wr( i )-wkr )+
443 $ abs( wi( i )-wki ).LT.eps3 ) THEN
444 wkr = wkr + eps3
445 GO TO 60
446 END IF
447 70 CONTINUE
448 wr( k ) = wkr
449*
450 pair = wki.NE.zero
451 IF( pair ) THEN
452 ksi = ksr + 1
453 ELSE
454 ksi = ksr
455 END IF
456 IF( leftv ) THEN
457*
458* Compute left eigenvector.
459*
460 CALL slaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,
461 $ wkr, wki, vl( kl, ksr ), vl( kl, ksi ),
462 $ work, ldwork, work( n*n+n+1 ), eps3, smlnum,
463 $ bignum, iinfo )
464 IF( iinfo.GT.0 ) THEN
465 IF( pair ) THEN
466 info = info + 2
467 ELSE
468 info = info + 1
469 END IF
470 ifaill( ksr ) = k
471 ifaill( ksi ) = k
472 ELSE
473 ifaill( ksr ) = 0
474 ifaill( ksi ) = 0
475 END IF
476 DO 80 i = 1, kl - 1
477 vl( i, ksr ) = zero
478 80 CONTINUE
479 IF( pair ) THEN
480 DO 90 i = 1, kl - 1
481 vl( i, ksi ) = zero
482 90 CONTINUE
483 END IF
484 END IF
485 IF( rightv ) THEN
486*
487* Compute right eigenvector.
488*
489 CALL slaein( .true., noinit, kr, h, ldh, wkr, wki,
490 $ vr( 1, ksr ), vr( 1, ksi ), work, ldwork,
491 $ work( n*n+n+1 ), eps3, smlnum, bignum,
492 $ iinfo )
493 IF( iinfo.GT.0 ) THEN
494 IF( pair ) THEN
495 info = info + 2
496 ELSE
497 info = info + 1
498 END IF
499 ifailr( ksr ) = k
500 ifailr( ksi ) = k
501 ELSE
502 ifailr( ksr ) = 0
503 ifailr( ksi ) = 0
504 END IF
505 DO 100 i = kr + 1, n
506 vr( i, ksr ) = zero
507 100 CONTINUE
508 IF( pair ) THEN
509 DO 110 i = kr + 1, n
510 vr( i, ksi ) = zero
511 110 CONTINUE
512 END IF
513 END IF
514*
515 IF( pair ) THEN
516 ksr = ksr + 2
517 ELSE
518 ksr = ksr + 1
519 END IF
520 END IF
521 120 CONTINUE
522*
523 RETURN
524*
525* End of SHSEIN
526*
527 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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:263
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:172