LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zhsein ( character  SIDE,
character  EIGSRC,
character  INITV,
logical, dimension( * )  SELECT,
integer  N,
complex*16, dimension( ldh, * )  H,
integer  LDH,
complex*16, dimension( * )  W,
complex*16, dimension( ldvl, * )  VL,
integer  LDVL,
complex*16, dimension( ldvr, * )  VR,
integer  LDVR,
integer  MM,
integer  M,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer, dimension( * )  IFAILL,
integer, dimension( * )  IFAILR,
integer  INFO 
)

ZHSEIN

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

Purpose:
 ZHSEIN uses inverse iteration to find specified right and/or left
 eigenvectors of a complex 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 W:
          = 'Q': the eigenvalues were found using ZHSEQR; 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 ZHSEIN 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, ZHSEIN 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]SELECT
          SELECT is LOGICAL array, dimension (N)
          Specifies the eigenvectors to be computed. To select the
          eigenvector corresponding to the eigenvalue W(j),
          SELECT(j) must be set to .TRUE..
[in]N
          N is INTEGER
          The order of the matrix H.  N >= 0.
[in]H
          H is COMPLEX*16 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]W
          W is COMPLEX*16 array, dimension (N)
          On entry, the eigenvalues of H.
          On exit, the real parts of W may have been altered since
          close eigenvalues are perturbed slightly in searching for
          independent eigenvectors.
[in,out]VL
          VL is COMPLEX*16 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 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.
          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 COMPLEX*16 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 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.
          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 (= the number of .TRUE. elements in
          SELECT).
[out]WORK
          WORK is COMPLEX*16 array, dimension (N*N)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (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 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 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.
Date
November 2013
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 247 of file zhsein.f.

247 *
248 * -- LAPACK computational routine (version 3.5.0) --
249 * -- LAPACK is a software package provided by Univ. of Tennessee, --
250 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
251 * November 2013
252 *
253 * .. Scalar Arguments ..
254  CHARACTER eigsrc, initv, side
255  INTEGER info, ldh, ldvl, ldvr, m, mm, n
256 * ..
257 * .. Array Arguments ..
258  LOGICAL select( * )
259  INTEGER ifaill( * ), ifailr( * )
260  DOUBLE PRECISION rwork( * )
261  COMPLEX*16 h( ldh, * ), vl( ldvl, * ), vr( ldvr, * ),
262  $ w( * ), work( * )
263 * ..
264 *
265 * =====================================================================
266 *
267 * .. Parameters ..
268  COMPLEX*16 zero
269  parameter ( zero = ( 0.0d+0, 0.0d+0 ) )
270  DOUBLE PRECISION rzero
271  parameter ( rzero = 0.0d+0 )
272 * ..
273 * .. Local Scalars ..
274  LOGICAL bothv, fromqr, leftv, noinit, rightv
275  INTEGER i, iinfo, k, kl, kln, kr, ks, ldwork
276  DOUBLE PRECISION eps3, hnorm, smlnum, ulp, unfl
277  COMPLEX*16 cdum, wk
278 * ..
279 * .. External Functions ..
280  LOGICAL lsame, disnan
281  DOUBLE PRECISION dlamch, zlanhs
282  EXTERNAL lsame, dlamch, zlanhs, disnan
283 * ..
284 * .. External Subroutines ..
285  EXTERNAL xerbla, zlaein
286 * ..
287 * .. Intrinsic Functions ..
288  INTRINSIC abs, dble, dimag, max
289 * ..
290 * .. Statement Functions ..
291  DOUBLE PRECISION cabs1
292 * ..
293 * .. Statement Function definitions ..
294  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
295 * ..
296 * .. Executable Statements ..
297 *
298 * Decode and test the input parameters.
299 *
300  bothv = lsame( side, 'B' )
301  rightv = lsame( side, 'R' ) .OR. bothv
302  leftv = lsame( side, 'L' ) .OR. bothv
303 *
304  fromqr = lsame( eigsrc, 'Q' )
305 *
306  noinit = lsame( initv, 'N' )
307 *
308 * Set M to the number of columns required to store the selected
309 * eigenvectors.
310 *
311  m = 0
312  DO 10 k = 1, n
313  IF( SELECT( k ) )
314  $ m = m + 1
315  10 CONTINUE
316 *
317  info = 0
318  IF( .NOT.rightv .AND. .NOT.leftv ) THEN
319  info = -1
320  ELSE IF( .NOT.fromqr .AND. .NOT.lsame( eigsrc, 'N' ) ) THEN
321  info = -2
322  ELSE IF( .NOT.noinit .AND. .NOT.lsame( initv, 'U' ) ) THEN
323  info = -3
324  ELSE IF( n.LT.0 ) THEN
325  info = -5
326  ELSE IF( ldh.LT.max( 1, n ) ) THEN
327  info = -7
328  ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
329  info = -10
330  ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
331  info = -12
332  ELSE IF( mm.LT.m ) THEN
333  info = -13
334  END IF
335  IF( info.NE.0 ) THEN
336  CALL xerbla( 'ZHSEIN', -info )
337  RETURN
338  END IF
339 *
340 * Quick return if possible.
341 *
342  IF( n.EQ.0 )
343  $ RETURN
344 *
345 * Set machine-dependent constants.
346 *
347  unfl = dlamch( 'Safe minimum' )
348  ulp = dlamch( 'Precision' )
349  smlnum = unfl*( n / ulp )
350 *
351  ldwork = n
352 *
353  kl = 1
354  kln = 0
355  IF( fromqr ) THEN
356  kr = 0
357  ELSE
358  kr = n
359  END IF
360  ks = 1
361 *
362  DO 100 k = 1, n
363  IF( SELECT( k ) ) THEN
364 *
365 * Compute eigenvector(s) corresponding to W(K).
366 *
367  IF( fromqr ) THEN
368 *
369 * If affiliation of eigenvalues is known, check whether
370 * the matrix splits.
371 *
372 * Determine KL and KR such that 1 <= KL <= K <= KR <= N
373 * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
374 * KR = N).
375 *
376 * Then inverse iteration can be performed with the
377 * submatrix H(KL:N,KL:N) for a left eigenvector, and with
378 * the submatrix H(1:KR,1:KR) for a right eigenvector.
379 *
380  DO 20 i = k, kl + 1, -1
381  IF( h( i, i-1 ).EQ.zero )
382  $ GO TO 30
383  20 CONTINUE
384  30 CONTINUE
385  kl = i
386  IF( k.GT.kr ) THEN
387  DO 40 i = k, n - 1
388  IF( h( i+1, i ).EQ.zero )
389  $ GO TO 50
390  40 CONTINUE
391  50 CONTINUE
392  kr = i
393  END IF
394  END IF
395 *
396  IF( kl.NE.kln ) THEN
397  kln = kl
398 *
399 * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
400 * has not ben computed before.
401 *
402  hnorm = zlanhs( 'I', kr-kl+1, h( kl, kl ), ldh, rwork )
403  IF( disnan( hnorm ) ) THEN
404  info = -6
405  RETURN
406  ELSE IF( hnorm.GT.rzero ) THEN
407  eps3 = hnorm*ulp
408  ELSE
409  eps3 = smlnum
410  END IF
411  END IF
412 *
413 * Perturb eigenvalue if it is close to any previous
414 * selected eigenvalues affiliated to the submatrix
415 * H(KL:KR,KL:KR). Close roots are modified by EPS3.
416 *
417  wk = w( k )
418  60 CONTINUE
419  DO 70 i = k - 1, kl, -1
420  IF( SELECT( i ) .AND. cabs1( w( i )-wk ).LT.eps3 ) THEN
421  wk = wk + eps3
422  GO TO 60
423  END IF
424  70 CONTINUE
425  w( k ) = wk
426 *
427  IF( leftv ) THEN
428 *
429 * Compute left eigenvector.
430 *
431  CALL zlaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,
432  $ wk, vl( kl, ks ), work, ldwork, rwork, eps3,
433  $ smlnum, iinfo )
434  IF( iinfo.GT.0 ) THEN
435  info = info + 1
436  ifaill( ks ) = k
437  ELSE
438  ifaill( ks ) = 0
439  END IF
440  DO 80 i = 1, kl - 1
441  vl( i, ks ) = zero
442  80 CONTINUE
443  END IF
444  IF( rightv ) THEN
445 *
446 * Compute right eigenvector.
447 *
448  CALL zlaein( .true., noinit, kr, h, ldh, wk, vr( 1, ks ),
449  $ work, ldwork, rwork, eps3, smlnum, iinfo )
450  IF( iinfo.GT.0 ) THEN
451  info = info + 1
452  ifailr( ks ) = k
453  ELSE
454  ifailr( ks ) = 0
455  END IF
456  DO 90 i = kr + 1, n
457  vr( i, ks ) = zero
458  90 CONTINUE
459  END IF
460  ks = ks + 1
461  END IF
462  100 CONTINUE
463 *
464  RETURN
465 *
466 * End of ZHSEIN
467 *
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
subroutine zlaein(RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, EPS3, SMLNUM, INFO)
ZLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iterat...
Definition: zlaein.f:151
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function zlanhs(NORM, N, A, LDA, WORK)
ZLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlanhs.f:111
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: