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

◆ 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.
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 240 of file chsein.f.

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