ztgevc.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 658 行 · 第 1/3 页

         COMPL = .TRUE.
         COMPR = .TRUE.
      ELSE
         ISIDE = -1
      END IF
<span class="comment">*</span><span class="comment">
</span>      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( LDS.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
         INFO = -8
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.221"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZTGEVC.221"></a><a href="ztgevc.f.html#ZTGEVC.1">ZTGEVC</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Count the number of eigenvectors
</span><span class="comment">*</span><span class="comment">
</span>      IF( .NOT.ILALL ) THEN
         IM = 0
         DO 10 J = 1, N
            IF( SELECT( J ) )
     $         IM = IM + 1
   10    CONTINUE
      ELSE
         IM = N
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check diagonal of B
</span><span class="comment">*</span><span class="comment">
</span>      ILBBAD = .FALSE.
      DO 20 J = 1, N
         IF( DIMAG( P( J, J ) ).NE.ZERO )
     $      ILBBAD = .TRUE.
   20 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      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 <a name="XERBLA.255"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZTGEVC.255"></a><a href="ztgevc.f.html#ZTGEVC.1">ZTGEVC</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      M = IM
      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Machine Constants
</span><span class="comment">*</span><span class="comment">
</span>      SAFMIN = <a name="DLAMCH.267"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
      BIG = ONE / SAFMIN
      CALL <a name="DLABAD.269"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>( SAFMIN, BIG )
      ULP = <a name="DLAMCH.270"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> )*<a name="DLAMCH.270"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Base'</span> )
      SMALL = SAFMIN*N / ULP
      BIG = ONE / SMALL
      BIGNUM = ONE / ( SAFMIN*N )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the 1-norm of each column of the strictly upper triangular
</span><span class="comment">*</span><span class="comment">     part of A and B to check for possible overflow in the triangular
</span><span class="comment">*</span><span class="comment">     solver.
</span><span class="comment">*</span><span class="comment">
</span>      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( S( I, J ) )
            RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) )
   30    CONTINUE
         ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) )
         BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) )
   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      ASCALE = ONE / MAX( ANORM, SAFMIN )
      BSCALE = ONE / MAX( BNORM, SAFMIN )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Left eigenvectors
</span><span class="comment">*</span><span class="comment">
</span>      IF( COMPL ) THEN
         IEIG = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Main loop over eigenvalues
</span><span class="comment">*</span><span class="comment">
</span>         DO 140 JE = 1, N
            IF( ILALL ) THEN
               ILCOMP = .TRUE.
            ELSE
               ILCOMP = SELECT( JE )
            END IF
            IF( ILCOMP ) THEN
               IEIG = IEIG + 1
<span class="comment">*</span><span class="comment">
</span>               IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
     $             ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Singular matrix pencil -- return unit eigenvector
</span><span class="comment">*</span><span class="comment">
</span>                  DO 50 JR = 1, N
                     VL( JR, IEIG ) = CZERO
   50             CONTINUE
                  VL( IEIG, IEIG ) = CONE
                  GO TO 140
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Non-singular eigenvalue:
</span><span class="comment">*</span><span class="comment">              Compute coefficients  a  and  b  in
</span><span class="comment">*</span><span class="comment">                   H
</span><span class="comment">*</span><span class="comment">                 y  ( a A - b B ) = 0
</span><span class="comment">*</span><span class="comment">
</span>               TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
     $                ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
               SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
               SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
               ACOEFF = SBETA*ASCALE
               BCOEFF = SALPHA*BSCALE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Scale to avoid underflow
</span><span class="comment">*</span><span class="comment">
</span>               LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
               LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
     $               SMALL
<span class="comment">*</span><span class="comment">
</span>               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
               END IF
<span class="comment">*</span><span class="comment">
</span>               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 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                                              H
</span><span class="comment">*</span><span class="comment">              Triangular solve of  (a A - b B)  y = 0
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                                      H
</span><span class="comment">*</span><span class="comment">              (rowwise in  (a A - b B) , or columnwise in a A - b B)
</span><span class="comment">*</span><span class="comment">
</span>               DO 100 J = JE + 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Compute
</span><span class="comment">*</span><span class="comment">                       j-1
</span><span class="comment">*</span><span class="comment">                 SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
</span><span class="comment">*</span><span class="comment">                       k=je
</span><span class="comment">*</span><span class="comment">                 (Scale if necessary)
</span><span class="comment">*</span><span class="comment">
</span>                  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
                  END IF
                  SUMA = CZERO
                  SUMB = CZERO
<span class="comment">*</span><span class="comment">
</span>                  DO 80 JR = JE, J - 1
                     SUMA = SUMA + DCONJG( S( JR, J ) )*WORK( JR )
                     SUMB = SUMB + DCONJG( P( JR, J ) )*WORK( JR )
   80             CONTINUE
                  SUM = ACOEFF*SUMA - DCONJG( BCOEFF )*SUMB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 with scaling and perturbation of the denominator
</span><span class="comment">*</span><span class="comment">
</span>                  D = DCONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
                  IF( ABS1( D ).LE.DMIN )
     $               D = DCMPLX( DMIN )
<span class="comment">*</span><span class="comment">
</span>                  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
                     END IF