LAPACK 3.12.0
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 260 of file shsein.f.

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*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:59
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
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:108
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: