dtrevc.f.html

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

                  END IF
<span class="comment">*</span><span class="comment">
</span>                  EMAX = ZERO
                  DO 120 K = 1, N
                     EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
     $                      ABS( VR( K, KI ) ) )
  120             CONTINUE
                  REMAX = ONE / EMAX
                  CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
                  CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            IS = IS - 1
            IF( IP.NE.0 )
     $         IS = IS - 1
  130       CONTINUE
            IF( IP.EQ.1 )
     $         IP = 0
            IF( IP.EQ.-1 )
     $         IP = 1
  140    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( LEFTV ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute left eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span>         IP = 0
         IS = 1
         DO 260 KI = 1, N
<span class="comment">*</span><span class="comment">
</span>            IF( IP.EQ.-1 )
     $         GO TO 250
            IF( KI.EQ.N )
     $         GO TO 150
            IF( T( KI+1, KI ).EQ.ZERO )
     $         GO TO 150
            IP = 1
<span class="comment">*</span><span class="comment">
</span>  150       CONTINUE
            IF( SOMEV ) THEN
               IF( .NOT.SELECT( KI ) )
     $            GO TO 250
            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 left 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 160 K = KI + 1, N
                  WORK( K+N ) = -T( KI, K )
  160          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Solve the quasi-triangular system:
</span><span class="comment">*</span><span class="comment">                 (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK
</span><span class="comment">*</span><span class="comment">
</span>               VMAX = ONE
               VCRIT = BIGNUM
<span class="comment">*</span><span class="comment">
</span>               JNXT = KI + 1
               DO 170 J = KI + 1, N
                  IF( J.LT.JNXT )
     $               GO TO 170
                  J1 = J
                  J2 = J
                  JNXT = J + 1
                  IF( J.LT.N ) THEN
                     IF( T( J+1, J ).NE.ZERO ) THEN
                        J2 = 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><span class="comment">*</span><span class="comment">                    Scale if necessary to avoid overflow when forming
</span><span class="comment">*</span><span class="comment">                    the right-hand side.
</span><span class="comment">*</span><span class="comment">
</span>                     IF( WORK( J ).GT.VCRIT ) THEN
                        REC = ONE / VMAX
                        CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
                        VMAX = ONE
                        VCRIT = BIGNUM
                     END IF
<span class="comment">*</span><span class="comment">
</span>                     WORK( J+N ) = WORK( J+N ) -
     $                             DDOT( J-KI-1, T( KI+1, J ), 1,
     $                             WORK( KI+1+N ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Solve (T(J,J)-WR)'*X = WORK
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="DLALN2.692"></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 if necessary
</span><span class="comment">*</span><span class="comment">
</span>                     IF( SCALE.NE.ONE )
     $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
                     WORK( J+N ) = X( 1, 1 )
                     VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
                     VCRIT = BIGNUM / VMAX
<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><span class="comment">*</span><span class="comment">                    Scale if necessary to avoid overflow when forming
</span><span class="comment">*</span><span class="comment">                    the right-hand side.
</span><span class="comment">*</span><span class="comment">
</span>                     BETA = MAX( WORK( J ), WORK( J+1 ) )
                     IF( BETA.GT.VCRIT ) THEN
                        REC = ONE / VMAX
                        CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
                        VMAX = ONE
                        VCRIT = BIGNUM
                     END IF
<span class="comment">*</span><span class="comment">
</span>                     WORK( J+N ) = WORK( J+N ) -
     $                             DDOT( J-KI-1, T( KI+1, J ), 1,
     $                             WORK( KI+1+N ), 1 )
<span class="comment">*</span><span class="comment">
</span>                     WORK( J+1+N ) = WORK( J+1+N ) -
     $                               DDOT( J-KI-1, T( KI+1, J+1 ), 1,
     $                               WORK( KI+1+N ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Solve
</span><span class="comment">*</span><span class="comment">                      [T(J,J)-WR   T(J,J+1)     ]'* X = SCALE*( WORK1 )
</span><span class="comment">*</span><span class="comment">                      [T(J+1,J)    T(J+1,J+1)-WR]             ( WORK2 )
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="DLALN2.731"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a>( .TRUE., 2, 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 if necessary
</span><span class="comment">*</span><span class="comment">
</span>                     IF( SCALE.NE.ONE )
     $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
                     WORK( J+N ) = X( 1, 1 )
                     WORK( J+1+N ) = X( 2, 1 )
<span class="comment">*</span><span class="comment">
</span>                     VMAX = MAX( ABS( WORK( J+N ) ),
     $                      ABS( WORK( J+1+N ) ), VMAX )
                     VCRIT = BIGNUM / VMAX
<span class="comment">*</span><span class="comment">
</span>                  END IF
  170          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy the vector x or Q*x to VL and normalize.
</span><span class="comment">*</span><span class="comment">
</span>               IF( .NOT.OVER ) THEN
                  CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
                  REMAX = ONE / ABS( VL( II, IS ) )
                  CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  DO 180 K = 1, KI - 1
                     VL( K, IS ) = ZERO
  180             CONTINUE
<span class="comment">*</span><span class="comment">
</span>               ELSE
<span class="comment">*</span><span class="comment">
</span>                  IF( KI.LT.N )
     $               CALL DGEMV( <span class="string">'N'</span>, N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
     $                           WORK( KI+1+N ), 1, WORK( KI+N ),
     $                           VL( 1, KI ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  II = IDAMAX( N, VL( 1, KI ), 1 )
                  REMAX = ONE / ABS( VL( II, KI ) )
                  CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
<span class="comment">*</span><span class="comment">
</span>               END IF
<span class="comment">*</span><span class="comment">
</span>            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Complex left eigenvector.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">               Initial solve:
</span><span class="comment">*</span><span class="comment">                 ((T(KI,KI)    T(KI,KI+1) )' - (WR - I* WI))*X = 0.
</span><span class="comment">*</span><span class="comment">                 ((T(KI+1,KI) T(KI+1,KI+1))                )
</span><span class="comment">*</span><span class="comment">
</span>               IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
                  WORK( KI+N ) = WI / T( KI, KI+1 )
                  WORK( KI+1+N2 ) = ONE