SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO ) * * -- LAPACK routine (instrumented to count operations, version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * * .. Scalar Arguments .. CHARACTER HOWMNY, SIDE INTEGER INFO, LDA, LDB, LDVL, LDVR, M, MM, N * .. * .. Array Arguments .. LOGICAL SELECT( * ) REAL RWORK( * ) COMPLEX A( LDA, * ), B( LDB, * ), VL( LDVL, * ), $ VR( LDVR, * ), WORK( * ) * .. * * ----------------------- Begin Timing Code ------------------------ * Common block to return operation count and iteration count * ITCNT is initialized to 0, OPS is only incremented * OPST is used to accumulate small contributions to OPS * to avoid roundoff error * .. Common blocks .. COMMON / LATIME / OPS, ITCNT * .. * .. Scalars in Common .. REAL ITCNT, OPS * .. * ------------------------ End Timing Code ------------------------- * * * 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 * * 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. * * Arguments * ========= * * SIDE (input) CHARACTER*1 * = 'R': compute right eigenvectors only; * = 'L': compute left eigenvectors only; * = 'B': compute both right and left eigenvectors. * * 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; * = '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.. * * 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). * * B (input) COMPLEX array, dimension (LDB,N) * The upper triangular matrix B. B must have real diagonal * elements. * * LDB (input) INTEGER * The leading dimension of array B. LDB >= 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 = 'B', the matrix Q*Y; * if HOWMNY = 'S', the left eigenvectors of (A,B) specified by * SELECT, stored consecutively in the columns of * VL, in the same order as their eigenvalues. * If SIDE = 'R', VL is not referenced. * * LDVL (input) INTEGER * The leading dimension of array VL. * LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. * * 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 = 'B', the matrix Z*X; * if HOWMNY = 'S', the right eigenvectors of (A,B) specified by * SELECT, stored consecutively in the columns of * VR, in the same order as their eigenvalues. * If SIDE = 'L', VR is not referenced. * * LDVR (input) INTEGER * The leading dimension of the array VR. * LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. * * MM (input) INTEGER * The leading dimension of the array VR. * LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. * * M (output) 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. * * WORK (workspace) COMPLEX array, dimension (2*N) * * RWORK (workspace) REAL array, dimension (2*N) * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * * ===================================================================== * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) COMPLEX CZERO, CONE PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), $ CONE = ( 1.0E+0, 0.0E+0 ) ) * .. * .. Local Scalars .. LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP, $ LSA, LSB INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, IOPST, ISIDE, $ ISRC, J, JE, JR REAL ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG, $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA, $ SCALE, SMALL, TEMP, ULP, XMAX COMPLEX BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X * .. * .. External Functions .. LOGICAL LSAME REAL SLAMCH COMPLEX CLADIV EXTERNAL CLADIV, LSAME, SLAMCH * .. * .. External Subroutines .. EXTERNAL CGEMV, SLABAD, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL * .. * .. Statement Functions .. REAL ABS1 * .. * .. Statement Function definitions .. ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( 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' ) .OR. LSAME( HOWMNY, 'T' ) ) 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( LDA.LT.MAX( 1, N ) ) THEN INFO = -6 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CTGEVC', -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( AIMAG( B( 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( 'CTGEVC', -INFO ) RETURN END IF * * Quick return if possible * M = IM IF( N.EQ.0 ) $ RETURN * * Machine Constants * SAFMIN = SLAMCH( 'Safe minimum' ) BIG = ONE / SAFMIN CALL SLABAD( SAFMIN, BIG ) ULP = SLAMCH( 'Epsilon' )*SLAMCH( '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( A( 1, 1 ) ) BNORM = ABS1( B( 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 ) ) 30 CONTINUE ANORM = MAX( ANORM, RWORK( J )+ABS1( A( J, J ) ) ) BNORM = MAX( BNORM, RWORK( N+J )+ABS1( B( J, J ) ) ) 40 CONTINUE * ASCALE = ONE / MAX( ANORM, SAFMIN ) BSCALE = ONE / MAX( BNORM, SAFMIN ) * ---------------------- Begin Timing Code ------------------------- OPS = OPS + REAL( 2*N**2+2*N+6 ) * ----------------------- End Timing Code -------------------------- * * 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( A( JE, JE ) ).LE.SAFMIN .AND. $ ABS( REAL( B( 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( A( JE, JE ) )*ASCALE, $ ABS( REAL( B( JE, JE ) ) )*BSCALE, SAFMIN ) SALPHA = ( TEMP*A( JE, JE ) )*ASCALE SBETA = ( TEMP*REAL( B( 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 * ----------------- Begin Timing Code ------------------ * Calculation of SALPHA through DMIN IOPST = 34 ELSE IOPST = 20 * ------------------ End Timing Code ------------------- 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*A(k,j) - b*B(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 * ---------------- Begin Timing Code ---------------- IOPST = IOPST + 2*( J-JE ) * ----------------- End Timing Code ----------------- END IF SUMA = CZERO SUMB = CZERO * DO 80 JR = JE, J - 1 SUMA = SUMA + CONJG( A( JR, J ) )*WORK( JR ) SUMB = SUMB + CONJG( B( 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) ) * * with scaling and perturbation of the denominator * D = CONJG( ACOEFF*A( J, J )-BCOEFF*B( J, J ) ) IF( ABS1( D ).LE.DMIN ) $ D = CMPLX( 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 * -------------- Begin Timing Code --------------- IOPST = IOPST + 2*( J-JE ) + 5 * --------------- End Timing Code ---------------- END IF END IF WORK( J ) = CLADIV( -SUM, D ) XMAX = MAX( XMAX, ABS1( WORK( J ) ) ) 100 CONTINUE * * Back transform eigenvector if HOWMNY='B'. * IF( ILBACK ) THEN CALL CGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL, $ WORK( JE ), 1, CZERO, WORK( N+1 ), 1 ) ISRC = 2 IBEG = 1 * ----------------- Begin Timing Code ------------------ IOPST = IOPST + 8*N*( N+1-JE ) * ------------------ End Timing Code ------------------- 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 * * ------------------- Begin Timing Code ------------------- OPS = OPS + ( REAL( 36*( N-JE )+8*( N-JE )*( N+1-JE )+3* $ ( N+1-IBEG )+1 )+REAL( IOPST ) ) * -------------------- End Timing Code -------------------- 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( A( JE, JE ) ).LE.SAFMIN .AND. $ ABS( REAL( B( 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( A( JE, JE ) )*ASCALE, $ ABS( REAL( B( JE, JE ) ) )*BSCALE, SAFMIN ) SALPHA = ( TEMP*A( JE, JE ) )*ASCALE SBETA = ( TEMP*REAL( B( 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 * ----------------- Begin Timing Code ------------------ * Calculation of SALPHA through DMIN IOPST = 34 ELSE IOPST = 20 * ------------------ End Timing Code ------------------- 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*A( JR, JE ) - BCOEFF*B( 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*A( J, J ) - BCOEFF*B( J, J ) IF( ABS1( D ).LE.DMIN ) $ D = CMPLX( 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 * -------------- Begin Timing Code --------------- IOPST = IOPST + 2*JE + 5 ELSE IOPST = IOPST + 3 * --------------- End Timing Code ---------------- END IF END IF * WORK( J ) = CLADIV( -WORK( J ), D ) * IF( J.GT.1 ) THEN * * w = w + x(j)*(a A(*,j) - b B(*,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 * ------------- Begin Timing Code ------------- IOPST = IOPST + 2*JE + 6 ELSE IOPST = IOPST + 6 * -------------- End Timing Code -------------- END IF END IF * 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 ) 200 CONTINUE END IF 210 CONTINUE * * Back transform eigenvector if HOWMNY='B'. * IF( ILBACK ) THEN CALL CGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1, $ CZERO, WORK( N+1 ), 1 ) ISRC = 2 IEND = N * ----------------- Begin Timing Code ------------------ IOPST = IOPST + 8*N*JE * ------------------ End Timing Code ------------------- 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 * * ------------------- Begin Timing Code ------------------- OPS = OPS + ( REAL( 30*( JE-2 )+8*( JE-1 )*( JE-2 )+3* $ IEND+22 )+REAL( IOPST ) ) * -------------------- End Timing Code -------------------- END IF 250 CONTINUE END IF * RETURN * * End of CTGEVC * END