ctrevc.f

Go to the documentation of this file.

*> \brief \b CTREVC
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download CTREVC + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrevc.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrevc.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrevc.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
*                          LDVR, MM, M, WORK, RWORK, INFO )
* 
*       .. Scalar Arguments ..
*       CHARACTER          HOWMNY, SIDE
*       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
*       ..
*       .. Array Arguments ..
*       LOGICAL            SELECT( * )
*       REAL               RWORK( * )
*       COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
*      $                   WORK( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CTREVC computes some or all of the right and/or left eigenvectors of
*> a complex upper triangular matrix T.
*> Matrices of this type are produced by the Schur factorization of
*> a complex general matrix:  A = Q*T*Q**H, as computed by CHSEQR.
*> 
*> The right eigenvector x and the left eigenvector y of T corresponding
*> to an eigenvalue w are defined by:
*> 
*>              T*x = w*x,     (y**H)*T = w*(y**H)
*> 
*> where y**H denotes the conjugate transpose of the vector y.
*> The eigenvalues are not input to this routine, but are read directly
*> from the diagonal of T.
*> 
*> This routine returns the matrices X and/or Y of right and left
*> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
*> input matrix.  If Q is the unitary factor that reduces a matrix A to
*> Schur form T, then Q*X and Q*Y are the matrices of right and left
*> eigenvectors of A.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] SIDE
*> \verbatim
*>          SIDE is CHARACTER*1
*>          = 'R':  compute right eigenvectors only;
*>          = 'L':  compute left eigenvectors only;
*>          = 'B':  compute both right and left eigenvectors.
*> \endverbatim
*>
*> \param[in] HOWMNY
*> \verbatim
*>          HOWMNY is CHARACTER*1
*>          = 'A':  compute all right and/or left eigenvectors;
*>          = 'B':  compute all right and/or left eigenvectors,
*>                  backtransformed using the matrices supplied in
*>                  VR and/or VL;
*>          = 'S':  compute selected right and/or left eigenvectors,
*>                  as indicated by the logical array SELECT.
*> \endverbatim
*>
*> \param[in] SELECT
*> \verbatim
*>          SELECT is LOGICAL array, dimension (N)
*>          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
*>          computed.
*>          The eigenvector corresponding to the j-th eigenvalue is
*>          computed if SELECT(j) = .TRUE..
*>          Not referenced if HOWMNY = 'A' or 'B'.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix T. N >= 0.
*> \endverbatim
*>
*> \param[in,out] T
*> \verbatim
*>          T is COMPLEX array, dimension (LDT,N)
*>          The upper triangular matrix T.  T is modified, but restored
*>          on exit.
*> \endverbatim
*>
*> \param[in] LDT
*> \verbatim
*>          LDT is INTEGER
*>          The leading dimension of the array T. LDT >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] VL
*> \verbatim
*>          VL is COMPLEX array, dimension (LDVL,MM)
*>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
*>          contain an N-by-N matrix Q (usually the unitary matrix Q of
*>          Schur vectors returned by CHSEQR).
*>          On exit, if SIDE = 'L' or 'B', VL contains:
*>          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
*>          if HOWMNY = 'B', the matrix Q*Y;
*>          if HOWMNY = 'S', the left eigenvectors of T specified by
*>                           SELECT, stored consecutively in the columns
*>                           of VL, in the same order as their
*>                           eigenvalues.
*>          Not referenced if SIDE = 'R'.
*> \endverbatim
*>
*> \param[in] LDVL
*> \verbatim
*>          LDVL is INTEGER
*>          The leading dimension of the array VL.  LDVL >= 1, and if
*>          SIDE = 'L' or 'B', LDVL >= N.
*> \endverbatim
*>
*> \param[in,out] VR
*> \verbatim
*>          VR is COMPLEX array, dimension (LDVR,MM)
*>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
*>          contain an N-by-N matrix Q (usually the unitary matrix Q of
*>          Schur vectors returned by CHSEQR).
*>          On exit, if SIDE = 'R' or 'B', VR contains:
*>          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
*>          if HOWMNY = 'B', the matrix Q*X;
*>          if HOWMNY = 'S', the right eigenvectors of T specified by
*>                           SELECT, stored consecutively in the columns
*>                           of VR, in the same order as their
*>                           eigenvalues.
*>          Not referenced if SIDE = 'L'.
*> \endverbatim
*>
*> \param[in] LDVR
*> \verbatim
*>          LDVR is INTEGER
*>          The leading dimension of the array VR.  LDVR >= 1, and if
*>          SIDE = 'R' or 'B'; LDVR >= N.
*> \endverbatim
*>
*> \param[in] MM
*> \verbatim
*>          MM is INTEGER
*>          The number of columns in the arrays VL and/or VR. MM >= M.
*> \endverbatim
*>
*> \param[out] M
*> \verbatim
*>          M is INTEGER
*>          The number of columns in the arrays VL and/or VR actually
*>          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
*>          is set to N.  Each selected eigenvector occupies one
*>          column.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX array, dimension (2*N)
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is REAL array, dimension (N)
*> \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. 
*
*> \date November 2011
*
*> \ingroup complexOTHERcomputational
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  The algorithm used in this program is basically backward (forward)
*>  substitution, with scaling to make the the code robust against
*>  possible overflow.
*>
*>  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|.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE ctrevc( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
     $                   ldvr, mm, m, work, rwork, info )
*
*  -- LAPACK computational routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      CHARACTER          HOWMNY, SIDE
      INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
*     ..
*     .. Array Arguments ..
      LOGICAL            SELECT( * )
      REAL               RWORK( * )
      COMPLEX            T( ldt, * ), VL( ldvl, * ), VR( ldvr, * ),
     $                   work( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter                ( zero = 0.0e+0, one = 1.0e+0 )
      COMPLEX            CMZERO, CMONE
      parameter                ( cmzero = ( 0.0e+0, 0.0e+0 ),
     $                   cmone = ( 1.0e+0, 0.0e+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
      INTEGER            I, II, IS, J, K, KI
      REAL               OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
      COMPLEX            CDUM
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICAMAX
      REAL               SCASUM, SLAMCH
      EXTERNAL           lsame, icamax, scasum, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           ccopy, cgemv, clatrs, csscal, slabad, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, aimag, cmplx, conjg, 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
*
      bothv = lsame( side, 'B' )
      rightv = lsame( side, 'R' ) .OR. bothv
      leftv = lsame( side, 'L' ) .OR. bothv
*
      allv = lsame( howmny, 'A' )
      over = lsame( howmny, 'B' )
      somev = lsame( howmny, 'S' )
*
*     Set M to the number of columns required to store the selected
*     eigenvectors.
*
      IF( somev ) 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.rightv .AND. .NOT.leftv ) THEN
         info = -1
      ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev ) 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. ( leftv .AND. ldvl.LT.n ) ) THEN
         info = -8
      ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
         info = -10
      ELSE IF( mm.LT.m ) THEN
         info = -11
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CTREVC', -info )
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( n.EQ.0 )
     $   RETURN
*
*     Set the constants to control overflow.
*
      unfl = slamch( 'Safe minimum' )
      ovfl = one / unfl
      CALL slabad( unfl, ovfl )
      ulp = slamch( 'Precision' )
      smlnum = unfl*( n / ulp )
*
*     Store the diagonal elements of T in working array WORK.
*
      DO 20 i = 1, n
         work( i+n ) = t( i, i )
   20 CONTINUE
*
*     Compute 1-norm of each column of strictly upper triangular
*     part of T to control overflow in triangular solver.
*
      rwork( 1 ) = zero
      DO 30 j = 2, n
         rwork( j ) = scasum( j-1, t( 1, j ), 1 )
   30 CONTINUE
*
      IF( rightv ) THEN
*
*        Compute right eigenvectors.
*
         is = m
         DO 80 ki = n, 1, -1
*
            IF( somev ) THEN
               IF( .NOT.SELECT( ki ) )
     $            GO TO 80
            END IF
            smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
*
            work( 1 ) = cmone
*
*           Form right-hand side.
*
            DO 40 k = 1, ki - 1
               work( k ) = -t( k, ki )
   40       CONTINUE
*
*           Solve the triangular system:
*              (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
*
            DO 50 k = 1, ki - 1
               t( k, k ) = t( k, k ) - t( ki, ki )
               IF( cabs1( t( k, k ) ).LT.smin )
     $            t( k, k ) = smin
   50       CONTINUE
*
            IF( ki.GT.1 ) THEN
               CALL clatrs( 'Upper', 'No transpose', 'Non-unit', 'Y',
     $                      ki-1, t, ldt, work( 1 ), scale, rwork,
     $                      info )
               work( ki ) = scale
            END IF
*
*           Copy the vector x or Q*x to VR and normalize.
*
            IF( .NOT.over ) THEN
               CALL ccopy( ki, work( 1 ), 1, vr( 1, is ), 1 )
*
               ii = icamax( ki, vr( 1, is ), 1 )
               remax = one / cabs1( vr( ii, is ) )
               CALL csscal( ki, remax, vr( 1, is ), 1 )
*
               DO 60 k = ki + 1, n
                  vr( k, is ) = cmzero
   60          CONTINUE
            ELSE
               IF( ki.GT.1 )
     $            CALL cgemv( 'N', n, ki-1, cmone, vr, ldvr, work( 1 ),
     $                        1, cmplx( scale ), vr( 1, ki ), 1 )
*
               ii = icamax( n, vr( 1, ki ), 1 )
               remax = one / cabs1( vr( ii, ki ) )
               CALL csscal( n, remax, vr( 1, ki ), 1 )
            END IF
*
*           Set back the original diagonal elements of T.
*
            DO 70 k = 1, ki - 1
               t( k, k ) = work( k+n )
   70       CONTINUE
*
            is = is - 1
   80    CONTINUE
      END IF
*
      IF( leftv ) THEN
*
*        Compute left eigenvectors.
*
         is = 1
         DO 130 ki = 1, n
*
            IF( somev ) THEN
               IF( .NOT.SELECT( ki ) )
     $            GO TO 130
            END IF
            smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
*
            work( n ) = cmone
*
*           Form right-hand side.
*
            DO 90 k = ki + 1, n
               work( k ) = -conjg( t( ki, k ) )
   90       CONTINUE
*
*           Solve the triangular system:
*              (T(KI+1:N,KI+1:N) - T(KI,KI))**H*X = SCALE*WORK.
*
            DO 100 k = ki + 1, n
               t( k, k ) = t( k, k ) - t( ki, ki )
               IF( cabs1( t( k, k ) ).LT.smin )
     $            t( k, k ) = smin
  100       CONTINUE
*
            IF( ki.LT.n ) THEN
               CALL clatrs( 'Upper', 'Conjugate transpose', 'Non-unit',
     $                      'Y', n-ki, t( ki+1, ki+1 ), ldt,
     $                      work( ki+1 ), scale, rwork, info )
               work( ki ) = scale
            END IF
*
*           Copy the vector x or Q*x to VL and normalize.
*
            IF( .NOT.over ) THEN
               CALL ccopy( n-ki+1, work( ki ), 1, vl( ki, is ), 1 )
*
               ii = icamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
               remax = one / cabs1( vl( ii, is ) )
               CALL csscal( n-ki+1, remax, vl( ki, is ), 1 )
*
               DO 110 k = 1, ki - 1
                  vl( k, is ) = cmzero
  110          CONTINUE
            ELSE
               IF( ki.LT.n )
     $            CALL cgemv( 'N', n, n-ki, cmone, vl( 1, ki+1 ), ldvl,
     $                        work( ki+1 ), 1, cmplx( scale ),
     $                        vl( 1, ki ), 1 )
*
               ii = icamax( n, vl( 1, ki ), 1 )
               remax = one / cabs1( vl( ii, ki ) )
               CALL csscal( n, remax, vl( 1, ki ), 1 )
            END IF
*
*           Set back the original diagonal elements of T.
*
            DO 120 k = ki + 1, n
               t( k, k ) = work( k+n )
  120       CONTINUE
*
            is = is + 1
  130    CONTINUE
      END IF
*
      RETURN
*
*     End of CTREVC
*
      END