dtgevc.f.html

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

     $             ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Singular matrix pencil -- unit eigenvector
</span><span class="comment">*</span><span class="comment">
</span>                  IEIG = IEIG - 1
                  DO 230 JR = 1, N
                     VR( JR, IEIG ) = ZERO
  230             CONTINUE
                  VR( IEIG, IEIG ) = ONE
                  GO TO 500
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Clear vector
</span><span class="comment">*</span><span class="comment">
</span>            DO 250 JW = 0, NW - 1
               DO 240 JR = 1, N
                  WORK( ( JW+2 )*N+JR ) = ZERO
  240          CONTINUE
  250       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute coefficients in  ( a A - b B ) x = 0
</span><span class="comment">*</span><span class="comment">              a  is  ACOEF
</span><span class="comment">*</span><span class="comment">              b  is  BCOEFR + i*BCOEFI
</span><span class="comment">*</span><span class="comment">
</span>            IF( .NOT.ILCPLX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Real eigenvalue
</span><span class="comment">*</span><span class="comment">
</span>               TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
     $                ABS( P( JE, JE ) )*BSCALE, SAFMIN )
               SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
               SBETA = ( TEMP*P( JE, JE ) )*BSCALE
               ACOEF = SBETA*ASCALE
               BCOEFR = SALFAR*BSCALE
               BCOEFI = ZERO
<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>               SCALE = ONE
               LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
               LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
     $               SMALL
               IF( LSA )
     $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
               IF( LSB )
     $            SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
     $                    MIN( BNORM, BIG ) )
               IF( LSA .OR. LSB ) THEN
                  SCALE = MIN( SCALE, ONE /
     $                    ( SAFMIN*MAX( ONE, ABS( ACOEF ),
     $                    ABS( BCOEFR ) ) ) )
                  IF( LSA ) THEN
                     ACOEF = ASCALE*( SCALE*SBETA )
                  ELSE
                     ACOEF = SCALE*ACOEF
                  END IF
                  IF( LSB ) THEN
                     BCOEFR = BSCALE*( SCALE*SALFAR )
                  ELSE
                     BCOEFR = SCALE*BCOEFR
                  END IF
               END IF
               ACOEFA = ABS( ACOEF )
               BCOEFA = ABS( BCOEFR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              First component is 1
</span><span class="comment">*</span><span class="comment">
</span>               WORK( 2*N+JE ) = ONE
               XMAX = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute contribution from column JE of A and B to sum
</span><span class="comment">*</span><span class="comment">              (See &quot;Further Details&quot;, above.)
</span><span class="comment">*</span><span class="comment">
</span>               DO 260 JR = 1, JE - 1
                  WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) -
     $                             ACOEF*S( JR, JE )
  260          CONTINUE
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Complex eigenvalue
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLAG2.896"></a><a href="dlag2.f.html#DLAG2.1">DLAG2</a>( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP,
     $                     SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
     $                     BCOEFI )
               IF( BCOEFI.EQ.ZERO ) THEN
                  INFO = JE - 1
                  RETURN
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Scale to avoid over/underflow
</span><span class="comment">*</span><span class="comment">
</span>               ACOEFA = ABS( ACOEF )
               BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
               SCALE = ONE
               IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
     $            SCALE = ( SAFMIN / ULP ) / ACOEFA
               IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
     $            SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
               IF( SAFMIN*ACOEFA.GT.ASCALE )
     $            SCALE = ASCALE / ( SAFMIN*ACOEFA )
               IF( SAFMIN*BCOEFA.GT.BSCALE )
     $            SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
               IF( SCALE.NE.ONE ) THEN
                  ACOEF = SCALE*ACOEF
                  ACOEFA = ABS( ACOEF )
                  BCOEFR = SCALE*BCOEFR
                  BCOEFI = SCALE*BCOEFI
                  BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute first two components of eigenvector
</span><span class="comment">*</span><span class="comment">              and contribution to sums
</span><span class="comment">*</span><span class="comment">
</span>               TEMP = ACOEF*S( JE, JE-1 )
               TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
               TEMP2I = -BCOEFI*P( JE, JE )
               IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
                  WORK( 2*N+JE ) = ONE
                  WORK( 3*N+JE ) = ZERO
                  WORK( 2*N+JE-1 ) = -TEMP2R / TEMP
                  WORK( 3*N+JE-1 ) = -TEMP2I / TEMP
               ELSE
                  WORK( 2*N+JE-1 ) = ONE
                  WORK( 3*N+JE-1 ) = ZERO
                  TEMP = ACOEF*S( JE-1, JE )
                  WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF*
     $                             S( JE-1, JE-1 ) ) / TEMP
                  WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP
               END IF
<span class="comment">*</span><span class="comment">
</span>               XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
     $                ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute contribution from columns JE and JE-1
</span><span class="comment">*</span><span class="comment">              of A and B to the sums.
</span><span class="comment">*</span><span class="comment">
</span>               CREALA = ACOEF*WORK( 2*N+JE-1 )
               CIMAGA = ACOEF*WORK( 3*N+JE-1 )
               CREALB = BCOEFR*WORK( 2*N+JE-1 ) -
     $                  BCOEFI*WORK( 3*N+JE-1 )
               CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) +
     $                  BCOEFR*WORK( 3*N+JE-1 )
               CRE2A = ACOEF*WORK( 2*N+JE )
               CIM2A = ACOEF*WORK( 3*N+JE )
               CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE )
               CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE )
               DO 270 JR = 1, JE - 2
                  WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) +
     $                             CREALB*P( JR, JE-1 ) -
     $                             CRE2A*S( JR, JE ) + CRE2B*P( JR, JE )
                  WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) +
     $                             CIMAGB*P( JR, JE-1 ) -
     $                             CIM2A*S( JR, JE ) + CIM2B*P( JR, JE )
  270          CONTINUE
            END IF
<span class="comment">*</span><span class="comment">
</span>            DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Columnwise triangular solve of  (a A - b B)  x = 0
</span><span class="comment">*</span><span class="comment">
</span>            IL2BY2 = .FALSE.
            DO 370 J = JE - NW, 1, -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              If a 2-by-2 block, is in position j-1:j, wait until
</span><span class="comment">*</span><span class="comment">              next iteration to process it (when it will be j:j+1)
</span><span class="comment">*</span><span class="comment">
</span>               IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN
                  IF( S( J, J-1 ).NE.ZERO ) THEN
                     IL2BY2 = .TRUE.
                     GO TO 370
                  END IF
               END IF
               BDIAG( 1 ) = P( J, J )
               IF( IL2BY2 ) THEN
                  NA = 2
                  BDIAG( 2 ) = P( J+1, J+1 )
               ELSE
                  NA = 1
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute x(j) (and x(j+1), if 2-by-2 block)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLALN2.997"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a>( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ),
     $                      LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ),
     $                      N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP,
     $                      IINFO )
               IF( SCALE.LT.ONE ) THEN
<span class="comment">*</span><span class="comment">
</span>                  DO 290 JW = 0, NW - 1
                     DO 280 JR = 1, JE
                        WORK( ( JW+2 )*N+JR ) = SCALE*
     $                     WORK( ( JW+2 )*N+JR )
  280                CONTINUE
  290             CONTINUE
               END IF
               XMAX = MAX( SCALE*XMAX, TEMP )
<span class="comment">*</span><span class="comment">
</span>               DO 310 JW = 1, NW
                  DO 300 JA = 1, NA
                     WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW )
  300             CONTINUE
  310          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
</sp