dtrevc.f.html

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

<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
         INFO = -1
      ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
         INFO = -8
      ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
         INFO = -10
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set M to the number of columns required to store the selected
</span><span class="comment">*</span><span class="comment">        eigenvectors, standardize the array SELECT if necessary, and
</span><span class="comment">*</span><span class="comment">        test MM.
</span><span class="comment">*</span><span class="comment">
</span>         IF( SOMEV ) THEN
            M = 0
            PAIR = .FALSE.
            DO 10 J = 1, N
               IF( PAIR ) THEN
                  PAIR = .FALSE.
                  SELECT( J ) = .FALSE.
               ELSE
                  IF( J.LT.N ) THEN
                     IF( T( J+1, J ).EQ.ZERO ) THEN
                        IF( SELECT( J ) )
     $                     M = M + 1
                     ELSE
                        PAIR = .TRUE.
                        IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
                           SELECT( J ) = .TRUE.
                           M = M + 2
                        END IF
                     END IF
                  ELSE
                     IF( SELECT( N ) )
     $                  M = M + 1
                  END IF
               END IF
   10       CONTINUE
         ELSE
            M = N
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( MM.LT.M ) THEN
            INFO = -11
         END IF
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.237"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DTREVC.237"></a><a href="dtrevc.f.html#DTREVC.1">DTREVC</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>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set the constants to control overflow.
</span><span class="comment">*</span><span class="comment">
</span>      UNFL = <a name="DLAMCH.248"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
      OVFL = ONE / UNFL
      CALL <a name="DLABAD.250"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>( UNFL, OVFL )
      ULP = <a name="DLAMCH.251"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
      SMLNUM = UNFL*( N / ULP )
      BIGNUM = ( ONE-ULP ) / SMLNUM
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute 1-norm of each column of strictly upper triangular
</span><span class="comment">*</span><span class="comment">     part of T to control overflow in triangular solver.
</span><span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = ZERO
      DO 30 J = 2, N
         WORK( J ) = ZERO
         DO 20 I = 1, J - 1
            WORK( J ) = WORK( J ) + ABS( T( I, J ) )
   20    CONTINUE
   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Index IP is used to specify the real or complex eigenvalue:
</span><span class="comment">*</span><span class="comment">       IP = 0, real eigenvalue,
</span><span class="comment">*</span><span class="comment">            1, first of conjugate complex pair: (wr,wi)
</span><span class="comment">*</span><span class="comment">           -1, second of conjugate complex pair: (wr,wi)
</span><span class="comment">*</span><span class="comment">
</span>      N2 = 2*N
<span class="comment">*</span><span class="comment">
</span>      IF( RIGHTV ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute right eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span>         IP = 0
         IS = M
         DO 140 KI = N, 1, -1
<span class="comment">*</span><span class="comment">
</span>            IF( IP.EQ.1 )
     $         GO TO 130
            IF( KI.EQ.1 )
     $         GO TO 40
            IF( T( KI, KI-1 ).EQ.ZERO )
     $         GO TO 40
            IP = -1
<span class="comment">*</span><span class="comment">
</span>   40       CONTINUE
            IF( SOMEV ) THEN
               IF( IP.EQ.0 ) THEN
                  IF( .NOT.SELECT( KI ) )
     $               GO TO 130
               ELSE
                  IF( .NOT.SELECT( KI-1 ) )
     $               GO TO 130
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute the KI-th eigenvalue (WR,WI).
</span><span class="comment">*</span><span class="comment">
</span>            WR = T( KI, KI )
            WI = ZERO
            IF( IP.NE.0 )
     $         WI = SQRT( ABS( T( KI, KI-1 ) ) )*
     $              SQRT( ABS( T( KI-1, KI ) ) )
            SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
<span class="comment">*</span><span class="comment">
</span>            IF( IP.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Real right eigenvector
</span><span class="comment">*</span><span class="comment">
</span>               WORK( KI+N ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Form right-hand side
</span><span class="comment">*</span><span class="comment">
</span>               DO 50 K = 1, KI - 1
                  WORK( K+N ) = -T( K, KI )
   50          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Solve the upper quasi-triangular system:
</span><span class="comment">*</span><span class="comment">                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
</span><span class="comment">*</span><span class="comment">
</span>               JNXT = KI - 1
               DO 60 J = KI - 1, 1, -1
                  IF( J.GT.JNXT )
     $               GO TO 60
                  J1 = J
                  J2 = J
                  JNXT = J - 1
                  IF( J.GT.1 ) THEN
                     IF( T( J, J-1 ).NE.ZERO ) THEN
                        J1 = J - 1
                        JNXT = J - 2
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( J1.EQ.J2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    1-by-1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="DLALN2.342"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a>( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
     $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
     $                            ZERO, X, 2, SCALE, XNORM, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Scale X(1,1) to avoid overflow when updating
</span><span class="comment">*</span><span class="comment">                    the right-hand side.
</span><span class="comment">*</span><span class="comment">
</span>                     IF( XNORM.GT.ONE ) THEN
                        IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
                           X( 1, 1 ) = X( 1, 1 ) / XNORM
                           SCALE = SCALE / XNORM
                        END IF
                     END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Scale if necessary
</span><span class="comment">*</span><span class="comment">
</span>                     IF( SCALE.NE.ONE )
     $                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
                     WORK( J+N ) = X( 1, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Update right-hand side
</span><span class="comment">*</span><span class="comment">
</span>                     CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
     $                           WORK( 1+N ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    2-by-2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="DLALN2.371"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a>( .FALSE., 2, 1, SMIN, ONE,
     $                            T( J-1, J-1 ), LDT, ONE, ONE,
     $                            WORK( J-1+N ), N, WR, ZERO, X, 2,
     $                            SCALE, XNORM, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Scale X(1,1) and X(2,1) to avoid overflow when
</span><span class="comment">*</span><span class="comment">                    updating the right-hand side.
</span><span class="comment">*</span><span class="comment">
</span>                     IF( XNORM.GT.ONE ) THEN
                        BETA = MAX( WORK( J-1 ), WORK( J ) )
                        IF( BETA.GT.BIGNUM / XNORM ) THEN
                           X( 1, 1 ) = X( 1, 1 ) / XNORM
                           X( 2, 1 ) = X( 2, 1 ) / XNORM