3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download ZTGEVC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgevc.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgevc.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgevc.f">
21 * SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
22 * LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
24 * .. Scalar Arguments ..
25 * CHARACTER HOWMNY, SIDE
26 * INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
28 * .. Array Arguments ..
30 * DOUBLE PRECISION RWORK( * )
31 * COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
32 * $ VR( LDVR, * ), WORK( * )
42 *> ZTGEVC computes some or all of the right and/or left eigenvectors of
43 *> a pair of complex matrices (S,P), where S and P are upper triangular.
44 *> Matrix pairs of this type are produced by the generalized Schur
45 *> factorization of a complex matrix pair (A,B):
47 *> A = Q*S*Z**H, B = Q*P*Z**H
49 *> as computed by ZGGHRD + ZHGEQZ.
51 *> The right eigenvector x and the left eigenvector y of (S,P)
52 *> corresponding to an eigenvalue w are defined by:
54 *> S*x = w*P*x, (y**H)*S = w*(y**H)*P,
56 *> where y**H denotes the conjugate tranpose of y.
57 *> The eigenvalues are not input to this routine, but are computed
58 *> directly from the diagonal elements of S and P.
60 *> This routine returns the matrices X and/or Y of right and left
61 *> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
62 *> where Z and Q are input matrices.
63 *> If Q and Z are the unitary factors from the generalized Schur
64 *> factorization of a matrix pair (A,B), then Z*X and Q*Y
65 *> are the matrices of right and left eigenvectors of (A,B).
73 *> SIDE is CHARACTER*1
74 *> = 'R': compute right eigenvectors only;
75 *> = 'L': compute left eigenvectors only;
76 *> = 'B': compute both right and left eigenvectors.
81 *> HOWMNY is CHARACTER*1
82 *> = 'A': compute all right and/or left eigenvectors;
83 *> = 'B': compute all right and/or left eigenvectors,
84 *> backtransformed by the matrices in VR and/or VL;
85 *> = 'S': compute selected right and/or left eigenvectors,
86 *> specified by the logical array SELECT.
91 *> SELECT is LOGICAL array, dimension (N)
92 *> If HOWMNY='S', SELECT specifies the eigenvectors to be
93 *> computed. The eigenvector corresponding to the j-th
94 *> eigenvalue is computed if SELECT(j) = .TRUE..
95 *> Not referenced if HOWMNY = 'A' or 'B'.
101 *> The order of the matrices S and P. N >= 0.
106 *> S is COMPLEX*16 array, dimension (LDS,N)
107 *> The upper triangular matrix S from a generalized Schur
108 *> factorization, as computed by ZHGEQZ.
114 *> The leading dimension of array S. LDS >= max(1,N).
119 *> P is COMPLEX*16 array, dimension (LDP,N)
120 *> The upper triangular matrix P from a generalized Schur
121 *> factorization, as computed by ZHGEQZ. P must have real
122 *> diagonal elements.
128 *> The leading dimension of array P. LDP >= max(1,N).
133 *> VL is COMPLEX*16 array, dimension (LDVL,MM)
134 *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
135 *> contain an N-by-N matrix Q (usually the unitary matrix Q
136 *> of left Schur vectors returned by ZHGEQZ).
137 *> On exit, if SIDE = 'L' or 'B', VL contains:
138 *> if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
139 *> if HOWMNY = 'B', the matrix Q*Y;
140 *> if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
141 *> SELECT, stored consecutively in the columns of
142 *> VL, in the same order as their eigenvalues.
143 *> Not referenced if SIDE = 'R'.
149 *> The leading dimension of array VL. LDVL >= 1, and if
150 *> SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
155 *> VR is COMPLEX*16 array, dimension (LDVR,MM)
156 *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
157 *> contain an N-by-N matrix Q (usually the unitary matrix Z
158 *> of right Schur vectors returned by ZHGEQZ).
159 *> On exit, if SIDE = 'R' or 'B', VR contains:
160 *> if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
161 *> if HOWMNY = 'B', the matrix Z*X;
162 *> if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
163 *> SELECT, stored consecutively in the columns of
164 *> VR, in the same order as their eigenvalues.
165 *> Not referenced if SIDE = 'L'.
171 *> The leading dimension of the array VR. LDVR >= 1, and if
172 *> SIDE = 'R' or 'B', LDVR >= N.
178 *> The number of columns in the arrays VL and/or VR. MM >= M.
184 *> The number of columns in the arrays VL and/or VR actually
185 *> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
186 *> is set to N. Each selected eigenvector occupies one column.
191 *> WORK is COMPLEX*16 array, dimension (2*N)
196 *> RWORK is DOUBLE PRECISION array, dimension (2*N)
202 *> = 0: successful exit.
203 *> < 0: if INFO = -i, the i-th argument had an illegal value.
209 *> \author Univ. of Tennessee
210 *> \author Univ. of California Berkeley
211 *> \author Univ. of Colorado Denver
214 *> \date November 2011
216 *> \ingroup complex16GEcomputational
218 * =====================================================================
219 SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
220 $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
222 * -- LAPACK computational routine (version 3.4.0) --
223 * -- LAPACK is a software package provided by Univ. of Tennessee, --
224 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
227 * .. Scalar Arguments ..
228 CHARACTER HOWMNY, SIDE
229 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
231 * .. Array Arguments ..
233 DOUBLE PRECISION RWORK( * )
234 COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
235 $ VR( LDVR, * ), WORK( * )
239 * =====================================================================
242 DOUBLE PRECISION ZERO, ONE
243 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
244 COMPLEX*16 CZERO, CONE
245 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
246 $ CONE = ( 1.0D+0, 0.0D+0 ) )
248 * .. Local Scalars ..
249 LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
251 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
253 DOUBLE PRECISION ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
254 $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
255 $ SCALE, SMALL, TEMP, ULP, XMAX
256 COMPLEX*16 BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
258 * .. External Functions ..
260 DOUBLE PRECISION DLAMCH
262 EXTERNAL LSAME, DLAMCH, ZLADIV
264 * .. External Subroutines ..
265 EXTERNAL DLABAD, XERBLA, ZGEMV
267 * .. Intrinsic Functions ..
268 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
270 * .. Statement Functions ..
271 DOUBLE PRECISION ABS1
273 * .. Statement Function definitions ..
274 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
276 * .. Executable Statements ..
278 * Decode and Test the input parameters
280 IF( LSAME( HOWMNY, 'A' ) ) THEN
284 ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
288 ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
296 IF( LSAME( SIDE, 'R' ) ) THEN
300 ELSE IF( LSAME( SIDE, 'L' ) ) THEN
304 ELSE IF( LSAME( SIDE, 'B' ) ) THEN
313 IF( ISIDE.LT.0 ) THEN
315 ELSE IF( IHWMNY.LT.0 ) THEN
317 ELSE IF( N.LT.0 ) THEN
319 ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
321 ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
325 CALL XERBLA( 'ZTGEVC', -INFO )
329 * Count the number of eigenvectors
331 IF( .NOT.ILALL ) THEN
341 * Check diagonal of B
345 IF( DIMAG( P( J, J ) ).NE.ZERO )
351 ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
353 ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
355 ELSE IF( MM.LT.IM ) THEN
359 CALL XERBLA( 'ZTGEVC', -INFO )
363 * Quick return if possible
371 SAFMIN = DLAMCH( 'Safe minimum' )
373 CALL DLABAD( SAFMIN, BIG )
374 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
375 SMALL = SAFMIN*N / ULP
377 BIGNUM = ONE / ( SAFMIN*N )
379 * Compute the 1-norm of each column of the strictly upper triangular
380 * part of A and B to check for possible overflow in the triangular
383 ANORM = ABS1( S( 1, 1 ) )
384 BNORM = ABS1( P( 1, 1 ) )
391 RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) )
392 RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) )
394 ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) )
395 BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) )
398 ASCALE = ONE / MAX( ANORM, SAFMIN )
399 BSCALE = ONE / MAX( BNORM, SAFMIN )
406 * Main loop over eigenvalues
412 ILCOMP = SELECT( JE )
417 IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
418 $ ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
420 * Singular matrix pencil -- return unit eigenvector
423 VL( JR, IEIG ) = CZERO
425 VL( IEIG, IEIG ) = CONE
429 * Non-singular eigenvalue:
430 * Compute coefficients a and b in
432 * y ( a A - b B ) = 0
434 TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
435 $ ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
436 SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
437 SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
438 ACOEFF = SBETA*ASCALE
439 BCOEFF = SALPHA*BSCALE
441 * Scale to avoid underflow
443 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
444 LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
449 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
451 $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
452 $ MIN( BNORM, BIG ) )
453 IF( LSA .OR. LSB ) THEN
454 SCALE = MIN( SCALE, ONE /
455 $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
456 $ ABS1( BCOEFF ) ) ) )
458 ACOEFF = ASCALE*( SCALE*SBETA )
460 ACOEFF = SCALE*ACOEFF
463 BCOEFF = BSCALE*( SCALE*SALPHA )
465 BCOEFF = SCALE*BCOEFF
469 ACOEFA = ABS( ACOEFF )
470 BCOEFA = ABS1( BCOEFF )
476 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
479 * Triangular solve of (a A - b B) y = 0
482 * (rowwise in (a A - b B) , or columnwise in a A - b B)
488 * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
490 * (Scale if necessary)
493 IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM*
496 WORK( JR ) = TEMP*WORK( JR )
504 SUMA = SUMA + DCONJG( S( JR, J ) )*WORK( JR )
505 SUMB = SUMB + DCONJG( P( JR, J ) )*WORK( JR )
507 SUM = ACOEFF*SUMA - DCONJG( BCOEFF )*SUMB
509 * Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
511 * with scaling and perturbation of the denominator
513 D = DCONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
514 IF( ABS1( D ).LE.DMIN )
517 IF( ABS1( D ).LT.ONE ) THEN
518 IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN
519 TEMP = ONE / ABS1( SUM )
521 WORK( JR ) = TEMP*WORK( JR )
527 WORK( J ) = ZLADIV( -SUM, D )
528 XMAX = MAX( XMAX, ABS1( WORK( J ) ) )
531 * Back transform eigenvector if HOWMNY='B'.
534 CALL ZGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL,
535 $ WORK( JE ), 1, CZERO, WORK( N+1 ), 1 )
543 * Copy and scale eigenvector into column of VL
547 XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
550 IF( XMAX.GT.SAFMIN ) THEN
553 VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
559 DO 130 JR = 1, IBEG - 1
560 VL( JR, IEIG ) = CZERO
572 * Main loop over eigenvalues
578 ILCOMP = SELECT( JE )
583 IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
584 $ ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
586 * Singular matrix pencil -- return unit eigenvector
589 VR( JR, IEIG ) = CZERO
591 VR( IEIG, IEIG ) = CONE
595 * Non-singular eigenvalue:
596 * Compute coefficients a and b in
598 * ( a A - b B ) x = 0
600 TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
601 $ ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
602 SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
603 SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
604 ACOEFF = SBETA*ASCALE
605 BCOEFF = SALPHA*BSCALE
607 * Scale to avoid underflow
609 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
610 LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
615 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
617 $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
618 $ MIN( BNORM, BIG ) )
619 IF( LSA .OR. LSB ) THEN
620 SCALE = MIN( SCALE, ONE /
621 $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
622 $ ABS1( BCOEFF ) ) ) )
624 ACOEFF = ASCALE*( SCALE*SBETA )
626 ACOEFF = SCALE*ACOEFF
629 BCOEFF = BSCALE*( SCALE*SALPHA )
631 BCOEFF = SCALE*BCOEFF
635 ACOEFA = ABS( ACOEFF )
636 BCOEFA = ABS1( BCOEFF )
642 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
644 * Triangular solve of (a A - b B) x = 0 (columnwise)
646 * WORK(1:j-1) contains sums w,
647 * WORK(j+1:JE) contains x
649 DO 170 JR = 1, JE - 1
650 WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
654 DO 210 J = JE - 1, 1, -1
656 * Form x(j) := - w(j) / d
657 * with scaling and perturbation of the denominator
659 D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
660 IF( ABS1( D ).LE.DMIN )
663 IF( ABS1( D ).LT.ONE ) THEN
664 IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN
665 TEMP = ONE / ABS1( WORK( J ) )
667 WORK( JR ) = TEMP*WORK( JR )
672 WORK( J ) = ZLADIV( -WORK( J ), D )
676 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
678 IF( ABS1( WORK( J ) ).GT.ONE ) THEN
679 TEMP = ONE / ABS1( WORK( J ) )
680 IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE.
683 WORK( JR ) = TEMP*WORK( JR )
688 CA = ACOEFF*WORK( J )
689 CB = BCOEFF*WORK( J )
691 WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
697 * Back transform eigenvector if HOWMNY='B'.
700 CALL ZGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1,
701 $ CZERO, WORK( N+1 ), 1 )
709 * Copy and scale eigenvector into column of VR
713 XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
716 IF( XMAX.GT.SAFMIN ) THEN
719 VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
725 DO 240 JR = IEND + 1, N
726 VR( JR, IEIG ) = CZERO