--- src/ctgevc.f	Thu Dec  5 20:18:07 2002
+++ debian/updates/src/ctgevc.f	Wed Dec  4 15:23:42 2002
@@ -1,19 +1,19 @@
-      SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
+      SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
      $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
 *
 *  -- LAPACK routine (version 3.0) --
 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 *     Courant Institute, Argonne National Lab, and Rice University
-*     June 30, 1999
+*     May 4, 2001
 *
 *     .. Scalar Arguments ..
       CHARACTER          HOWMNY, SIDE
-      INTEGER            INFO, LDA, LDB, LDVL, LDVR, M, MM, N
+      INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
 *     ..
 *     .. Array Arguments ..
       LOGICAL            SELECT( * )
       REAL               RWORK( * )
-      COMPLEX            A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
+      COMPLEX            P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
      $                   VR( LDVR, * ), WORK( * )
 *     ..
 *
@@ -21,28 +21,30 @@
 *  Purpose
 *  =======
 *
-*  CTGEVC computes some or all of the right and/or left generalized
-*  eigenvectors of a pair of complex upper triangular matrices (A,B).
-*
-*  The right generalized eigenvector x and the left generalized
-*  eigenvector y of (A,B) corresponding to a generalized eigenvalue
-*  w are defined by:
-*
-*          (A - wB) * x = 0  and  y**H * (A - wB) = 0
-*
+*  CTGEVC 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 CGGHRD + CHGEQZ.
+*  
+*  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.
-*
-*  If an eigenvalue w is determined by zero diagonal elements of both A
-*  and B, a unit vector is returned as the corresponding eigenvector.
-*
-*  If all eigenvectors are requested, the routine may either return
-*  the matrices X and/or Y of right or left eigenvectors of (A,B), or
-*  the products Z*X and/or Q*Y, where Z and Q are input unitary
-*  matrices.  If (A,B) was obtained from the generalized Schur
-*  factorization of an original pair of matrices
-*     (A0,B0) = (Q*A*Z**H,Q*B*Z**H),
-*  then Z*X and Q*Y are the matrices of right or left eigenvectors of
-*  A.
+*  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).
 *
 *  Arguments
 *  =========
@@ -54,66 +56,66 @@
 *
 *  HOWMNY  (input) CHARACTER*1
 *          = 'A': compute all right and/or left eigenvectors;
-*          = 'B': compute all right and/or left eigenvectors, and
-*                 backtransform them using the input matrices supplied
-*                 in VR and/or VL;
+*          = '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.
 *
 *  SELECT  (input) LOGICAL array, dimension (N)
 *          If HOWMNY='S', SELECT specifies the eigenvectors to be
-*          computed.
-*          If HOWMNY='A' or 'B', SELECT is not referenced.
-*          To select the eigenvector corresponding to the j-th
-*          eigenvalue, SELECT(j) must be set to .TRUE..
+*          computed.  The eigenvector corresponding to the j-th
+*          eigenvalue is computed if SELECT(j) = .TRUE..
+*          Not referenced if HOWMNY = 'A' or 'B'.
 *
 *  N       (input) INTEGER
-*          The order of the matrices A and B.  N >= 0.
-*
-*  A       (input) COMPLEX array, dimension (LDA,N)
-*          The upper triangular matrix A.
-*
-*  LDA     (input) INTEGER
-*          The leading dimension of array A.  LDA >= max(1,N).
+*          The order of the matrices S and P.  N >= 0.
 *
-*  B       (input) COMPLEX array, dimension (LDB,N)
-*          The upper triangular matrix B.  B must have real diagonal
-*          elements.
+*  S       (input) COMPLEX array, dimension (LDS,N)
+*          The upper triangular matrix S from a generalized Schur
+*          factorization, as computed by CHGEQZ.
+*
+*  LDS     (input) INTEGER
+*          The leading dimension of array S.  LDS >= max(1,N).
+*
+*  P       (input) COMPLEX array, dimension (LDP,N)
+*          The upper triangular matrix P from a generalized Schur
+*          factorization, as computed by CHGEQZ.  P must have real
+*          diagonal elements.
 *
-*  LDB     (input) INTEGER
-*          The leading dimension of array B.  LDB >= max(1,N).
+*  LDP     (input) INTEGER
+*          The leading dimension of array P.  LDP >= max(1,N).
 *
 *  VL      (input/output) 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 left Schur vectors returned by CHGEQZ).
 *          On exit, if SIDE = 'L' or 'B', VL contains:
-*          if HOWMNY = 'A', the matrix Y of left eigenvectors of (A,B);
+*          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 (A,B) specified by
+*          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.
-*          If SIDE = 'R', VL is not referenced.
+*          Not referenced if SIDE = 'R'.
 *
 *  LDVL    (input) INTEGER
-*          The leading dimension of array VL.
-*          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
+*          The leading dimension of array VL.  LDVL >= 1, and if
+*          SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
 *
 *  VR      (input/output) 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 Z
 *          of right Schur vectors returned by CHGEQZ).
 *          On exit, if SIDE = 'R' or 'B', VR contains:
-*          if HOWMNY = 'A', the matrix X of right eigenvectors of (A,B);
+*          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 (A,B) specified by
+*          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.
-*          If SIDE = 'L', VR is not referenced.
+*          Not referenced if SIDE = 'L'.
 *
 *  LDVR    (input) INTEGER
-*          The leading dimension of the array VR.
-*          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
+*          The leading dimension of the array VR.  LDVR >= 1, and if
+*          SIDE = 'R' or 'B', LDVR >= N.
 *
 *  MM      (input) INTEGER
 *          The number of columns in the arrays VL and/or VR. MM >= M.
@@ -180,7 +182,7 @@
          IHWMNY = 2
          ILALL = .FALSE.
          ILBACK = .FALSE.
-      ELSE IF( LSAME( HOWMNY, 'B' ) .OR. LSAME( HOWMNY, 'T' ) ) THEN
+      ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
          IHWMNY = 3
          ILALL = .TRUE.
          ILBACK = .TRUE.
@@ -211,9 +213,9 @@
          INFO = -2
       ELSE IF( N.LT.0 ) THEN
          INFO = -4
-      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+      ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
          INFO = -6
-      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+      ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
          INFO = -8
       END IF
       IF( INFO.NE.0 ) THEN
@@ -237,7 +239,7 @@
 *
       ILBBAD = .FALSE.
       DO 20 J = 1, N
-         IF( AIMAG( B( J, J ) ).NE.ZERO )
+         IF( AIMAG( P( J, J ) ).NE.ZERO )
      $      ILBBAD = .TRUE.
    20 CONTINUE
 *
@@ -275,19 +277,19 @@
 *     part of A and B to check for possible overflow in the triangular
 *     solver.
 *
-      ANORM = ABS1( A( 1, 1 ) )
-      BNORM = ABS1( B( 1, 1 ) )
+      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( A( I, J ) )
-            RWORK( N+J ) = RWORK( N+J ) + ABS1( B( I, J ) )
+            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( A( J, J ) ) )
-         BNORM = MAX( BNORM, RWORK( N+J )+ABS1( B( J, J ) ) )
+         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 )
@@ -309,8 +311,8 @@
             IF( ILCOMP ) THEN
                IEIG = IEIG + 1
 *
-               IF( ABS1( A( JE, JE ) ).LE.SAFMIN .AND.
-     $             ABS( REAL( B( JE, JE ) ) ).LE.SAFMIN ) THEN
+               IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
+     $             ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN
 *
 *                 Singular matrix pencil -- return unit eigenvector
 *
@@ -326,10 +328,10 @@
 *                   H
 *                 y  ( a A - b B ) = 0
 *
-               TEMP = ONE / MAX( ABS1( A( JE, JE ) )*ASCALE,
-     $                ABS( REAL( B( JE, JE ) ) )*BSCALE, SAFMIN )
-               SALPHA = ( TEMP*A( JE, JE ) )*ASCALE
-               SBETA = ( TEMP*REAL( B( JE, JE ) ) )*BSCALE
+               TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
+     $                ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN )
+               SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
+               SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE
                ACOEFF = SBETA*ASCALE
                BCOEFF = SALPHA*BSCALE
 *
@@ -380,7 +382,7 @@
 *
 *                 Compute
 *                       j-1
-*                 SUM = sum  conjg( a*A(k,j) - b*B(k,j) )*x(k)
+*                 SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
 *                       k=je
 *                 (Scale if necessary)
 *
@@ -396,16 +398,16 @@
                   SUMB = CZERO
 *
                   DO 80 JR = JE, J - 1
-                     SUMA = SUMA + CONJG( A( JR, J ) )*WORK( JR )
-                     SUMB = SUMB + CONJG( B( JR, J ) )*WORK( JR )
+                     SUMA = SUMA + CONJG( S( JR, J ) )*WORK( JR )
+                     SUMB = SUMB + CONJG( P( JR, J ) )*WORK( JR )
    80             CONTINUE
                   SUM = ACOEFF*SUMA - CONJG( BCOEFF )*SUMB
 *
-*                 Form x(j) = - SUM / conjg( a*A(j,j) - b*B(j,j) )
+*                 Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
 *
 *                 with scaling and perturbation of the denominator
 *
-                  D = CONJG( ACOEFF*A( J, J )-BCOEFF*B( J, J ) )
+                  D = CONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
                   IF( ABS1( D ).LE.DMIN )
      $               D = CMPLX( DMIN )
 *
@@ -475,8 +477,8 @@
             IF( ILCOMP ) THEN
                IEIG = IEIG - 1
 *
-               IF( ABS1( A( JE, JE ) ).LE.SAFMIN .AND.
-     $             ABS( REAL( B( JE, JE ) ) ).LE.SAFMIN ) THEN
+               IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
+     $             ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN
 *
 *                 Singular matrix pencil -- return unit eigenvector
 *
@@ -492,10 +494,10 @@
 *
 *              ( a A - b B ) x  = 0
 *
-               TEMP = ONE / MAX( ABS1( A( JE, JE ) )*ASCALE,
-     $                ABS( REAL( B( JE, JE ) ) )*BSCALE, SAFMIN )
-               SALPHA = ( TEMP*A( JE, JE ) )*ASCALE
-               SBETA = ( TEMP*REAL( B( JE, JE ) ) )*BSCALE
+               TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
+     $                ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN )
+               SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
+               SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE
                ACOEFF = SBETA*ASCALE
                BCOEFF = SALPHA*BSCALE
 *
@@ -542,7 +544,7 @@
 *              WORK(j+1:JE) contains x
 *
                DO 170 JR = 1, JE - 1
-                  WORK( JR ) = ACOEFF*A( JR, JE ) - BCOEFF*B( JR, JE )
+                  WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
   170          CONTINUE
                WORK( JE ) = CONE
 *
@@ -551,7 +553,7 @@
 *                 Form x(j) := - w(j) / d
 *                 with scaling and perturbation of the denominator
 *
-                  D = ACOEFF*A( J, J ) - BCOEFF*B( J, J )
+                  D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
                   IF( ABS1( D ).LE.DMIN )
      $               D = CMPLX( DMIN )
 *
@@ -568,7 +570,7 @@
 *
                   IF( J.GT.1 ) THEN
 *
-*                    w = w + x(j)*(a A(*,j) - b B(*,j) ) with scaling
+*                    w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
 *
                      IF( ABS1( WORK( J ) ).GT.ONE ) THEN
                         TEMP = ONE / ABS1( WORK( J ) )
@@ -583,8 +585,8 @@
                      CA = ACOEFF*WORK( J )
                      CB = BCOEFF*WORK( J )
                      DO 200 JR = 1, J - 1
-                        WORK( JR ) = WORK( JR ) + CA*A( JR, J ) -
-     $                               CB*B( JR, J )
+                        WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
+     $                               CB*P( JR, J )
   200                CONTINUE
                   END IF
   210          CONTINUE
