ctrsna()

subroutine ctrsna	(	character	job,
		character	howmny,
		logical, dimension( * )	select,
		integer	n,
		complex, dimension( ldt, * )	t,
		integer	ldt,
		complex, dimension( ldvl, * )	vl,
		integer	ldvl,
		complex, dimension( ldvr, * )	vr,
		integer	ldvr,
		real, dimension( * )	s,
		real, dimension( * )	sep,
		integer	mm,
		integer	m,
		complex, dimension( ldwork, * )	work,
		integer	ldwork,
		real, dimension( * )	rwork,
		integer	info )

CTRSNA

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

Purpose:

!>
!> CTRSNA estimates reciprocal condition numbers for specified
!> eigenvalues and/or right eigenvectors of a complex upper triangular
!> matrix T (or of any matrix Q*T*Q**H with Q unitary).
!> 

Parameters

[in]	JOB	!> JOB is CHARACTER*1 !> Specifies whether condition numbers are required for !> eigenvalues (S) or eigenvectors (SEP): !> = 'E': for eigenvalues only (S); !> = 'V': for eigenvectors only (SEP); !> = 'B': for both eigenvalues and eigenvectors (S and SEP). !>
[in]	HOWMNY	!> HOWMNY is CHARACTER*1 !> = 'A': compute condition numbers for all eigenpairs; !> = 'S': compute condition numbers for selected eigenpairs !> specified by the array SELECT. !>
[in]	SELECT	!> SELECT is LOGICAL array, dimension (N) !> If HOWMNY = 'S', SELECT specifies the eigenpairs for which !> condition numbers are required. To select condition numbers !> for the j-th eigenpair, SELECT(j) must be set to .TRUE.. !> If HOWMNY = 'A', SELECT is not referenced. !>
[in]	N	!> N is INTEGER !> The order of the matrix T. N >= 0. !>
[in]	T	!> T is COMPLEX array, dimension (LDT,N) !> The upper triangular matrix T. !>
[in]	LDT	!> LDT is INTEGER !> The leading dimension of the array T. LDT >= max(1,N). !>
[in]	VL	!> VL is COMPLEX array, dimension (LDVL,M) !> If JOB = 'E' or 'B', VL must contain left eigenvectors of T !> (or of any Q*T*Q**H with Q unitary), corresponding to the !> eigenpairs specified by HOWMNY and SELECT. The eigenvectors !> must be stored in consecutive columns of VL, as returned by !> CHSEIN or CTREVC. !> If JOB = 'V', VL is not referenced. !>
[in]	LDVL	!> LDVL is INTEGER !> The leading dimension of the array VL. !> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N. !>
[in]	VR	!> VR is COMPLEX array, dimension (LDVR,M) !> If JOB = 'E' or 'B', VR must contain right eigenvectors of T !> (or of any Q*T*Q**H with Q unitary), corresponding to the !> eigenpairs specified by HOWMNY and SELECT. The eigenvectors !> must be stored in consecutive columns of VR, as returned by !> CHSEIN or CTREVC. !> If JOB = 'V', VR is not referenced. !>
[in]	LDVR	!> LDVR is INTEGER !> The leading dimension of the array VR. !> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N. !>
[out]	S	!> S is REAL array, dimension (MM) !> If JOB = 'E' or 'B', the reciprocal condition numbers of the !> selected eigenvalues, stored in consecutive elements of the !> array. Thus S(j), SEP(j), and the j-th columns of VL and VR !> all correspond to the same eigenpair (but not in general the !> j-th eigenpair, unless all eigenpairs are selected). !> If JOB = 'V', S is not referenced. !>
[out]	SEP	!> SEP is REAL array, dimension (MM) !> If JOB = 'V' or 'B', the estimated reciprocal condition !> numbers of the selected eigenvectors, stored in consecutive !> elements of the array. !> If JOB = 'E', SEP is not referenced. !>
[in]	MM	!> MM is INTEGER !> The number of elements in the arrays S (if JOB = 'E' or 'B') !> and/or SEP (if JOB = 'V' or 'B'). MM >= M. !>
[out]	M	!> M is INTEGER !> The number of elements of the arrays S and/or SEP actually !> used to store the estimated condition numbers. !> If HOWMNY = 'A', M is set to N. !>
[out]	WORK	!> WORK is COMPLEX array, dimension (LDWORK,N+6) !> If JOB = 'E', WORK is not referenced. !>
[in]	LDWORK	!> LDWORK is INTEGER !> The leading dimension of the array WORK. !> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N. !>
[out]	RWORK	!> RWORK is REAL array, dimension (N) !> If JOB = 'E', RWORK is not referenced. !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -i, the i-th argument had an illegal value !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  The reciprocal of the condition number of an eigenvalue lambda is
!>  defined as
!>
!>          S(lambda) = |v**H*u| / (norm(u)*norm(v))
!>
!>  where u and v are the right and left eigenvectors of T corresponding
!>  to lambda; v**H denotes the conjugate transpose of v, and norm(u)
!>  denotes the Euclidean norm. These reciprocal condition numbers always
!>  lie between zero (very badly conditioned) and one (very well
!>  conditioned). If n = 1, S(lambda) is defined to be 1.
!>
!>  An approximate error bound for a computed eigenvalue W(i) is given by
!>
!>                      EPS * norm(T) / S(i)
!>
!>  where EPS is the machine precision.
!>
!>  The reciprocal of the condition number of the right eigenvector u
!>  corresponding to lambda is defined as follows. Suppose
!>
!>              T = ( lambda  c  )
!>                  (   0    T22 )
!>
!>  Then the reciprocal condition number is
!>
!>          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
!>
!>  where sigma-min denotes the smallest singular value. We approximate
!>  the smallest singular value by the reciprocal of an estimate of the
!>  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
!>  defined to be abs(T(1,1)).
!>
!>  An approximate error bound for a computed right eigenvector VR(i)
!>  is given by
!>
!>                      EPS * norm(T) / SEP(i)
!> 

Definition at line 244 of file ctrsna.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          HOWMNY, JOB
      INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
*     ..
*     .. Array Arguments ..
      LOGICAL            SELECT( * )
      REAL               RWORK( * ), S( * ), SEP( * )
      COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
     $                   WORK( LDWORK, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            SOMCON, WANTBH, WANTS, WANTSP
      CHARACTER          NORMIN
      INTEGER            I, IERR, IX, J, K, KASE, KS
      REAL               BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
     $                   XNORM
      COMPLEX            CDUM, PROD
*     ..
*     .. Local Arrays ..
      INTEGER            ISAVE( 3 )
      COMPLEX            DUMMY( 1 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICAMAX
      REAL               SCNRM2, SLAMCH
      COMPLEX            CDOTC
      EXTERNAL           lsame, icamax, scnrm2, slamch,
     $                   cdotc
*     ..
*     .. External Subroutines ..
      EXTERNAL           clacn2, clacpy, clatrs, csrscl, ctrexc,
     $                   xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, aimag, max, real
*     ..
*     .. Statement Functions ..
      REAL               CABS1
*     ..
*     .. Statement Function definitions ..
      cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
*     ..
*     .. Executable Statements ..
*
*     Decode and test the input parameters
*
      wantbh = lsame( job, 'B' )
      wants = lsame( job, 'E' ) .OR. wantbh
      wantsp = lsame( job, 'V' ) .OR. wantbh
*
      somcon = lsame( howmny, 'S' )
*
*     Set M to the number of eigenpairs for which condition numbers are
*     to be computed.
*
      IF( somcon ) THEN
         m = 0
         DO 10 j = 1, n
            IF( SELECT( j ) )
     $         m = m + 1
   10    CONTINUE
      ELSE
         m = n
      END IF
*
      info = 0
      IF( .NOT.wants .AND. .NOT.wantsp ) THEN
         info = -1
      ELSE IF( .NOT.lsame( howmny, 'A' ) .AND. .NOT.somcon ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -4
      ELSE IF( ldt.LT.max( 1, n ) ) THEN
         info = -6
      ELSE IF( ldvl.LT.1 .OR. ( wants .AND. ldvl.LT.n ) ) THEN
         info = -8
      ELSE IF( ldvr.LT.1 .OR. ( wants .AND. ldvr.LT.n ) ) THEN
         info = -10
      ELSE IF( mm.LT.m ) THEN
         info = -13
      ELSE IF( ldwork.LT.1 .OR. ( wantsp .AND. ldwork.LT.n ) ) THEN
         info = -16
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CTRSNA', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 )
     $   RETURN
*
      IF( n.EQ.1 ) THEN
         IF( somcon ) THEN
            IF( .NOT.SELECT( 1 ) )
     $         RETURN
         END IF
         IF( wants )
     $      s( 1 ) = one
         IF( wantsp )
     $      sep( 1 ) = abs( t( 1, 1 ) )
         RETURN
      END IF
*
*     Get machine constants
*
      eps = slamch( 'P' )
      smlnum = slamch( 'S' ) / eps
      bignum = one / smlnum
*
      ks = 1
      DO 50 k = 1, n
*
         IF( somcon ) THEN
            IF( .NOT.SELECT( k ) )
     $         GO TO 50
         END IF
*
         IF( wants ) THEN
*
*           Compute the reciprocal condition number of the k-th
*           eigenvalue.
*
            prod = cdotc( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
            rnrm = scnrm2( n, vr( 1, ks ), 1 )
            lnrm = scnrm2( n, vl( 1, ks ), 1 )
            s( ks ) = abs( prod ) / ( rnrm*lnrm )
*
         END IF
*
         IF( wantsp ) THEN
*
*           Estimate the reciprocal condition number of the k-th
*           eigenvector.
*
*           Copy the matrix T to the array WORK and swap the k-th
*           diagonal element to the (1,1) position.
*
            CALL clacpy( 'Full', n, n, t, ldt, work, ldwork )
            CALL ctrexc( 'No Q', n, work, ldwork, dummy, 1, k, 1,
     $                   ierr )
*
*           Form  C = T22 - lambda*I in WORK(2:N,2:N).
*
            DO 20 i = 2, n
               work( i, i ) = work( i, i ) - work( 1, 1 )
   20       CONTINUE
*
*           Estimate a lower bound for the 1-norm of inv(C**H). The 1st
*           and (N+1)th columns of WORK are used to store work vectors.
*
            sep( ks ) = zero
            est = zero
            kase = 0
            normin = 'N'
   30       CONTINUE
            CALL clacn2( n-1, work( 1, n+1 ), work, est, kase,
     $                   isave )
*
            IF( kase.NE.0 ) THEN
               IF( kase.EQ.1 ) THEN
*
*                 Solve C**H*x = scale*b
*
                  CALL clatrs( 'Upper', 'Conjugate transpose',
     $                         'Nonunit', normin, n-1, work( 2, 2 ),
     $                         ldwork, work, scale, rwork, ierr )
               ELSE
*
*                 Solve C*x = scale*b
*
                  CALL clatrs( 'Upper', 'No transpose', 'Nonunit',
     $                         normin, n-1, work( 2, 2 ), ldwork, work,
     $                         scale, rwork, ierr )
               END IF
               normin = 'Y'
               IF( scale.NE.one ) THEN
*
*                 Multiply by 1/SCALE if doing so will not cause
*                 overflow.
*
                  ix = icamax( n-1, work, 1 )
                  xnorm = cabs1( work( ix, 1 ) )
                  IF( scale.LT.xnorm*smlnum .OR. scale.EQ.zero )
     $               GO TO 40
                  CALL csrscl( n, scale, work, 1 )
               END IF
               GO TO 30
            END IF
*
            sep( ks ) = one / max( est, smlnum )
         END IF
*
   40    CONTINUE
         ks = ks + 1
   50 CONTINUE
      RETURN
*
*     End of CTRSNA
*

Here is the call graph for this function:

Here is the caller graph for this function: