ztgevc.f

*> \brief \b ZTGEVC
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download ZTGEVC + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgevc.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgevc.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgevc.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
*                          LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          HOWMNY, SIDE
*       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
*       ..
*       .. Array Arguments ..
*       LOGICAL            SELECT( * )
*       DOUBLE PRECISION   RWORK( * )
*       COMPLEX*16         P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
*      $                   VR( LDVR, * ), WORK( * )
*       ..
*
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZTGEVC computes some or all of the right and/or left eigenvectors of
*> a pair of complex matrices (S,P), where S and P are upper triangular.
*> Matrix pairs of this type are produced by the generalized Schur
*> factorization of a complex matrix pair (A,B):
*>
*>    A = Q*S*Z**H,  B = Q*P*Z**H
*>
*> as computed by ZGGHRD + ZHGEQZ.
*>
*> The right eigenvector x and the left eigenvector y of (S,P)
*> corresponding to an eigenvalue w are defined by:
*>
*>    S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
*>
*> where y**H denotes the conjugate tranpose of y.
*> The eigenvalues are not input to this routine, but are computed
*> directly from the diagonal elements of S and P.
*>
*> This routine returns the matrices X and/or Y of right and left
*> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
*> where Z and Q are input matrices.
*> If Q and Z are the unitary factors from the generalized Schur
*> factorization of a matrix pair (A,B), then Z*X and Q*Y
*> are the matrices of right and left eigenvectors of (A,B).
*> \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 by the matrices in VR and/or VL;
*>          = 'S': compute selected right and/or left eigenvectors,
*>                 specified 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 matrices S and P.  N >= 0.
*> \endverbatim
*>
*> \param[in] S
*> \verbatim
*>          S is COMPLEX*16 array, dimension (LDS,N)
*>          The upper triangular matrix S from a generalized Schur
*>          factorization, as computed by ZHGEQZ.
*> \endverbatim
*>
*> \param[in] LDS
*> \verbatim
*>          LDS is INTEGER
*>          The leading dimension of array S.  LDS >= max(1,N).
*> \endverbatim
*>
*> \param[in] P
*> \verbatim
*>          P is COMPLEX*16 array, dimension (LDP,N)
*>          The upper triangular matrix P from a generalized Schur
*>          factorization, as computed by ZHGEQZ.  P must have real
*>          diagonal elements.
*> \endverbatim
*>
*> \param[in] LDP
*> \verbatim
*>          LDP is INTEGER
*>          The leading dimension of array P.  LDP >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] VL
*> \verbatim
*>          VL is COMPLEX*16 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 left Schur vectors returned by ZHGEQZ).
*>          On exit, if SIDE = 'L' or 'B', VL contains:
*>          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
*>          if HOWMNY = 'B', the matrix Q*Y;
*>          if HOWMNY = 'S', the left eigenvectors of (S,P) 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 array VL.  LDVL >= 1, and if
*>          SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
*> \endverbatim
*>
*> \param[in,out] VR
*> \verbatim
*>          VR is COMPLEX*16 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 Z
*>          of right Schur vectors returned by ZHGEQZ).
*>          On exit, if SIDE = 'R' or 'B', VR contains:
*>          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
*>          if HOWMNY = 'B', the matrix Z*X;
*>          if HOWMNY = 'S', the right eigenvectors of (S,P) 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*16 array, dimension (2*N)
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is DOUBLE PRECISION array, dimension (2*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.
*
*> \ingroup tgevc
*
*  =====================================================================
      SUBROUTINE ztgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
     $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, 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, SIDE
      INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
*     ..
*     .. Array Arguments ..
      LOGICAL            SELECT( * )
      DOUBLE PRECISION   RWORK( * )
      COMPLEX*16         P( ldp, * ), S( lds, * ), VL( ldvl, * ),
     $                   vr( ldvr, * ), work( * )
*     ..
*
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d+0, one = 1.0d+0 )
      COMPLEX*16         CZERO, CONE
      parameter( czero = ( 0.0d+0, 0.0d+0 ),
     $                   cone = ( 1.0d+0, 0.0d+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
     $                   lsa, lsb
      INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
     $                   j, je, jr
      DOUBLE PRECISION   ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
     $                   bignum, bnorm, bscale, dmin, safmin, sbeta,
     $                   scale, small, temp, ulp, xmax
      COMPLEX*16         BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH
      COMPLEX*16         ZLADIV
      EXTERNAL           lsame, dlamch, zladiv
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlabad, xerbla, zgemv
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, dcmplx, dconjg, dimag, max, min
*     ..
*     .. Statement Functions ..
      DOUBLE PRECISION   ABS1
*     ..
*     .. Statement Function definitions ..
      abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
*     ..
*     .. Executable Statements ..
*
*     Decode and Test the input parameters
*
      IF( lsame( howmny, 'A' ) ) THEN
         ihwmny = 1
         ilall = .true.
         ilback = .false.
      ELSE IF( lsame( howmny, 'S' ) ) THEN
         ihwmny = 2
         ilall = .false.
         ilback = .false.
      ELSE IF( lsame( howmny, 'B' ) ) THEN
         ihwmny = 3
         ilall = .true.
         ilback = .true.
      ELSE
         ihwmny = -1
      END IF
*
      IF( lsame( side, 'R' ) ) THEN
         iside = 1
         compl = .false.
         compr = .true.
      ELSE IF( lsame( side, 'L' ) ) THEN
         iside = 2
         compl = .true.
         compr = .false.
      ELSE IF( lsame( side, 'B' ) ) THEN
         iside = 3
         compl = .true.
         compr = .true.
      ELSE
         iside = -1
      END IF
*
      info = 0
      IF( iside.LT.0 ) THEN
         info = -1
      ELSE IF( ihwmny.LT.0 ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -4
      ELSE IF( lds.LT.max( 1, n ) ) THEN
         info = -6
      ELSE IF( ldp.LT.max( 1, n ) ) THEN
         info = -8
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZTGEVC', -info )
         RETURN
      END IF
*
*     Count the number of eigenvectors
*
      IF( .NOT.ilall ) THEN
         im = 0
         DO 10 j = 1, n
            IF( SELECT( j ) )
     $         im = im + 1
   10    CONTINUE
      ELSE
         im = n
      END IF
*
*     Check diagonal of B
*
      ilbbad = .false.
      DO 20 j = 1, n
         IF( dimag( p( j, j ) ).NE.zero )
     $      ilbbad = .true.
   20 CONTINUE
*
      IF( ilbbad ) THEN
         info = -7
      ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
         info = -10
      ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
         info = -12
      ELSE IF( mm.LT.im ) THEN
         info = -13
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZTGEVC', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      m = im
      IF( n.EQ.0 )
     $   RETURN
*
*     Machine Constants
*
      safmin = dlamch( 'Safe minimum' )
      big = one / safmin
      CALL dlabad( safmin, big )
      ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
      small = safmin*n / ulp
      big = one / small
      bignum = one / ( safmin*n )
*
*     Compute the 1-norm of each column of the strictly upper triangular
*     part of A and B to check for possible overflow in the triangular
*     solver.
*
      anorm = abs1( s( 1, 1 ) )
      bnorm = abs1( p( 1, 1 ) )
      rwork( 1 ) = zero
      rwork( n+1 ) = zero
      DO 40 j = 2, n
         rwork( j ) = zero
         rwork( n+j ) = zero
         DO 30 i = 1, j - 1
            rwork( j ) = rwork( j ) + abs1( s( i, j ) )
            rwork( n+j ) = rwork( n+j ) + abs1( p( i, j ) )
   30    CONTINUE
         anorm = max( anorm, rwork( j )+abs1( s( j, j ) ) )
         bnorm = max( bnorm, rwork( n+j )+abs1( p( j, j ) ) )
   40 CONTINUE
*
      ascale = one / max( anorm, safmin )
      bscale = one / max( bnorm, safmin )
*
*     Left eigenvectors
*
      IF( compl ) THEN
         ieig = 0
*
*        Main loop over eigenvalues
*
         DO 140 je = 1, n
            IF( ilall ) THEN
               ilcomp = .true.
            ELSE
               ilcomp = SELECT( je )
            END IF
            IF( ilcomp ) THEN
               ieig = ieig + 1
*
               IF( abs1( s( je, je ) ).LE.safmin .AND.
     $             abs( dble( p( je, je ) ) ).LE.safmin ) THEN
*
*                 Singular matrix pencil -- return unit eigenvector
*
                  DO 50 jr = 1, n
                     vl( jr, ieig ) = czero
   50             CONTINUE
                  vl( ieig, ieig ) = cone
                  GO TO 140
               END IF
*
*              Non-singular eigenvalue:
*              Compute coefficients  a  and  b  in
*                   H
*                 y  ( a A - b B ) = 0
*
               temp = one / max( abs1( s( je, je ) )*ascale,
     $                abs( dble( p( je, je ) ) )*bscale, safmin )
               salpha = ( temp*s( je, je ) )*ascale
               sbeta = ( temp*dble( p( je, je ) ) )*bscale
               acoeff = sbeta*ascale
               bcoeff = salpha*bscale
*
*              Scale to avoid underflow
*
               lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
               lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
     $               small
*
               scale = one
               IF( lsa )
     $            scale = ( small / abs( sbeta ) )*min( anorm, big )
               IF( lsb )
     $            scale = max( scale, ( small / abs1( salpha ) )*
     $                    min( bnorm, big ) )
               IF( lsa .OR. lsb ) THEN
                  scale = min( scale, one /
     $                    ( safmin*max( one, abs( acoeff ),
     $                    abs1( bcoeff ) ) ) )
                  IF( lsa ) THEN
                     acoeff = ascale*( scale*sbeta )
                  ELSE
                     acoeff = scale*acoeff
                  END IF
                  IF( lsb ) THEN
                     bcoeff = bscale*( scale*salpha )
                  ELSE
                     bcoeff = scale*bcoeff
                  END IF
               END IF
*
               acoefa = abs( acoeff )
               bcoefa = abs1( bcoeff )
               xmax = one
               DO 60 jr = 1, n
                  work( jr ) = czero
   60          CONTINUE
               work( je ) = cone
               dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
*
*                                              H
*              Triangular solve of  (a A - b B)  y = 0
*
*                                      H
*              (rowwise in  (a A - b B) , or columnwise in a A - b B)
*
               DO 100 j = je + 1, n
*
*                 Compute
*                       j-1
*                 SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
*                       k=je
*                 (Scale if necessary)
*
                  temp = one / xmax
                  IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GT.bignum*
     $                temp ) THEN
                     DO 70 jr = je, j - 1
                        work( jr ) = temp*work( jr )
   70                CONTINUE
                     xmax = one
                  END IF
                  suma = czero
                  sumb = czero
*
                  DO 80 jr = je, j - 1
                     suma = suma + dconjg( s( jr, j ) )*work( jr )
                     sumb = sumb + dconjg( p( jr, j ) )*work( jr )
   80             CONTINUE
                  sum = acoeff*suma - dconjg( bcoeff )*sumb
*
*                 Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
*
*                 with scaling and perturbation of the denominator
*
                  d = dconjg( acoeff*s( j, j )-bcoeff*p( j, j ) )
                  IF( abs1( d ).LE.dmin )
     $               d = dcmplx( dmin )
*
                  IF( abs1( d ).LT.one ) THEN
                     IF( abs1( sum ).GE.bignum*abs1( d ) ) THEN
                        temp = one / abs1( sum )
                        DO 90 jr = je, j - 1
                           work( jr ) = temp*work( jr )
   90                   CONTINUE
                        xmax = temp*xmax
                        sum = temp*sum
                     END IF
                  END IF
                  work( j ) = zladiv( -sum, d )
                  xmax = max( xmax, abs1( work( j ) ) )
  100          CONTINUE
*
*              Back transform eigenvector if HOWMNY='B'.
*
               IF( ilback ) THEN
                  CALL zgemv( 'N', n, n+1-je, cone, vl( 1, je ), ldvl,
     $                        work( je ), 1, czero, work( n+1 ), 1 )
                  isrc = 2
                  ibeg = 1
               ELSE
                  isrc = 1
                  ibeg = je
               END IF
*
*              Copy and scale eigenvector into column of VL
*
               xmax = zero
               DO 110 jr = ibeg, n
                  xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
  110          CONTINUE
*
               IF( xmax.GT.safmin ) THEN
                  temp = one / xmax
                  DO 120 jr = ibeg, n
                     vl( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
  120             CONTINUE
               ELSE
                  ibeg = n + 1
               END IF
*
               DO 130 jr = 1, ibeg - 1
                  vl( jr, ieig ) = czero
  130          CONTINUE
*
            END IF
  140    CONTINUE
      END IF
*
*     Right eigenvectors
*
      IF( compr ) THEN
         ieig = im + 1
*
*        Main loop over eigenvalues
*
         DO 250 je = n, 1, -1
            IF( ilall ) THEN
               ilcomp = .true.
            ELSE
               ilcomp = SELECT( je )
            END IF
            IF( ilcomp ) THEN
               ieig = ieig - 1
*
               IF( abs1( s( je, je ) ).LE.safmin .AND.
     $             abs( dble( p( je, je ) ) ).LE.safmin ) THEN
*
*                 Singular matrix pencil -- return unit eigenvector
*
                  DO 150 jr = 1, n
                     vr( jr, ieig ) = czero
  150             CONTINUE
                  vr( ieig, ieig ) = cone
                  GO TO 250
               END IF
*
*              Non-singular eigenvalue:
*              Compute coefficients  a  and  b  in
*
*              ( a A - b B ) x  = 0
*
               temp = one / max( abs1( s( je, je ) )*ascale,
     $                abs( dble( p( je, je ) ) )*bscale, safmin )
               salpha = ( temp*s( je, je ) )*ascale
               sbeta = ( temp*dble( p( je, je ) ) )*bscale
               acoeff = sbeta*ascale
               bcoeff = salpha*bscale
*
*              Scale to avoid underflow
*
               lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
               lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
     $               small
*
               scale = one
               IF( lsa )
     $            scale = ( small / abs( sbeta ) )*min( anorm, big )
               IF( lsb )
     $            scale = max( scale, ( small / abs1( salpha ) )*
     $                    min( bnorm, big ) )
               IF( lsa .OR. lsb ) THEN
                  scale = min( scale, one /
     $                    ( safmin*max( one, abs( acoeff ),
     $                    abs1( bcoeff ) ) ) )
                  IF( lsa ) THEN
                     acoeff = ascale*( scale*sbeta )
                  ELSE
                     acoeff = scale*acoeff
                  END IF
                  IF( lsb ) THEN
                     bcoeff = bscale*( scale*salpha )
                  ELSE
                     bcoeff = scale*bcoeff
                  END IF
               END IF
*
               acoefa = abs( acoeff )
               bcoefa = abs1( bcoeff )
               xmax = one
               DO 160 jr = 1, n
                  work( jr ) = czero
  160          CONTINUE
               work( je ) = cone
               dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
*
*              Triangular solve of  (a A - b B) x = 0  (columnwise)
*
*              WORK(1:j-1) contains sums w,
*              WORK(j+1:JE) contains x
*
               DO 170 jr = 1, je - 1
                  work( jr ) = acoeff*s( jr, je ) - bcoeff*p( jr, je )
  170          CONTINUE
               work( je ) = cone
*
               DO 210 j = je - 1, 1, -1
*
*                 Form x(j) := - w(j) / d
*                 with scaling and perturbation of the denominator
*
                  d = acoeff*s( j, j ) - bcoeff*p( j, j )
                  IF( abs1( d ).LE.dmin )
     $               d = dcmplx( dmin )
*
                  IF( abs1( d ).LT.one ) THEN
                     IF( abs1( work( j ) ).GE.bignum*abs1( d ) ) THEN
                        temp = one / abs1( work( j ) )
                        DO 180 jr = 1, je
                           work( jr ) = temp*work( jr )
  180                   CONTINUE
                     END IF
                  END IF
*
                  work( j ) = zladiv( -work( j ), d )
*
                  IF( j.GT.1 ) THEN
*
*                    w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
*
                     IF( abs1( work( j ) ).GT.one ) THEN
                        temp = one / abs1( work( j ) )
                        IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GE.
     $                      bignum*temp ) THEN
                           DO 190 jr = 1, je
                              work( jr ) = temp*work( jr )
  190                      CONTINUE
                        END IF
                     END IF
*
                     ca = acoeff*work( j )
                     cb = bcoeff*work( j )
                     DO 200 jr = 1, j - 1
                        work( jr ) = work( jr ) + ca*s( jr, j ) -
     $                               cb*p( jr, j )
  200                CONTINUE
                  END IF
  210          CONTINUE
*
*              Back transform eigenvector if HOWMNY='B'.
*
               IF( ilback ) THEN
                  CALL zgemv( 'N', n, je, cone, vr, ldvr, work, 1,
     $                        czero, work( n+1 ), 1 )
                  isrc = 2
                  iend = n
               ELSE
                  isrc = 1
                  iend = je
               END IF
*
*              Copy and scale eigenvector into column of VR
*
               xmax = zero
               DO 220 jr = 1, iend
                  xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
  220          CONTINUE
*
               IF( xmax.GT.safmin ) THEN
                  temp = one / xmax
                  DO 230 jr = 1, iend
                     vr( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
  230             CONTINUE
               ELSE
                  iend = 0
               END IF
*
               DO 240 jr = iend + 1, n
                  vr( jr, ieig ) = czero
  240          CONTINUE
*
            END IF
  250    CONTINUE
      END IF
*
      RETURN
*
*     End of ZTGEVC
*
      END