LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ chsein()

subroutine chsein ( character  SIDE,
character  EIGSRC,
character  INITV,
logical, dimension( * )  SELECT,
integer  N,
complex, dimension( ldh, * )  H,
integer  LDH,
complex, dimension( * )  W,
complex, dimension( ldvl, * )  VL,
integer  LDVL,
complex, dimension( ldvr, * )  VR,
integer  LDVR,
integer  MM,
integer  M,
complex, dimension( * )  WORK,
real, dimension( * )  RWORK,
integer, dimension( * )  IFAILL,
integer, dimension( * )  IFAILR,
integer  INFO 
)

CHSEIN

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

Purpose:
 CHSEIN 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 CHSEQR; 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 CHSEIN 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, CHSEIN 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 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 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 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 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 array, dimension (N*N)
[out]RWORK
          RWORK is REAL 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
December 2016
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 chsein.f.

247 *
248 * -- LAPACK computational routine (version 3.7.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 * December 2016
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  REAL rwork( * )
261  COMPLEX h( ldh, * ), vl( ldvl, * ), vr( ldvr, * ),
262  $ w( * ), work( * )
263 * ..
264 *
265 * =====================================================================
266 *
267 * .. Parameters ..
268  COMPLEX zero
269  parameter( zero = ( 0.0e+0, 0.0e+0 ) )
270  REAL rzero
271  parameter( rzero = 0.0e+0 )
272 * ..
273 * .. Local Scalars ..
274  LOGICAL bothv, fromqr, leftv, noinit, rightv
275  INTEGER i, iinfo, k, kl, kln, kr, ks, ldwork
276  REAL eps3, hnorm, smlnum, ulp, unfl
277  COMPLEX cdum, wk
278 * ..
279 * .. External Functions ..
280  LOGICAL lsame, sisnan
281  REAL clanhs, slamch
282  EXTERNAL lsame, clanhs, slamch, sisnan
283 * ..
284 * .. External Subroutines ..
285  EXTERNAL claein, xerbla
286 * ..
287 * .. Intrinsic Functions ..
288  INTRINSIC abs, aimag, max, real
289 * ..
290 * .. Statement Functions ..
291  REAL cabs1
292 * ..
293 * .. Statement Function definitions ..
294  cabs1( cdum ) = abs( REAL( CDUM ) ) + abs( aimag( 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( 'CHSEIN', -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 = slamch( 'Safe minimum' )
348  ulp = slamch( '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 = clanhs( 'I', kr-kl+1, h( kl, kl ), ldh, rwork )
403  IF( sisnan( 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 claein( .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 claein( .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 CHSEIN
467 *
subroutine claein(RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, EPS3, SMLNUM, INFO)
CLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iterat...
Definition: claein.f:151
real function clanhs(NORM, N, A, LDA, WORK)
CLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clanhs.f:111
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function sisnan(SIN)
SISNAN tests input for NaN.
Definition: sisnan.f:61
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
Here is the call graph for this function:
Here is the caller graph for this function: