LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ shsein()

subroutine shsein ( character side,
character eigsrc,
character initv,
logical, dimension( * ) select,
integer n,
real, dimension( ldh, * ) h,
integer ldh,
real, dimension( * ) wr,
real, dimension( * ) wi,
real, dimension( ldvl, * ) vl,
integer ldvl,
real, dimension( ldvr, * ) vr,
integer ldvr,
integer mm,
integer m,
real, dimension( * ) work,
integer, dimension( * ) ifaill,
integer, dimension( * ) ifailr,
integer info )

SHSEIN

Download SHSEIN + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SHSEIN uses inverse iteration to find specified right and/or left
!> eigenvectors of a real upper Hessenberg matrix H.
!>
!> The right eigenvector x and the left eigenvector y of the matrix H
!> corresponding to an eigenvalue w are defined by:
!>
!>              H * x = w * x,     y**h * H = w * y**h
!>
!> where y**h denotes the conjugate transpose of the vector y.
!>
Parameters
[in]SIDE
!>          SIDE is CHARACTER*1
!>          = 'R': compute right eigenvectors only;
!>          = 'L': compute left eigenvectors only;
!>          = 'B': compute both right and left eigenvectors.
!>
[in]EIGSRC
!>          EIGSRC is CHARACTER*1
!>          Specifies the source of eigenvalues supplied in (WR,WI):
!>          = 'Q': the eigenvalues were found using SHSEQR; thus, if
!>                 H has zero subdiagonal elements, and so is
!>                 block-triangular, then the j-th eigenvalue can be
!>                 assumed to be an eigenvalue of the block containing
!>                 the j-th row/column.  This property allows SHSEIN to
!>                 perform inverse iteration on just one diagonal block.
!>          = 'N': no assumptions are made on the correspondence
!>                 between eigenvalues and diagonal blocks.  In this
!>                 case, SHSEIN must always perform inverse iteration
!>                 using the whole matrix H.
!>
[in]INITV
!>          INITV is CHARACTER*1
!>          = 'N': no initial vectors are supplied;
!>          = 'U': user-supplied initial vectors are stored in the arrays
!>                 VL and/or VR.
!>
[in,out]SELECT
!>          SELECT is LOGICAL array, dimension (N)
!>          Specifies the eigenvectors to be computed. To select the
!>          real eigenvector corresponding to a real eigenvalue WR(j),
!>          SELECT(j) must be set to .TRUE.. To select the complex
!>          eigenvector corresponding to a complex eigenvalue
!>          (WR(j),WI(j)), with complex conjugate (WR(j+1),WI(j+1)),
!>          either SELECT(j) or SELECT(j+1) or both must be set to
!>          .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is
!>          .FALSE..
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix H.  N >= 0.
!>
[in]H
!>          H is REAL array, dimension (LDH,N)
!>          The upper Hessenberg matrix H.
!>          If a NaN is detected in H, the routine will return with INFO=-6.
!>
[in]LDH
!>          LDH is INTEGER
!>          The leading dimension of the array H.  LDH >= max(1,N).
!>
[in,out]WR
!>          WR is REAL array, dimension (N)
!>
[in]WI
!>          WI is REAL array, dimension (N)
!>
!>          On entry, the real and imaginary parts of the eigenvalues of
!>          H; a complex conjugate pair of eigenvalues must be stored in
!>          consecutive elements of WR and WI.
!>          On exit, WR may have been altered since close eigenvalues
!>          are perturbed slightly in searching for independent
!>          eigenvectors.
!>
[in,out]VL
!>          VL is REAL array, dimension (LDVL,MM)
!>          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
!>          contain starting vectors for the inverse iteration for the
!>          left eigenvectors; the starting vector for each eigenvector
!>          must be in the same column(s) in which the eigenvector will
!>          be stored.
!>          On exit, if SIDE = 'L' or 'B', the left eigenvectors
!>          specified by SELECT will be stored consecutively in the
!>          columns of VL, in the same order as their eigenvalues. A
!>          complex eigenvector corresponding to a complex eigenvalue is
!>          stored in two consecutive columns, the first holding the real
!>          part and the second the imaginary part.
!>          If SIDE = 'R', VL is not referenced.
!>
[in]LDVL
!>          LDVL is INTEGER
!>          The leading dimension of the array VL.
!>          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
!>
[in,out]VR
!>          VR is REAL array, dimension (LDVR,MM)
!>          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
!>          contain starting vectors for the inverse iteration for the
!>          right eigenvectors; the starting vector for each eigenvector
!>          must be in the same column(s) in which the eigenvector will
!>          be stored.
!>          On exit, if SIDE = 'R' or 'B', the right eigenvectors
!>          specified by SELECT will be stored consecutively in the
!>          columns of VR, in the same order as their eigenvalues. A
!>          complex eigenvector corresponding to a complex eigenvalue is
!>          stored in two consecutive columns, the first holding the real
!>          part and the second the imaginary part.
!>          If SIDE = 'L', VR is not referenced.
!>
[in]LDVR
!>          LDVR is INTEGER
!>          The leading dimension of the array VR.
!>          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
!>
[in]MM
!>          MM is INTEGER
!>          The number of columns in the arrays VL and/or VR. MM >= M.
!>
[out]M
!>          M is INTEGER
!>          The number of columns in the arrays VL and/or VR required to
!>          store the eigenvectors; each selected real eigenvector
!>          occupies one column and each selected complex eigenvector
!>          occupies two columns.
!>
[out]WORK
!>          WORK is REAL array, dimension ((N+2)*N)
!>
[out]IFAILL
!>          IFAILL is INTEGER array, dimension (MM)
!>          If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
!>          eigenvector in the i-th column of VL (corresponding to the
!>          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
!>          eigenvector converged satisfactorily. If the i-th and (i+1)th
!>          columns of VL hold a complex eigenvector, then IFAILL(i) and
!>          IFAILL(i+1) are set to the same value.
!>          If SIDE = 'R', IFAILL is not referenced.
!>
[out]IFAILR
!>          IFAILR is INTEGER array, dimension (MM)
!>          If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
!>          eigenvector in the i-th column of VR (corresponding to the
!>          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
!>          eigenvector converged satisfactorily. If the i-th and (i+1)th
!>          columns of VR hold a complex eigenvector, then IFAILR(i) and
!>          IFAILR(i+1) are set to the same value.
!>          If SIDE = 'L', IFAILR is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, i is the number of eigenvectors which
!>                failed to converge; see IFAILL and IFAILR for further
!>                details.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Each eigenvector is normalized so that the element of largest
!>  magnitude has magnitude 1; here the magnitude of a complex number
!>  (x,y) is taken to be |x|+|y|.
!>

Definition at line 258 of file shsein.f.

262*
263* -- LAPACK computational routine --
264* -- LAPACK is a software package provided by Univ. of Tennessee, --
265* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
266*
267* .. Scalar Arguments ..
268 CHARACTER EIGSRC, INITV, SIDE
269 INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
270* ..
271* .. Array Arguments ..
272 LOGICAL SELECT( * )
273 INTEGER IFAILL( * ), IFAILR( * )
274 REAL H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
275 $ WI( * ), WORK( * ), WR( * )
276* ..
277*
278* =====================================================================
279*
280* .. Parameters ..
281 REAL ZERO, ONE
282 parameter( zero = 0.0e+0, one = 1.0e+0 )
283* ..
284* .. Local Scalars ..
285 LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, PAIR, RIGHTV
286 INTEGER I, IINFO, K, KL, KLN, KR, KSI, KSR, LDWORK
287 REAL BIGNUM, EPS3, HNORM, SMLNUM, ULP, UNFL, WKI,
288 $ WKR
289* ..
290* .. External Functions ..
291 LOGICAL LSAME, SISNAN
292 REAL SLAMCH, SLANHS
293 EXTERNAL lsame, slamch, slanhs, sisnan
294* ..
295* .. External Subroutines ..
296 EXTERNAL slaein, xerbla
297* ..
298* .. Intrinsic Functions ..
299 INTRINSIC abs, max
300* ..
301* .. Executable Statements ..
302*
303* Decode and test the input parameters.
304*
305 bothv = lsame( side, 'B' )
306 rightv = lsame( side, 'R' ) .OR. bothv
307 leftv = lsame( side, 'L' ) .OR. bothv
308*
309 fromqr = lsame( eigsrc, 'Q' )
310*
311 noinit = lsame( initv, 'N' )
312*
313* Set M to the number of columns required to store the selected
314* eigenvectors, and standardize the array SELECT.
315*
316 m = 0
317 pair = .false.
318 DO 10 k = 1, n
319 IF( pair ) THEN
320 pair = .false.
321 SELECT( k ) = .false.
322 ELSE
323 IF( wi( k ).EQ.zero ) THEN
324 IF( SELECT( k ) )
325 $ m = m + 1
326 ELSE
327 pair = .true.
328 IF( SELECT( k ) .OR. SELECT( k+1 ) ) THEN
329 SELECT( k ) = .true.
330 m = m + 2
331 END IF
332 END IF
333 END IF
334 10 CONTINUE
335*
336 info = 0
337 IF( .NOT.rightv .AND. .NOT.leftv ) THEN
338 info = -1
339 ELSE IF( .NOT.fromqr .AND. .NOT.lsame( eigsrc, 'N' ) ) THEN
340 info = -2
341 ELSE IF( .NOT.noinit .AND. .NOT.lsame( initv, 'U' ) ) THEN
342 info = -3
343 ELSE IF( n.LT.0 ) THEN
344 info = -5
345 ELSE IF( ldh.LT.max( 1, n ) ) THEN
346 info = -7
347 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
348 info = -11
349 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
350 info = -13
351 ELSE IF( mm.LT.m ) THEN
352 info = -14
353 END IF
354 IF( info.NE.0 ) THEN
355 CALL xerbla( 'SHSEIN', -info )
356 RETURN
357 END IF
358*
359* Quick return if possible.
360*
361 IF( n.EQ.0 )
362 $ RETURN
363*
364* Set machine-dependent constants.
365*
366 unfl = slamch( 'Safe minimum' )
367 ulp = slamch( 'Precision' )
368 smlnum = unfl*( real( n ) / ulp )
369 bignum = ( one-ulp ) / smlnum
370*
371 ldwork = n + 1
372*
373 kl = 1
374 kln = 0
375 IF( fromqr ) THEN
376 kr = 0
377 ELSE
378 kr = n
379 END IF
380 ksr = 1
381*
382 DO 120 k = 1, n
383 IF( SELECT( k ) ) THEN
384*
385* Compute eigenvector(s) corresponding to W(K).
386*
387 IF( fromqr ) THEN
388*
389* If affiliation of eigenvalues is known, check whether
390* the matrix splits.
391*
392* Determine KL and KR such that 1 <= KL <= K <= KR <= N
393* and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
394* KR = N).
395*
396* Then inverse iteration can be performed with the
397* submatrix H(KL:N,KL:N) for a left eigenvector, and with
398* the submatrix H(1:KR,1:KR) for a right eigenvector.
399*
400 DO 20 i = k, kl + 1, -1
401 IF( h( i, i-1 ).EQ.zero )
402 $ GO TO 30
403 20 CONTINUE
404 30 CONTINUE
405 kl = i
406 IF( k.GT.kr ) THEN
407 DO 40 i = k, n - 1
408 IF( h( i+1, i ).EQ.zero )
409 $ GO TO 50
410 40 CONTINUE
411 50 CONTINUE
412 kr = i
413 END IF
414 END IF
415*
416 IF( kl.NE.kln ) THEN
417 kln = kl
418*
419* Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
420* has not ben computed before.
421*
422 hnorm = slanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work )
423 IF( sisnan( hnorm ) ) THEN
424 info = -6
425 RETURN
426 ELSE 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 ),
460 $ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sisnan
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:57
slaein
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:171
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slanhs
real function slanhs(norm, n, a, lda, work)
SLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slanhs.f:106
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: