strsna.f

*> \brief \b STRSNA
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download STRSNA + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strsna.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strsna.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strsna.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE STRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
*                          LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
*                          INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          HOWMNY, JOB
*       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
*       ..
*       .. Array Arguments ..
*       LOGICAL            SELECT( * )
*       INTEGER            IWORK( * )
*       REAL               S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
*      $                   VR( LDVR, * ), WORK( LDWORK, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> STRSNA estimates reciprocal condition numbers for specified
*> eigenvalues and/or right eigenvectors of a real upper
*> quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
*> orthogonal).
*>
*> T must be in Schur canonical form (as returned by SHSEQR), that is,
*> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
*> 2-by-2 diagonal block has its diagonal elements equal and its
*> off-diagonal elements of opposite sign.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOB
*> \verbatim
*>          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).
*> \endverbatim
*>
*> \param[in] HOWMNY
*> \verbatim
*>          HOWMNY is CHARACTER*1
*>          = 'A': compute condition numbers for all eigenpairs;
*>          = 'S': compute condition numbers for selected eigenpairs
*>                 specified by the array SELECT.
*> \endverbatim
*>
*> \param[in] SELECT
*> \verbatim
*>          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 eigenpair corresponding to a real eigenvalue w(j),
*>          SELECT(j) must be set to .TRUE.. To select condition numbers
*>          corresponding to a complex conjugate pair of eigenvalues w(j)
*>          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
*>          set to .TRUE..
*>          If HOWMNY = 'A', SELECT is not referenced.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix T. N >= 0.
*> \endverbatim
*>
*> \param[in] T
*> \verbatim
*>          T is REAL array, dimension (LDT,N)
*>          The upper quasi-triangular matrix T, in Schur canonical form.
*> \endverbatim
*>
*> \param[in] LDT
*> \verbatim
*>          LDT is INTEGER
*>          The leading dimension of the array T. LDT >= max(1,N).
*> \endverbatim
*>
*> \param[in] VL
*> \verbatim
*>          VL is REAL array, dimension (LDVL,M)
*>          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
*>          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
*>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
*>          must be stored in consecutive columns of VL, as returned by
*>          SHSEIN or STREVC.
*>          If JOB = 'V', VL is not referenced.
*> \endverbatim
*>
*> \param[in] LDVL
*> \verbatim
*>          LDVL is INTEGER
*>          The leading dimension of the array VL.
*>          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
*> \endverbatim
*>
*> \param[in] VR
*> \verbatim
*>          VR is REAL array, dimension (LDVR,M)
*>          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
*>          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
*>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
*>          must be stored in consecutive columns of VR, as returned by
*>          SHSEIN or STREVC.
*>          If JOB = 'V', VR is not referenced.
*> \endverbatim
*>
*> \param[in] LDVR
*> \verbatim
*>          LDVR is INTEGER
*>          The leading dimension of the array VR.
*>          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*>          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. For a complex conjugate pair of eigenvalues two
*>          consecutive elements of S are set to the same value. 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.
*> \endverbatim
*>
*> \param[out] SEP
*> \verbatim
*>          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. For a complex eigenvector two
*>          consecutive elements of SEP are set to the same value. If
*>          the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
*>          is set to 0; this can only occur when the true value would be
*>          very small anyway.
*>          If JOB = 'E', SEP is not referenced.
*> \endverbatim
*>
*> \param[in] MM
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[out] M
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (LDWORK,N+6)
*>          If JOB = 'E', WORK is not referenced.
*> \endverbatim
*>
*> \param[in] LDWORK
*> \verbatim
*>          LDWORK is INTEGER
*>          The leading dimension of the array WORK.
*>          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (2*(N-1))
*>          If JOB = 'E', IWORK is not referenced.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0: successful exit
*>          < 0: if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup trsna
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  The reciprocal of the condition number of an eigenvalue lambda is
*>  defined as
*>
*>          S(lambda) = |v**T*u| / (norm(u)*norm(v))
*>
*>  where u and v are the right and left eigenvectors of T corresponding
*>  to lambda; v**T denotes the 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)
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE strsna( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
     $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
     $                   INFO )
*
*  -- 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( * )
      INTEGER            IWORK( * )
      REAL               S( * ), SEP( * ), T( ldt, * ), VL( ldvl, * ),
     $                   vr( ldvr, * ), work( ldwork, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, TWO
      parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            PAIR, SOMCON, WANTBH, WANTS, WANTSP
      INTEGER            I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
      REAL               BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
     $                   mu, prod, prod1, prod2, rnrm, scale, smlnum, sn
*     ..
*     .. Local Arrays ..
      INTEGER            ISAVE( 3 )
      REAL               DUMMY( 1 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SDOT, SLAMCH, SLAPY2, SNRM2
      EXTERNAL           lsame, sdot, slamch, slapy2, snrm2
*     ..
*     .. External Subroutines ..
      EXTERNAL           slabad, slacn2, slacpy, slaqtr, strexc, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, sqrt
*     ..
*     .. 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' )
*
      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
*
*        Set M to the number of eigenpairs for which condition numbers
*        are required, and test MM.
*
         IF( somcon ) THEN
            m = 0
            pair = .false.
            DO 10 k = 1, n
               IF( pair ) THEN
                  pair = .false.
               ELSE
                  IF( k.LT.n ) THEN
                     IF( t( k+1, k ).EQ.zero ) THEN
                        IF( SELECT( k ) )
     $                     m = m + 1
                     ELSE
                        pair = .true.
                        IF( SELECT( k ) .OR. SELECT( k+1 ) )
     $                     m = m + 2
                     END IF
                  ELSE
                     IF( SELECT( n ) )
     $                  m = m + 1
                  END IF
               END IF
   10       CONTINUE
         ELSE
            m = n
         END IF
*
         IF( mm.LT.m ) THEN
            info = -13
         ELSE IF( ldwork.LT.1 .OR. ( wantsp .AND. ldwork.LT.n ) ) THEN
            info = -16
         END IF
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'STRSNA', -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
      CALL slabad( smlnum, bignum )
*
      ks = 0
      pair = .false.
      DO 60 k = 1, n
*
*        Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
*
         IF( pair ) THEN
            pair = .false.
            GO TO 60
         ELSE
            IF( k.LT.n )
     $         pair = t( k+1, k ).NE.zero
         END IF
*
*        Determine whether condition numbers are required for the k-th
*        eigenpair.
*
         IF( somcon ) THEN
            IF( pair ) THEN
               IF( .NOT.SELECT( k ) .AND. .NOT.SELECT( k+1 ) )
     $            GO TO 60
            ELSE
               IF( .NOT.SELECT( k ) )
     $            GO TO 60
            END IF
         END IF
*
         ks = ks + 1
*
         IF( wants ) THEN
*
*           Compute the reciprocal condition number of the k-th
*           eigenvalue.
*
            IF( .NOT.pair ) THEN
*
*              Real eigenvalue.
*
               prod = sdot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
               rnrm = snrm2( n, vr( 1, ks ), 1 )
               lnrm = snrm2( n, vl( 1, ks ), 1 )
               s( ks ) = abs( prod ) / ( rnrm*lnrm )
            ELSE
*
*              Complex eigenvalue.
*
               prod1 = sdot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
               prod1 = prod1 + sdot( n, vr( 1, ks+1 ), 1, vl( 1, ks+1 ),
     $                 1 )
               prod2 = sdot( n, vl( 1, ks ), 1, vr( 1, ks+1 ), 1 )
               prod2 = prod2 - sdot( n, vl( 1, ks+1 ), 1, vr( 1, ks ),
     $                 1 )
               rnrm = slapy2( snrm2( n, vr( 1, ks ), 1 ),
     $                snrm2( n, vr( 1, ks+1 ), 1 ) )
               lnrm = slapy2( snrm2( n, vl( 1, ks ), 1 ),
     $                snrm2( n, vl( 1, ks+1 ), 1 ) )
               cond = slapy2( prod1, prod2 ) / ( rnrm*lnrm )
               s( ks ) = cond
               s( ks+1 ) = cond
            END IF
         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 diagonal
*           block beginning at T(k,k) to the (1,1) position.
*
            CALL slacpy( 'Full', n, n, t, ldt, work, ldwork )
            ifst = k
            ilst = 1
            CALL strexc( 'No Q', n, work, ldwork, dummy, 1, ifst, ilst,
     $                   work( 1, n+1 ), ierr )
*
            IF( ierr.EQ.1 .OR. ierr.EQ.2 ) THEN
*
*              Could not swap because blocks not well separated
*
               scale = one
               est = bignum
            ELSE
*
*              Reordering successful
*
               IF( work( 2, 1 ).EQ.zero ) THEN
*
*                 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
                  n2 = 1
                  nn = n - 1
               ELSE
*
*                 Triangularize the 2 by 2 block by unitary
*                 transformation U = [  cs   i*ss ]
*                                    [ i*ss   cs  ].
*                 such that the (1,1) position of WORK is complex
*                 eigenvalue lambda with positive imaginary part. (2,2)
*                 position of WORK is the complex eigenvalue lambda
*                 with negative imaginary  part.
*
                  mu = sqrt( abs( work( 1, 2 ) ) )*
     $                 sqrt( abs( work( 2, 1 ) ) )
                  delta = slapy2( mu, work( 2, 1 ) )
                  cs = mu / delta
                  sn = -work( 2, 1 ) / delta
*
*                 Form
*
*                 C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
*                                          [   mu                     ]
*                                          [         ..               ]
*                                          [             ..           ]
*                                          [                  mu      ]
*                 where C**T is transpose of matrix C,
*                 and RWORK is stored starting in the N+1-st column of
*                 WORK.
*
                  DO 30 j = 3, n
                     work( 2, j ) = cs*work( 2, j )
                     work( j, j ) = work( j, j ) - work( 1, 1 )
   30             CONTINUE
                  work( 2, 2 ) = zero
*
                  work( 1, n+1 ) = two*mu
                  DO 40 i = 2, n - 1
                     work( i, n+1 ) = sn*work( 1, i+1 )
   40             CONTINUE
                  n2 = 2
                  nn = 2*( n-1 )
               END IF
*
*              Estimate norm(inv(C**T))
*
               est = zero
               kase = 0
   50          CONTINUE
               CALL slacn2( nn, work( 1, n+2 ), work( 1, n+4 ), iwork,
     $                      est, kase, isave )
               IF( kase.NE.0 ) THEN
                  IF( kase.EQ.1 ) THEN
                     IF( n2.EQ.1 ) THEN
*
*                       Real eigenvalue: solve C**T*x = scale*c.
*
                        CALL slaqtr( .true., .true., n-1, work( 2, 2 ),
     $                               ldwork, dummy, dumm, scale,
     $                               work( 1, n+4 ), work( 1, n+6 ),
     $                               ierr )
                     ELSE
*
*                       Complex eigenvalue: solve
*                       C**T*(p+iq) = scale*(c+id) in real arithmetic.
*
                        CALL slaqtr( .true., .false., n-1, work( 2, 2 ),
     $                               ldwork, work( 1, n+1 ), mu, scale,
     $                               work( 1, n+4 ), work( 1, n+6 ),
     $                               ierr )
                     END IF
                  ELSE
                     IF( n2.EQ.1 ) THEN
*
*                       Real eigenvalue: solve C*x = scale*c.
*
                        CALL slaqtr( .false., .true., n-1, work( 2, 2 ),
     $                               ldwork, dummy, dumm, scale,
     $                               work( 1, n+4 ), work( 1, n+6 ),
     $                               ierr )
                     ELSE
*
*                       Complex eigenvalue: solve
*                       C*(p+iq) = scale*(c+id) in real arithmetic.
*
                        CALL slaqtr( .false., .false., n-1,
     $                               work( 2, 2 ), ldwork,
     $                               work( 1, n+1 ), mu, scale,
     $                               work( 1, n+4 ), work( 1, n+6 ),
     $                               ierr )
*
                     END IF
                  END IF
*
                  GO TO 50
               END IF
            END IF
*
            sep( ks ) = scale / max( est, smlnum )
            IF( pair )
     $         sep( ks+1 ) = sep( ks )
         END IF
*
         IF( pair )
     $      ks = ks + 1
*
   60 CONTINUE
      RETURN
*
*     End of STRSNA
*
      END