stgsna.f.html

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

      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0,
     $                   FOUR = 4.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS
      INTEGER            I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2
      REAL               ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND,
     $                   EPS, LNRM, RNRM, ROOT1, ROOT2, SCALE, SMLNUM,
     $                   TMPII, TMPIR, TMPRI, TMPRR, UHAV, UHAVI, UHBV,
     $                   UHBVI
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      REAL               DUMMY( 1 ), DUMMY1( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.301"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      REAL               SDOT, <a name="SLAMCH.302"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLAPY2.302"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>, SNRM2
      EXTERNAL           <a name="LSAME.303"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, SDOT, <a name="SLAMCH.303"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLAPY2.303"></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           SGEMV, <a name="SLACPY.306"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>, <a name="SLAG2.306"></a><a href="slag2.f.html#SLAG2.1">SLAG2</a>, <a name="STGEXC.306"></a><a href="stgexc.f.html#STGEXC.1">STGEXC</a>, <a name="STGSYL.306"></a><a href="stgsyl.f.html#STGSYL.1">STGSYL</a>, <a name="XERBLA.306"></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          MAX, MIN, 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.315"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> )
      WANTS = <a name="LSAME.316"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) .OR. WANTBH
      WANTDF = <a name="LSAME.317"></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.319"></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
      LQUERY = ( LWORK.EQ.-1 )
<span class="comment">*</span><span class="comment">
</span>      IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
         INFO = -1
      ELSE IF( .NOT.<a name="LSAME.326"></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( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -8
      ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
         INFO = -10
      ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
         INFO = -12
      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( A( 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( N.EQ.0 ) THEN
            LWMIN = 1
         ELSE IF( <a name="LSAME.371"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'V'</span> ) .OR. <a name="LSAME.371"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> ) ) THEN
            LWMIN = 2*N*( N + 2 ) + 16
         ELSE
            LWMIN = N
         END IF
         WORK( 1 ) = LWMIN
<span class="comment">*</span><span class="comment">
</span>         IF( MM.LT.M ) THEN
            INFO = -15
         ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -18
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.386"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="STGSNA.386"></a><a href="stgsna.f.html#STGSNA.1">STGSNA</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         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">     Get machine constants
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.399"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = <a name="SLAMCH.400"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> ) / EPS
      KS = 0
      PAIR = .FALSE.
<span class="comment">*</span><span class="comment">
</span>      DO 20 K = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine whether A(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 20
         ELSE
            IF( K.LT.N )
     $         PAIR = A( 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 20
            ELSE
               IF( .NOT.SELECT( K ) )
     $            GO TO 20
            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( PAIR ) THEN
<span class="comment">*</span><span class="comment">