strsna.f.html

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

</span><span class="comment">*</span><span class="comment">  The reciprocal of the condition number of the right eigenvector u
</span><span class="comment">*</span><span class="comment">  corresponding to lambda is defined as follows. Suppose
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              T = ( lambda  c  )
</span><span class="comment">*</span><span class="comment">                  (   0    T22 )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Then the reciprocal condition number is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where sigma-min denotes the smallest singular value. We approximate
</span><span class="comment">*</span><span class="comment">  the smallest singular value by the reciprocal of an estimate of the
</span><span class="comment">*</span><span class="comment">  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
</span><span class="comment">*</span><span class="comment">  defined to be abs(T(1,1)).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  An approximate error bound for a computed right eigenvector VR(i)
</span><span class="comment">*</span><span class="comment">  is given by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                      EPS * norm(T) / SEP(i)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            PAIR, SOMCON, WANTBH, WANTS, WANTSP
      INTEGER            I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
      REAL               BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
     $                   MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            ISAVE( 3 )
      REAL               DUMMY( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.193"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      REAL               SDOT, <a name="SLAMCH.194"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLAPY2.194"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>, SNRM2
      EXTERNAL           <a name="LSAME.195"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, SDOT, <a name="SLAMCH.195"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLAPY2.195"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>, SNRM2
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="SLABAD.198"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>, <a name="SLACN2.198"></a><a href="slacn2.f.html#SLACN2.1">SLACN2</a>, <a name="SLACPY.198"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>, <a name="SLAQTR.198"></a><a href="slaqtr.f.html#SLAQTR.1">SLAQTR</a>, <a name="STREXC.198"></a><a href="strexc.f.html#STREXC.1">STREXC</a>, <a name="XERBLA.198"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode and test the input parameters
</span><span class="comment">*</span><span class="comment">
</span>      WANTBH = <a name="LSAME.207"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> )
      WANTS = <a name="LSAME.208"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) .OR. WANTBH
      WANTSP = <a name="LSAME.209"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'V'</span> ) .OR. WANTBH
<span class="comment">*</span><span class="comment">
</span>      SOMCON = <a name="LSAME.211"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( HOWMNY, <span class="string">'S'</span> )
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
         INFO = -1
      ELSE IF( .NOT.<a name="LSAME.216"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( HOWMNY, <span class="string">'A'</span> ) .AND. .NOT.SOMCON ) 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. ( WANTS .AND. LDVL.LT.N ) ) THEN
         INFO = -8
      ELSE IF( LDVR.LT.1 .OR. ( WANTS .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 eigenpairs for which condition numbers
</span><span class="comment">*</span><span class="comment">        are required, and test MM.
</span><span class="comment">*</span><span class="comment">
</span>         IF( SOMCON ) THEN
            M = 0
            PAIR = .FALSE.
            DO 10 K = 1, N
               IF( PAIR ) THEN
                  PAIR = .FALSE.
               ELSE
                  IF( K.LT.N ) THEN
                     IF( T( K+1, K ).EQ.ZERO ) THEN
                        IF( SELECT( K ) )
     $                     M = M + 1
                     ELSE
                        PAIR = .TRUE.
                        IF( SELECT( K ) .OR. SELECT( K+1 ) )
     $                     M = M + 2
                     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 = -13
         ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
            INFO = -16
         END IF
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.264"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="STRSNA.264"></a><a href="strsna.f.html#STRSNA.1">STRSNA</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>      IF( N.EQ.1 ) THEN
         IF( SOMCON ) THEN
            IF( .NOT.SELECT( 1 ) )
     $         RETURN
         END IF
         IF( WANTS )
     $      S( 1 ) = ONE
         IF( WANTSP )
     $      SEP( 1 ) = ABS( T( 1, 1 ) )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.287"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = <a name="SLAMCH.288"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> ) / EPS
      BIGNUM = ONE / SMLNUM
      CALL <a name="SLABAD.290"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SMLNUM, BIGNUM )
<span class="comment">*</span><span class="comment">
</span>      KS = 0
      PAIR = .FALSE.
      DO 60 K = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
</span><span class="comment">*</span><span class="comment">
</span>         IF( PAIR ) THEN
            PAIR = .FALSE.
            GO TO 60
         ELSE
            IF( K.LT.N )
     $         PAIR = T( K+1, K ).NE.ZERO
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine whether condition numbers are required for the k-th
</span><span class="comment">*</span><span class="comment">        eigenpair.
</span><span class="comment">*</span><span class="comment">
</span>         IF( SOMCON ) THEN
            IF( PAIR ) THEN
               IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
     $            GO TO 60
            ELSE
               IF( .NOT.SELECT( K ) )
     $            GO TO 60
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         KS = KS + 1
<span class="comment">*</span><span class="comment">
</span>         IF( WANTS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute the reciprocal condition number of the k-th
</span><span class="comment">*</span><span class="comment">           eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span>            IF( .NOT.PAIR ) 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">