Annotation of rpl/lapack/lapack/ztgevc.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b ZTGEVC
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZTGEVC + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgevc.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgevc.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgevc.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
! 22: * LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER HOWMNY, SIDE
! 26: * INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * LOGICAL SELECT( * )
! 30: * DOUBLE PRECISION RWORK( * )
! 31: * COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
! 32: * $ VR( LDVR, * ), WORK( * )
! 33: * ..
! 34: *
! 35: *
! 36: *
! 37: *> \par Purpose:
! 38: * =============
! 39: *>
! 40: *> \verbatim
! 41: *>
! 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):
! 46: *>
! 47: *> A = Q*S*Z**H, B = Q*P*Z**H
! 48: *>
! 49: *> as computed by ZGGHRD + ZHGEQZ.
! 50: *>
! 51: *> The right eigenvector x and the left eigenvector y of (S,P)
! 52: *> corresponding to an eigenvalue w are defined by:
! 53: *>
! 54: *> S*x = w*P*x, (y**H)*S = w*(y**H)*P,
! 55: *>
! 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.
! 59: *>
! 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).
! 66: *> \endverbatim
! 67: *
! 68: * Arguments:
! 69: * ==========
! 70: *
! 71: *> \param[in] SIDE
! 72: *> \verbatim
! 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.
! 77: *> \endverbatim
! 78: *>
! 79: *> \param[in] HOWMNY
! 80: *> \verbatim
! 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.
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] SELECT
! 90: *> \verbatim
! 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'.
! 96: *> \endverbatim
! 97: *>
! 98: *> \param[in] N
! 99: *> \verbatim
! 100: *> N is INTEGER
! 101: *> The order of the matrices S and P. N >= 0.
! 102: *> \endverbatim
! 103: *>
! 104: *> \param[in] S
! 105: *> \verbatim
! 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.
! 109: *> \endverbatim
! 110: *>
! 111: *> \param[in] LDS
! 112: *> \verbatim
! 113: *> LDS is INTEGER
! 114: *> The leading dimension of array S. LDS >= max(1,N).
! 115: *> \endverbatim
! 116: *>
! 117: *> \param[in] P
! 118: *> \verbatim
! 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.
! 123: *> \endverbatim
! 124: *>
! 125: *> \param[in] LDP
! 126: *> \verbatim
! 127: *> LDP is INTEGER
! 128: *> The leading dimension of array P. LDP >= max(1,N).
! 129: *> \endverbatim
! 130: *>
! 131: *> \param[in,out] VL
! 132: *> \verbatim
! 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'.
! 144: *> \endverbatim
! 145: *>
! 146: *> \param[in] LDVL
! 147: *> \verbatim
! 148: *> LDVL is INTEGER
! 149: *> The leading dimension of array VL. LDVL >= 1, and if
! 150: *> SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
! 151: *> \endverbatim
! 152: *>
! 153: *> \param[in,out] VR
! 154: *> \verbatim
! 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'.
! 166: *> \endverbatim
! 167: *>
! 168: *> \param[in] LDVR
! 169: *> \verbatim
! 170: *> LDVR is INTEGER
! 171: *> The leading dimension of the array VR. LDVR >= 1, and if
! 172: *> SIDE = 'R' or 'B', LDVR >= N.
! 173: *> \endverbatim
! 174: *>
! 175: *> \param[in] MM
! 176: *> \verbatim
! 177: *> MM is INTEGER
! 178: *> The number of columns in the arrays VL and/or VR. MM >= M.
! 179: *> \endverbatim
! 180: *>
! 181: *> \param[out] M
! 182: *> \verbatim
! 183: *> M is INTEGER
! 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.
! 187: *> \endverbatim
! 188: *>
! 189: *> \param[out] WORK
! 190: *> \verbatim
! 191: *> WORK is COMPLEX*16 array, dimension (2*N)
! 192: *> \endverbatim
! 193: *>
! 194: *> \param[out] RWORK
! 195: *> \verbatim
! 196: *> RWORK is DOUBLE PRECISION array, dimension (2*N)
! 197: *> \endverbatim
! 198: *>
! 199: *> \param[out] INFO
! 200: *> \verbatim
! 201: *> INFO is INTEGER
! 202: *> = 0: successful exit.
! 203: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 204: *> \endverbatim
! 205: *
! 206: * Authors:
! 207: * ========
! 208: *
! 209: *> \author Univ. of Tennessee
! 210: *> \author Univ. of California Berkeley
! 211: *> \author Univ. of Colorado Denver
! 212: *> \author NAG Ltd.
! 213: *
! 214: *> \date November 2011
! 215: *
! 216: *> \ingroup complex16GEcomputational
! 217: *
! 218: * =====================================================================
1.1 bertrand 219: SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
220: $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
221: *
1.8 ! bertrand 222: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 223: * -- LAPACK is a software package provided by Univ. of Tennessee, --
224: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 225: * November 2011
1.1 bertrand 226: *
227: * .. Scalar Arguments ..
228: CHARACTER HOWMNY, SIDE
229: INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
230: * ..
231: * .. Array Arguments ..
232: LOGICAL SELECT( * )
233: DOUBLE PRECISION RWORK( * )
234: COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
235: $ VR( LDVR, * ), WORK( * )
236: * ..
237: *
238: *
239: * =====================================================================
240: *
241: * .. Parameters ..
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 ) )
247: * ..
248: * .. Local Scalars ..
249: LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
250: $ LSA, LSB
251: INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
252: $ J, JE, JR
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
257: * ..
258: * .. External Functions ..
259: LOGICAL LSAME
260: DOUBLE PRECISION DLAMCH
261: COMPLEX*16 ZLADIV
262: EXTERNAL LSAME, DLAMCH, ZLADIV
263: * ..
264: * .. External Subroutines ..
265: EXTERNAL DLABAD, XERBLA, ZGEMV
266: * ..
267: * .. Intrinsic Functions ..
268: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
269: * ..
270: * .. Statement Functions ..
271: DOUBLE PRECISION ABS1
272: * ..
273: * .. Statement Function definitions ..
274: ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
275: * ..
276: * .. Executable Statements ..
277: *
278: * Decode and Test the input parameters
279: *
280: IF( LSAME( HOWMNY, 'A' ) ) THEN
281: IHWMNY = 1
282: ILALL = .TRUE.
283: ILBACK = .FALSE.
284: ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
285: IHWMNY = 2
286: ILALL = .FALSE.
287: ILBACK = .FALSE.
288: ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
289: IHWMNY = 3
290: ILALL = .TRUE.
291: ILBACK = .TRUE.
292: ELSE
293: IHWMNY = -1
294: END IF
295: *
296: IF( LSAME( SIDE, 'R' ) ) THEN
297: ISIDE = 1
298: COMPL = .FALSE.
299: COMPR = .TRUE.
300: ELSE IF( LSAME( SIDE, 'L' ) ) THEN
301: ISIDE = 2
302: COMPL = .TRUE.
303: COMPR = .FALSE.
304: ELSE IF( LSAME( SIDE, 'B' ) ) THEN
305: ISIDE = 3
306: COMPL = .TRUE.
307: COMPR = .TRUE.
308: ELSE
309: ISIDE = -1
310: END IF
311: *
312: INFO = 0
313: IF( ISIDE.LT.0 ) THEN
314: INFO = -1
315: ELSE IF( IHWMNY.LT.0 ) THEN
316: INFO = -2
317: ELSE IF( N.LT.0 ) THEN
318: INFO = -4
319: ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
320: INFO = -6
321: ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
322: INFO = -8
323: END IF
324: IF( INFO.NE.0 ) THEN
325: CALL XERBLA( 'ZTGEVC', -INFO )
326: RETURN
327: END IF
328: *
329: * Count the number of eigenvectors
330: *
331: IF( .NOT.ILALL ) THEN
332: IM = 0
333: DO 10 J = 1, N
334: IF( SELECT( J ) )
335: $ IM = IM + 1
336: 10 CONTINUE
337: ELSE
338: IM = N
339: END IF
340: *
341: * Check diagonal of B
342: *
343: ILBBAD = .FALSE.
344: DO 20 J = 1, N
345: IF( DIMAG( P( J, J ) ).NE.ZERO )
346: $ ILBBAD = .TRUE.
347: 20 CONTINUE
348: *
349: IF( ILBBAD ) THEN
350: INFO = -7
351: ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
352: INFO = -10
353: ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
354: INFO = -12
355: ELSE IF( MM.LT.IM ) THEN
356: INFO = -13
357: END IF
358: IF( INFO.NE.0 ) THEN
359: CALL XERBLA( 'ZTGEVC', -INFO )
360: RETURN
361: END IF
362: *
363: * Quick return if possible
364: *
365: M = IM
366: IF( N.EQ.0 )
367: $ RETURN
368: *
369: * Machine Constants
370: *
371: SAFMIN = DLAMCH( 'Safe minimum' )
372: BIG = ONE / SAFMIN
373: CALL DLABAD( SAFMIN, BIG )
374: ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
375: SMALL = SAFMIN*N / ULP
376: BIG = ONE / SMALL
377: BIGNUM = ONE / ( SAFMIN*N )
378: *
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
381: * solver.
382: *
383: ANORM = ABS1( S( 1, 1 ) )
384: BNORM = ABS1( P( 1, 1 ) )
385: RWORK( 1 ) = ZERO
386: RWORK( N+1 ) = ZERO
387: DO 40 J = 2, N
388: RWORK( J ) = ZERO
389: RWORK( N+J ) = ZERO
390: DO 30 I = 1, J - 1
391: RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) )
392: RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) )
393: 30 CONTINUE
394: ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) )
395: BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) )
396: 40 CONTINUE
397: *
398: ASCALE = ONE / MAX( ANORM, SAFMIN )
399: BSCALE = ONE / MAX( BNORM, SAFMIN )
400: *
401: * Left eigenvectors
402: *
403: IF( COMPL ) THEN
404: IEIG = 0
405: *
406: * Main loop over eigenvalues
407: *
408: DO 140 JE = 1, N
409: IF( ILALL ) THEN
410: ILCOMP = .TRUE.
411: ELSE
412: ILCOMP = SELECT( JE )
413: END IF
414: IF( ILCOMP ) THEN
415: IEIG = IEIG + 1
416: *
417: IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
418: $ ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
419: *
420: * Singular matrix pencil -- return unit eigenvector
421: *
422: DO 50 JR = 1, N
423: VL( JR, IEIG ) = CZERO
424: 50 CONTINUE
425: VL( IEIG, IEIG ) = CONE
426: GO TO 140
427: END IF
428: *
429: * Non-singular eigenvalue:
430: * Compute coefficients a and b in
431: * H
432: * y ( a A - b B ) = 0
433: *
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
440: *
441: * Scale to avoid underflow
442: *
443: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
444: LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
445: $ SMALL
446: *
447: SCALE = ONE
448: IF( LSA )
449: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
450: IF( LSB )
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 ) ) ) )
457: IF( LSA ) THEN
458: ACOEFF = ASCALE*( SCALE*SBETA )
459: ELSE
460: ACOEFF = SCALE*ACOEFF
461: END IF
462: IF( LSB ) THEN
463: BCOEFF = BSCALE*( SCALE*SALPHA )
464: ELSE
465: BCOEFF = SCALE*BCOEFF
466: END IF
467: END IF
468: *
469: ACOEFA = ABS( ACOEFF )
470: BCOEFA = ABS1( BCOEFF )
471: XMAX = ONE
472: DO 60 JR = 1, N
473: WORK( JR ) = CZERO
474: 60 CONTINUE
475: WORK( JE ) = CONE
476: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
477: *
478: * H
479: * Triangular solve of (a A - b B) y = 0
480: *
481: * H
482: * (rowwise in (a A - b B) , or columnwise in a A - b B)
483: *
484: DO 100 J = JE + 1, N
485: *
486: * Compute
487: * j-1
488: * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
489: * k=je
490: * (Scale if necessary)
491: *
492: TEMP = ONE / XMAX
493: IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM*
494: $ TEMP ) THEN
495: DO 70 JR = JE, J - 1
496: WORK( JR ) = TEMP*WORK( JR )
497: 70 CONTINUE
498: XMAX = ONE
499: END IF
500: SUMA = CZERO
501: SUMB = CZERO
502: *
503: DO 80 JR = JE, J - 1
504: SUMA = SUMA + DCONJG( S( JR, J ) )*WORK( JR )
505: SUMB = SUMB + DCONJG( P( JR, J ) )*WORK( JR )
506: 80 CONTINUE
507: SUM = ACOEFF*SUMA - DCONJG( BCOEFF )*SUMB
508: *
509: * Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
510: *
511: * with scaling and perturbation of the denominator
512: *
513: D = DCONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
514: IF( ABS1( D ).LE.DMIN )
515: $ D = DCMPLX( DMIN )
516: *
517: IF( ABS1( D ).LT.ONE ) THEN
518: IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN
519: TEMP = ONE / ABS1( SUM )
520: DO 90 JR = JE, J - 1
521: WORK( JR ) = TEMP*WORK( JR )
522: 90 CONTINUE
523: XMAX = TEMP*XMAX
524: SUM = TEMP*SUM
525: END IF
526: END IF
527: WORK( J ) = ZLADIV( -SUM, D )
528: XMAX = MAX( XMAX, ABS1( WORK( J ) ) )
529: 100 CONTINUE
530: *
531: * Back transform eigenvector if HOWMNY='B'.
532: *
533: IF( ILBACK ) THEN
534: CALL ZGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL,
535: $ WORK( JE ), 1, CZERO, WORK( N+1 ), 1 )
536: ISRC = 2
537: IBEG = 1
538: ELSE
539: ISRC = 1
540: IBEG = JE
541: END IF
542: *
543: * Copy and scale eigenvector into column of VL
544: *
545: XMAX = ZERO
546: DO 110 JR = IBEG, N
547: XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
548: 110 CONTINUE
549: *
550: IF( XMAX.GT.SAFMIN ) THEN
551: TEMP = ONE / XMAX
552: DO 120 JR = IBEG, N
553: VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
554: 120 CONTINUE
555: ELSE
556: IBEG = N + 1
557: END IF
558: *
559: DO 130 JR = 1, IBEG - 1
560: VL( JR, IEIG ) = CZERO
561: 130 CONTINUE
562: *
563: END IF
564: 140 CONTINUE
565: END IF
566: *
567: * Right eigenvectors
568: *
569: IF( COMPR ) THEN
570: IEIG = IM + 1
571: *
572: * Main loop over eigenvalues
573: *
574: DO 250 JE = N, 1, -1
575: IF( ILALL ) THEN
576: ILCOMP = .TRUE.
577: ELSE
578: ILCOMP = SELECT( JE )
579: END IF
580: IF( ILCOMP ) THEN
581: IEIG = IEIG - 1
582: *
583: IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
584: $ ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
585: *
586: * Singular matrix pencil -- return unit eigenvector
587: *
588: DO 150 JR = 1, N
589: VR( JR, IEIG ) = CZERO
590: 150 CONTINUE
591: VR( IEIG, IEIG ) = CONE
592: GO TO 250
593: END IF
594: *
595: * Non-singular eigenvalue:
596: * Compute coefficients a and b in
597: *
598: * ( a A - b B ) x = 0
599: *
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
606: *
607: * Scale to avoid underflow
608: *
609: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
610: LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
611: $ SMALL
612: *
613: SCALE = ONE
614: IF( LSA )
615: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
616: IF( LSB )
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 ) ) ) )
623: IF( LSA ) THEN
624: ACOEFF = ASCALE*( SCALE*SBETA )
625: ELSE
626: ACOEFF = SCALE*ACOEFF
627: END IF
628: IF( LSB ) THEN
629: BCOEFF = BSCALE*( SCALE*SALPHA )
630: ELSE
631: BCOEFF = SCALE*BCOEFF
632: END IF
633: END IF
634: *
635: ACOEFA = ABS( ACOEFF )
636: BCOEFA = ABS1( BCOEFF )
637: XMAX = ONE
638: DO 160 JR = 1, N
639: WORK( JR ) = CZERO
640: 160 CONTINUE
641: WORK( JE ) = CONE
642: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
643: *
644: * Triangular solve of (a A - b B) x = 0 (columnwise)
645: *
646: * WORK(1:j-1) contains sums w,
647: * WORK(j+1:JE) contains x
648: *
649: DO 170 JR = 1, JE - 1
650: WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
651: 170 CONTINUE
652: WORK( JE ) = CONE
653: *
654: DO 210 J = JE - 1, 1, -1
655: *
656: * Form x(j) := - w(j) / d
657: * with scaling and perturbation of the denominator
658: *
659: D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
660: IF( ABS1( D ).LE.DMIN )
661: $ D = DCMPLX( DMIN )
662: *
663: IF( ABS1( D ).LT.ONE ) THEN
664: IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN
665: TEMP = ONE / ABS1( WORK( J ) )
666: DO 180 JR = 1, JE
667: WORK( JR ) = TEMP*WORK( JR )
668: 180 CONTINUE
669: END IF
670: END IF
671: *
672: WORK( J ) = ZLADIV( -WORK( J ), D )
673: *
674: IF( J.GT.1 ) THEN
675: *
676: * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
677: *
678: IF( ABS1( WORK( J ) ).GT.ONE ) THEN
679: TEMP = ONE / ABS1( WORK( J ) )
680: IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE.
681: $ BIGNUM*TEMP ) THEN
682: DO 190 JR = 1, JE
683: WORK( JR ) = TEMP*WORK( JR )
684: 190 CONTINUE
685: END IF
686: END IF
687: *
688: CA = ACOEFF*WORK( J )
689: CB = BCOEFF*WORK( J )
690: DO 200 JR = 1, J - 1
691: WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
692: $ CB*P( JR, J )
693: 200 CONTINUE
694: END IF
695: 210 CONTINUE
696: *
697: * Back transform eigenvector if HOWMNY='B'.
698: *
699: IF( ILBACK ) THEN
700: CALL ZGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1,
701: $ CZERO, WORK( N+1 ), 1 )
702: ISRC = 2
703: IEND = N
704: ELSE
705: ISRC = 1
706: IEND = JE
707: END IF
708: *
709: * Copy and scale eigenvector into column of VR
710: *
711: XMAX = ZERO
712: DO 220 JR = 1, IEND
713: XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
714: 220 CONTINUE
715: *
716: IF( XMAX.GT.SAFMIN ) THEN
717: TEMP = ONE / XMAX
718: DO 230 JR = 1, IEND
719: VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
720: 230 CONTINUE
721: ELSE
722: IEND = 0
723: END IF
724: *
725: DO 240 JR = IEND + 1, N
726: VR( JR, IEIG ) = CZERO
727: 240 CONTINUE
728: *
729: END IF
730: 250 CONTINUE
731: END IF
732: *
733: RETURN
734: *
735: * End of ZTGEVC
736: *
737: END
CVSweb interface <joel.bertrand@systella.fr>