ztrevc.f.html

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

            IF( SELECT( J ) )
     $         M = M + 1
   10    CONTINUE
      ELSE
         M = N
      END IF
<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 IF( MM.LT.M ) THEN
         INFO = -11
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.212"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZTREVC.212"></a><a href="ztrevc.f.html#ZTREVC.1">ZTREVC</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.223"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
      OVFL = ONE / UNFL
      CALL <a name="DLABAD.225"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>( UNFL, OVFL )
      ULP = <a name="DLAMCH.226"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
      SMLNUM = UNFL*( N / ULP )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Store the diagonal elements of T in working array WORK.
</span><span class="comment">*</span><span class="comment">
</span>      DO 20 I = 1, N
         WORK( I+N ) = T( I, I )
   20 CONTINUE
<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>      RWORK( 1 ) = ZERO
      DO 30 J = 2, N
         RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 )
   30 CONTINUE
<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>         IS = M
         DO 80 KI = N, 1, -1
<span class="comment">*</span><span class="comment">
</span>            IF( SOMEV ) THEN
               IF( .NOT.SELECT( KI ) )
     $            GO TO 80
            END IF
            SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
<span class="comment">*</span><span class="comment">
</span>            WORK( 1 ) = CMONE
<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 40 K = 1, KI - 1
               WORK( K ) = -T( K, KI )
   40       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve the triangular system:
</span><span class="comment">*</span><span class="comment">              (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
</span><span class="comment">*</span><span class="comment">
</span>            DO 50 K = 1, KI - 1
               T( K, K ) = T( K, K ) - T( KI, KI )
               IF( CABS1( T( K, K ) ).LT.SMIN )
     $            T( K, K ) = SMIN
   50       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( KI.GT.1 ) THEN
               CALL <a name="ZLATRS.274"></a><a href="zlatrs.f.html#ZLATRS.1">ZLATRS</a>( <span class="string">'Upper'</span>, <span class="string">'No transpose'</span>, <span class="string">'Non-unit'</span>, <span class="string">'Y'</span>,
     $                      KI-1, T, LDT, WORK( 1 ), SCALE, RWORK,
     $                      INFO )
               WORK( KI ) = SCALE
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy the vector x or Q*x to VR and normalize.
</span><span class="comment">*</span><span class="comment">
</span>            IF( .NOT.OVER ) THEN
               CALL ZCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 )
<span class="comment">*</span><span class="comment">
</span>               II = IZAMAX( KI, VR( 1, IS ), 1 )
               REMAX = ONE / CABS1( VR( II, IS ) )
               CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 )
<span class="comment">*</span><span class="comment">
</span>               DO 60 K = KI + 1, N
                  VR( K, IS ) = CMZERO
   60          CONTINUE
            ELSE
               IF( KI.GT.1 )
     $            CALL ZGEMV( <span class="string">'N'</span>, N, KI-1, CMONE, VR, LDVR, WORK( 1 ),
     $                        1, DCMPLX( SCALE ), VR( 1, KI ), 1 )
<span class="comment">*</span><span class="comment">
</span>               II = IZAMAX( N, VR( 1, KI ), 1 )
               REMAX = ONE / CABS1( VR( II, KI ) )
               CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Set back the original diagonal elements of T.
</span><span class="comment">*</span><span class="comment">
</span>            DO 70 K = 1, KI - 1
               T( K, K ) = WORK( K+N )
   70       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IS = IS - 1
   80    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>         IS = 1
         DO 130 KI = 1, N
<span class="comment">*</span><span class="comment">
</span>            IF( SOMEV ) THEN
               IF( .NOT.SELECT( KI ) )
     $            GO TO 130
            END IF
            SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
<span class="comment">*</span><span class="comment">
</span>            WORK( N ) = CMONE
<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 90 K = KI + 1, N
               WORK( K ) = -DCONJG( T( KI, K ) )
   90       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve the triangular system:
</span><span class="comment">*</span><span class="comment">              (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK.
</span><span class="comment">*</span><span class="comment">
</span>            DO 100 K = KI + 1, N
               T( K, K ) = T( K, K ) - T( KI, KI )
               IF( CABS1( T( K, K ) ).LT.SMIN )
     $            T( K, K ) = SMIN
  100       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( KI.LT.N ) THEN
               CALL <a name="ZLATRS.343"></a><a href="zlatrs.f.html#ZLATRS.1">ZLATRS</a>( <span class="string">'Upper'</span>, <span class="string">'Conjugate transpose'</span>, <span class="string">'Non-unit'</span>,
     $                      <span class="string">'Y'</span>, N-KI, T( KI+1, KI+1 ), LDT,
     $                      WORK( KI+1 ), SCALE, RWORK, INFO )
               WORK( KI ) = SCALE
            END IF
<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 ZCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 )
<span class="comment">*</span><span class="comment">
</span>               II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
               REMAX = ONE / CABS1( VL( II, IS ) )
               CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
<span class="comment">*</span><span class="comment">
</span>               DO 110 K = 1, KI - 1
                  VL( K, IS ) = CMZERO
  110          CONTINUE
            ELSE
               IF( KI.LT.N )
     $            CALL ZGEMV( <span class="string">'N'</span>, N, N-KI, CMONE, VL( 1, KI+1 ), LDVL,
     $                        WORK( KI+1 ), 1, DCMPLX( SCALE ),
     $                        VL( 1, KI ), 1 )
<span class="comment">*</span><span class="comment">
</span>               II = IZAMAX( N, VL( 1, KI ), 1 )
               REMAX = ONE / CABS1( VL( II, KI ) )
               CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Set back the original diagonal elements of T.
</span><span class="comment">*</span><span class="comment">
</span>            DO 120 K = KI + 1, N
               T( K, K ) = WORK( K+N )
  120       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IS = IS + 1
  130    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZTREVC.384"></a><a href="ztrevc.f.html#ZTREVC.1">ZTREVC</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>