File:  [local] / rpl / lapack / lapack / ztgevc.f
Revision 1.3: download - view: text, annotated - select for diffs - revision graph
Fri Aug 6 15:29:02 2010 UTC (13 years, 9 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Cohérence

    1:       SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
    2:      $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
    3: *
    4: *  -- LAPACK routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       CHARACTER          HOWMNY, SIDE
   11:       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
   12: *     ..
   13: *     .. Array Arguments ..
   14:       LOGICAL            SELECT( * )
   15:       DOUBLE PRECISION   RWORK( * )
   16:       COMPLEX*16         P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
   17:      $                   VR( LDVR, * ), WORK( * )
   18: *     ..
   19: *
   20: *
   21: *  Purpose
   22: *  =======
   23: *
   24: *  ZTGEVC computes some or all of the right and/or left eigenvectors of
   25: *  a pair of complex matrices (S,P), where S and P are upper triangular.
   26: *  Matrix pairs of this type are produced by the generalized Schur
   27: *  factorization of a complex matrix pair (A,B):
   28: *  
   29: *     A = Q*S*Z**H,  B = Q*P*Z**H
   30: *  
   31: *  as computed by ZGGHRD + ZHGEQZ.
   32: *  
   33: *  The right eigenvector x and the left eigenvector y of (S,P)
   34: *  corresponding to an eigenvalue w are defined by:
   35: *  
   36: *     S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
   37: *  
   38: *  where y**H denotes the conjugate tranpose of y.
   39: *  The eigenvalues are not input to this routine, but are computed
   40: *  directly from the diagonal elements of S and P.
   41: *  
   42: *  This routine returns the matrices X and/or Y of right and left
   43: *  eigenvectors of (S,P), or the products Z*X and/or Q*Y,
   44: *  where Z and Q are input matrices.
   45: *  If Q and Z are the unitary factors from the generalized Schur
   46: *  factorization of a matrix pair (A,B), then Z*X and Q*Y
   47: *  are the matrices of right and left eigenvectors of (A,B).
   48: *
   49: *  Arguments
   50: *  =========
   51: *
   52: *  SIDE    (input) CHARACTER*1
   53: *          = 'R': compute right eigenvectors only;
   54: *          = 'L': compute left eigenvectors only;
   55: *          = 'B': compute both right and left eigenvectors.
   56: *
   57: *  HOWMNY  (input) CHARACTER*1
   58: *          = 'A': compute all right and/or left eigenvectors;
   59: *          = 'B': compute all right and/or left eigenvectors,
   60: *                 backtransformed by the matrices in VR and/or VL;
   61: *          = 'S': compute selected right and/or left eigenvectors,
   62: *                 specified by the logical array SELECT.
   63: *
   64: *  SELECT  (input) LOGICAL array, dimension (N)
   65: *          If HOWMNY='S', SELECT specifies the eigenvectors to be
   66: *          computed.  The eigenvector corresponding to the j-th
   67: *          eigenvalue is computed if SELECT(j) = .TRUE..
   68: *          Not referenced if HOWMNY = 'A' or 'B'.
   69: *
   70: *  N       (input) INTEGER
   71: *          The order of the matrices S and P.  N >= 0.
   72: *
   73: *  S       (input) COMPLEX*16 array, dimension (LDS,N)
   74: *          The upper triangular matrix S from a generalized Schur
   75: *          factorization, as computed by ZHGEQZ.
   76: *
   77: *  LDS     (input) INTEGER
   78: *          The leading dimension of array S.  LDS >= max(1,N).
   79: *
   80: *  P       (input) COMPLEX*16 array, dimension (LDP,N)
   81: *          The upper triangular matrix P from a generalized Schur
   82: *          factorization, as computed by ZHGEQZ.  P must have real
   83: *          diagonal elements.
   84: *
   85: *  LDP     (input) INTEGER
   86: *          The leading dimension of array P.  LDP >= max(1,N).
   87: *
   88: *  VL      (input/output) COMPLEX*16 array, dimension (LDVL,MM)
   89: *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
   90: *          contain an N-by-N matrix Q (usually the unitary matrix Q
   91: *          of left Schur vectors returned by ZHGEQZ).
   92: *          On exit, if SIDE = 'L' or 'B', VL contains:
   93: *          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
   94: *          if HOWMNY = 'B', the matrix Q*Y;
   95: *          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
   96: *                      SELECT, stored consecutively in the columns of
   97: *                      VL, in the same order as their eigenvalues.
   98: *          Not referenced if SIDE = 'R'.
   99: *
  100: *  LDVL    (input) INTEGER
  101: *          The leading dimension of array VL.  LDVL >= 1, and if
  102: *          SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
  103: *
  104: *  VR      (input/output) COMPLEX*16 array, dimension (LDVR,MM)
  105: *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
  106: *          contain an N-by-N matrix Q (usually the unitary matrix Z
  107: *          of right Schur vectors returned by ZHGEQZ).
  108: *          On exit, if SIDE = 'R' or 'B', VR contains:
  109: *          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
  110: *          if HOWMNY = 'B', the matrix Z*X;
  111: *          if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
  112: *                      SELECT, stored consecutively in the columns of
  113: *                      VR, in the same order as their eigenvalues.
  114: *          Not referenced if SIDE = 'L'.
  115: *
  116: *  LDVR    (input) INTEGER
  117: *          The leading dimension of the array VR.  LDVR >= 1, and if
  118: *          SIDE = 'R' or 'B', LDVR >= N.
  119: *
  120: *  MM      (input) INTEGER
  121: *          The number of columns in the arrays VL and/or VR. MM >= M.
  122: *
  123: *  M       (output) INTEGER
  124: *          The number of columns in the arrays VL and/or VR actually
  125: *          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
  126: *          is set to N.  Each selected eigenvector occupies one column.
  127: *
  128: *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
  129: *
  130: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
  131: *
  132: *  INFO    (output) INTEGER
  133: *          = 0:  successful exit.
  134: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  135: *
  136: *  =====================================================================
  137: *
  138: *     .. Parameters ..
  139:       DOUBLE PRECISION   ZERO, ONE
  140:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  141:       COMPLEX*16         CZERO, CONE
  142:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
  143:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
  144: *     ..
  145: *     .. Local Scalars ..
  146:       LOGICAL            COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
  147:      $                   LSA, LSB
  148:       INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
  149:      $                   J, JE, JR
  150:       DOUBLE PRECISION   ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
  151:      $                   BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
  152:      $                   SCALE, SMALL, TEMP, ULP, XMAX
  153:       COMPLEX*16         BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
  154: *     ..
  155: *     .. External Functions ..
  156:       LOGICAL            LSAME
  157:       DOUBLE PRECISION   DLAMCH
  158:       COMPLEX*16         ZLADIV
  159:       EXTERNAL           LSAME, DLAMCH, ZLADIV
  160: *     ..
  161: *     .. External Subroutines ..
  162:       EXTERNAL           DLABAD, XERBLA, ZGEMV
  163: *     ..
  164: *     .. Intrinsic Functions ..
  165:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
  166: *     ..
  167: *     .. Statement Functions ..
  168:       DOUBLE PRECISION   ABS1
  169: *     ..
  170: *     .. Statement Function definitions ..
  171:       ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
  172: *     ..
  173: *     .. Executable Statements ..
  174: *
  175: *     Decode and Test the input parameters
  176: *
  177:       IF( LSAME( HOWMNY, 'A' ) ) THEN
  178:          IHWMNY = 1
  179:          ILALL = .TRUE.
  180:          ILBACK = .FALSE.
  181:       ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
  182:          IHWMNY = 2
  183:          ILALL = .FALSE.
  184:          ILBACK = .FALSE.
  185:       ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
  186:          IHWMNY = 3
  187:          ILALL = .TRUE.
  188:          ILBACK = .TRUE.
  189:       ELSE
  190:          IHWMNY = -1
  191:       END IF
  192: *
  193:       IF( LSAME( SIDE, 'R' ) ) THEN
  194:          ISIDE = 1
  195:          COMPL = .FALSE.
  196:          COMPR = .TRUE.
  197:       ELSE IF( LSAME( SIDE, 'L' ) ) THEN
  198:          ISIDE = 2
  199:          COMPL = .TRUE.
  200:          COMPR = .FALSE.
  201:       ELSE IF( LSAME( SIDE, 'B' ) ) THEN
  202:          ISIDE = 3
  203:          COMPL = .TRUE.
  204:          COMPR = .TRUE.
  205:       ELSE
  206:          ISIDE = -1
  207:       END IF
  208: *
  209:       INFO = 0
  210:       IF( ISIDE.LT.0 ) THEN
  211:          INFO = -1
  212:       ELSE IF( IHWMNY.LT.0 ) THEN
  213:          INFO = -2
  214:       ELSE IF( N.LT.0 ) THEN
  215:          INFO = -4
  216:       ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
  217:          INFO = -6
  218:       ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
  219:          INFO = -8
  220:       END IF
  221:       IF( INFO.NE.0 ) THEN
  222:          CALL XERBLA( 'ZTGEVC', -INFO )
  223:          RETURN
  224:       END IF
  225: *
  226: *     Count the number of eigenvectors
  227: *
  228:       IF( .NOT.ILALL ) THEN
  229:          IM = 0
  230:          DO 10 J = 1, N
  231:             IF( SELECT( J ) )
  232:      $         IM = IM + 1
  233:    10    CONTINUE
  234:       ELSE
  235:          IM = N
  236:       END IF
  237: *
  238: *     Check diagonal of B
  239: *
  240:       ILBBAD = .FALSE.
  241:       DO 20 J = 1, N
  242:          IF( DIMAG( P( J, J ) ).NE.ZERO )
  243:      $      ILBBAD = .TRUE.
  244:    20 CONTINUE
  245: *
  246:       IF( ILBBAD ) THEN
  247:          INFO = -7
  248:       ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
  249:          INFO = -10
  250:       ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
  251:          INFO = -12
  252:       ELSE IF( MM.LT.IM ) THEN
  253:          INFO = -13
  254:       END IF
  255:       IF( INFO.NE.0 ) THEN
  256:          CALL XERBLA( 'ZTGEVC', -INFO )
  257:          RETURN
  258:       END IF
  259: *
  260: *     Quick return if possible
  261: *
  262:       M = IM
  263:       IF( N.EQ.0 )
  264:      $   RETURN
  265: *
  266: *     Machine Constants
  267: *
  268:       SAFMIN = DLAMCH( 'Safe minimum' )
  269:       BIG = ONE / SAFMIN
  270:       CALL DLABAD( SAFMIN, BIG )
  271:       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
  272:       SMALL = SAFMIN*N / ULP
  273:       BIG = ONE / SMALL
  274:       BIGNUM = ONE / ( SAFMIN*N )
  275: *
  276: *     Compute the 1-norm of each column of the strictly upper triangular
  277: *     part of A and B to check for possible overflow in the triangular
  278: *     solver.
  279: *
  280:       ANORM = ABS1( S( 1, 1 ) )
  281:       BNORM = ABS1( P( 1, 1 ) )
  282:       RWORK( 1 ) = ZERO
  283:       RWORK( N+1 ) = ZERO
  284:       DO 40 J = 2, N
  285:          RWORK( J ) = ZERO
  286:          RWORK( N+J ) = ZERO
  287:          DO 30 I = 1, J - 1
  288:             RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) )
  289:             RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) )
  290:    30    CONTINUE
  291:          ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) )
  292:          BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) )
  293:    40 CONTINUE
  294: *
  295:       ASCALE = ONE / MAX( ANORM, SAFMIN )
  296:       BSCALE = ONE / MAX( BNORM, SAFMIN )
  297: *
  298: *     Left eigenvectors
  299: *
  300:       IF( COMPL ) THEN
  301:          IEIG = 0
  302: *
  303: *        Main loop over eigenvalues
  304: *
  305:          DO 140 JE = 1, N
  306:             IF( ILALL ) THEN
  307:                ILCOMP = .TRUE.
  308:             ELSE
  309:                ILCOMP = SELECT( JE )
  310:             END IF
  311:             IF( ILCOMP ) THEN
  312:                IEIG = IEIG + 1
  313: *
  314:                IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
  315:      $             ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
  316: *
  317: *                 Singular matrix pencil -- return unit eigenvector
  318: *
  319:                   DO 50 JR = 1, N
  320:                      VL( JR, IEIG ) = CZERO
  321:    50             CONTINUE
  322:                   VL( IEIG, IEIG ) = CONE
  323:                   GO TO 140
  324:                END IF
  325: *
  326: *              Non-singular eigenvalue:
  327: *              Compute coefficients  a  and  b  in
  328: *                   H
  329: *                 y  ( a A - b B ) = 0
  330: *
  331:                TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
  332:      $                ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
  333:                SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
  334:                SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
  335:                ACOEFF = SBETA*ASCALE
  336:                BCOEFF = SALPHA*BSCALE
  337: *
  338: *              Scale to avoid underflow
  339: *
  340:                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
  341:                LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
  342:      $               SMALL
  343: *
  344:                SCALE = ONE
  345:                IF( LSA )
  346:      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
  347:                IF( LSB )
  348:      $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
  349:      $                    MIN( BNORM, BIG ) )
  350:                IF( LSA .OR. LSB ) THEN
  351:                   SCALE = MIN( SCALE, ONE /
  352:      $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
  353:      $                    ABS1( BCOEFF ) ) ) )
  354:                   IF( LSA ) THEN
  355:                      ACOEFF = ASCALE*( SCALE*SBETA )
  356:                   ELSE
  357:                      ACOEFF = SCALE*ACOEFF
  358:                   END IF
  359:                   IF( LSB ) THEN
  360:                      BCOEFF = BSCALE*( SCALE*SALPHA )
  361:                   ELSE
  362:                      BCOEFF = SCALE*BCOEFF
  363:                   END IF
  364:                END IF
  365: *
  366:                ACOEFA = ABS( ACOEFF )
  367:                BCOEFA = ABS1( BCOEFF )
  368:                XMAX = ONE
  369:                DO 60 JR = 1, N
  370:                   WORK( JR ) = CZERO
  371:    60          CONTINUE
  372:                WORK( JE ) = CONE
  373:                DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
  374: *
  375: *                                              H
  376: *              Triangular solve of  (a A - b B)  y = 0
  377: *
  378: *                                      H
  379: *              (rowwise in  (a A - b B) , or columnwise in a A - b B)
  380: *
  381:                DO 100 J = JE + 1, N
  382: *
  383: *                 Compute
  384: *                       j-1
  385: *                 SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
  386: *                       k=je
  387: *                 (Scale if necessary)
  388: *
  389:                   TEMP = ONE / XMAX
  390:                   IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM*
  391:      $                TEMP ) THEN
  392:                      DO 70 JR = JE, J - 1
  393:                         WORK( JR ) = TEMP*WORK( JR )
  394:    70                CONTINUE
  395:                      XMAX = ONE
  396:                   END IF
  397:                   SUMA = CZERO
  398:                   SUMB = CZERO
  399: *
  400:                   DO 80 JR = JE, J - 1
  401:                      SUMA = SUMA + DCONJG( S( JR, J ) )*WORK( JR )
  402:                      SUMB = SUMB + DCONJG( P( JR, J ) )*WORK( JR )
  403:    80             CONTINUE
  404:                   SUM = ACOEFF*SUMA - DCONJG( BCOEFF )*SUMB
  405: *
  406: *                 Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
  407: *
  408: *                 with scaling and perturbation of the denominator
  409: *
  410:                   D = DCONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
  411:                   IF( ABS1( D ).LE.DMIN )
  412:      $               D = DCMPLX( DMIN )
  413: *
  414:                   IF( ABS1( D ).LT.ONE ) THEN
  415:                      IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN
  416:                         TEMP = ONE / ABS1( SUM )
  417:                         DO 90 JR = JE, J - 1
  418:                            WORK( JR ) = TEMP*WORK( JR )
  419:    90                   CONTINUE
  420:                         XMAX = TEMP*XMAX
  421:                         SUM = TEMP*SUM
  422:                      END IF
  423:                   END IF
  424:                   WORK( J ) = ZLADIV( -SUM, D )
  425:                   XMAX = MAX( XMAX, ABS1( WORK( J ) ) )
  426:   100          CONTINUE
  427: *
  428: *              Back transform eigenvector if HOWMNY='B'.
  429: *
  430:                IF( ILBACK ) THEN
  431:                   CALL ZGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL,
  432:      $                        WORK( JE ), 1, CZERO, WORK( N+1 ), 1 )
  433:                   ISRC = 2
  434:                   IBEG = 1
  435:                ELSE
  436:                   ISRC = 1
  437:                   IBEG = JE
  438:                END IF
  439: *
  440: *              Copy and scale eigenvector into column of VL
  441: *
  442:                XMAX = ZERO
  443:                DO 110 JR = IBEG, N
  444:                   XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
  445:   110          CONTINUE
  446: *
  447:                IF( XMAX.GT.SAFMIN ) THEN
  448:                   TEMP = ONE / XMAX
  449:                   DO 120 JR = IBEG, N
  450:                      VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
  451:   120             CONTINUE
  452:                ELSE
  453:                   IBEG = N + 1
  454:                END IF
  455: *
  456:                DO 130 JR = 1, IBEG - 1
  457:                   VL( JR, IEIG ) = CZERO
  458:   130          CONTINUE
  459: *
  460:             END IF
  461:   140    CONTINUE
  462:       END IF
  463: *
  464: *     Right eigenvectors
  465: *
  466:       IF( COMPR ) THEN
  467:          IEIG = IM + 1
  468: *
  469: *        Main loop over eigenvalues
  470: *
  471:          DO 250 JE = N, 1, -1
  472:             IF( ILALL ) THEN
  473:                ILCOMP = .TRUE.
  474:             ELSE
  475:                ILCOMP = SELECT( JE )
  476:             END IF
  477:             IF( ILCOMP ) THEN
  478:                IEIG = IEIG - 1
  479: *
  480:                IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
  481:      $             ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
  482: *
  483: *                 Singular matrix pencil -- return unit eigenvector
  484: *
  485:                   DO 150 JR = 1, N
  486:                      VR( JR, IEIG ) = CZERO
  487:   150             CONTINUE
  488:                   VR( IEIG, IEIG ) = CONE
  489:                   GO TO 250
  490:                END IF
  491: *
  492: *              Non-singular eigenvalue:
  493: *              Compute coefficients  a  and  b  in
  494: *
  495: *              ( a A - b B ) x  = 0
  496: *
  497:                TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
  498:      $                ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
  499:                SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
  500:                SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
  501:                ACOEFF = SBETA*ASCALE
  502:                BCOEFF = SALPHA*BSCALE
  503: *
  504: *              Scale to avoid underflow
  505: *
  506:                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
  507:                LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
  508:      $               SMALL
  509: *
  510:                SCALE = ONE
  511:                IF( LSA )
  512:      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
  513:                IF( LSB )
  514:      $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
  515:      $                    MIN( BNORM, BIG ) )
  516:                IF( LSA .OR. LSB ) THEN
  517:                   SCALE = MIN( SCALE, ONE /
  518:      $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
  519:      $                    ABS1( BCOEFF ) ) ) )
  520:                   IF( LSA ) THEN
  521:                      ACOEFF = ASCALE*( SCALE*SBETA )
  522:                   ELSE
  523:                      ACOEFF = SCALE*ACOEFF
  524:                   END IF
  525:                   IF( LSB ) THEN
  526:                      BCOEFF = BSCALE*( SCALE*SALPHA )
  527:                   ELSE
  528:                      BCOEFF = SCALE*BCOEFF
  529:                   END IF
  530:                END IF
  531: *
  532:                ACOEFA = ABS( ACOEFF )
  533:                BCOEFA = ABS1( BCOEFF )
  534:                XMAX = ONE
  535:                DO 160 JR = 1, N
  536:                   WORK( JR ) = CZERO
  537:   160          CONTINUE
  538:                WORK( JE ) = CONE
  539:                DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
  540: *
  541: *              Triangular solve of  (a A - b B) x = 0  (columnwise)
  542: *
  543: *              WORK(1:j-1) contains sums w,
  544: *              WORK(j+1:JE) contains x
  545: *
  546:                DO 170 JR = 1, JE - 1
  547:                   WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
  548:   170          CONTINUE
  549:                WORK( JE ) = CONE
  550: *
  551:                DO 210 J = JE - 1, 1, -1
  552: *
  553: *                 Form x(j) := - w(j) / d
  554: *                 with scaling and perturbation of the denominator
  555: *
  556:                   D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
  557:                   IF( ABS1( D ).LE.DMIN )
  558:      $               D = DCMPLX( DMIN )
  559: *
  560:                   IF( ABS1( D ).LT.ONE ) THEN
  561:                      IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN
  562:                         TEMP = ONE / ABS1( WORK( J ) )
  563:                         DO 180 JR = 1, JE
  564:                            WORK( JR ) = TEMP*WORK( JR )
  565:   180                   CONTINUE
  566:                      END IF
  567:                   END IF
  568: *
  569:                   WORK( J ) = ZLADIV( -WORK( J ), D )
  570: *
  571:                   IF( J.GT.1 ) THEN
  572: *
  573: *                    w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
  574: *
  575:                      IF( ABS1( WORK( J ) ).GT.ONE ) THEN
  576:                         TEMP = ONE / ABS1( WORK( J ) )
  577:                         IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE.
  578:      $                      BIGNUM*TEMP ) THEN
  579:                            DO 190 JR = 1, JE
  580:                               WORK( JR ) = TEMP*WORK( JR )
  581:   190                      CONTINUE
  582:                         END IF
  583:                      END IF
  584: *
  585:                      CA = ACOEFF*WORK( J )
  586:                      CB = BCOEFF*WORK( J )
  587:                      DO 200 JR = 1, J - 1
  588:                         WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
  589:      $                               CB*P( JR, J )
  590:   200                CONTINUE
  591:                   END IF
  592:   210          CONTINUE
  593: *
  594: *              Back transform eigenvector if HOWMNY='B'.
  595: *
  596:                IF( ILBACK ) THEN
  597:                   CALL ZGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1,
  598:      $                        CZERO, WORK( N+1 ), 1 )
  599:                   ISRC = 2
  600:                   IEND = N
  601:                ELSE
  602:                   ISRC = 1
  603:                   IEND = JE
  604:                END IF
  605: *
  606: *              Copy and scale eigenvector into column of VR
  607: *
  608:                XMAX = ZERO
  609:                DO 220 JR = 1, IEND
  610:                   XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
  611:   220          CONTINUE
  612: *
  613:                IF( XMAX.GT.SAFMIN ) THEN
  614:                   TEMP = ONE / XMAX
  615:                   DO 230 JR = 1, IEND
  616:                      VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
  617:   230             CONTINUE
  618:                ELSE
  619:                   IEND = 0
  620:                END IF
  621: *
  622:                DO 240 JR = IEND + 1, N
  623:                   VR( JR, IEIG ) = CZERO
  624:   240          CONTINUE
  625: *
  626:             END IF
  627:   250    CONTINUE
  628:       END IF
  629: *
  630:       RETURN
  631: *
  632: *     End of ZTGEVC
  633: *
  634:       END

CVSweb interface <joel.bertrand@systella.fr>